#!/usr/bin/env python

# LISA-specific ILE fork.
#
# This executable is intentionally separate from
# integrate_likelihood_extrinsic_batchmode to avoid conflicts with the normal
# ILE and the incoming calmarg-in-loop work.  It was seeded from Aasim Jan's
# LISA-RIFT lisa_rift_paper branch and should be treated as an experimental
# LISA test surface until the newer driver architecture has a stable hook.

# This code marginalizes over extrinsic parameters by precomputing terms, setting up a sampler (based on adaptive importance sampling) and creating a likelihood function which\
# is passed downstream to the sampler. 

# FOLLOWUP TEST (if known injection)
#    util_PrintInjection.py --inj mdc.xml.gz --event 0 --verbose
#   util_PrintInjection.py --inj maxpt_zero_noise.xml.gz.xml.gz --event 0
"""
Integrate the extrinsic parameters of the prefactored likelihood function.
"""
import functools
from optparse import OptionParser, OptionGroup
import sys
import numpy
import numpy as np
import os
print("###########################################################################################")
print("# Loading packages and samplers")
print("###########################################################################################")
try:
  import cupy
  xpy_default=cupy
  identity_convert = cupy.asnumpy
  identity_convert_togpu = cupy.asarray
  identity_convert_lnL= identity_convert # see later, will be assigned something different if cupy active AND AC  AND gpu
  junk_to_check_installed = cupy.array(5)  # this will fail if GPU not installed correctly
  if not('RIFT_LOWLATENCY' in os.environ):
    print(cupy.show_config())  # print provenance/debugging information

    # Check memory allocation
    mem_info = cupy.cuda.Device().mem_info  # memory available in bytes
    n_mb_total = mem_info[1]/1024./1024.
    n_mb_free = mem_info[0]/1024/1024.
    print(  " cupy memory [total, available] {} {} Mb".format(n_mb_total,n_mb_free))
    # Add failure mode test if memory insufficient!
    if n_mb_free < 3000:
      print( "  - low cupy memory - ")
  cupy_success=True
except:
  print( ' no cupy')
#  import numpy as cupy  # will automatically replace cupy calls with numpy!
  xpy_default=numpy  # just in case, to make replacement clear and to enable override
  identity_convert = lambda x: x  # trivial return itself
  identity_convert_togpu = lambda x: x
  identity_convert_lnL= lambda x:x # see later
  cupy_success=False

import lal
from igwn_ligolw import utils, lsctables, table, ligolw
from igwn_ligolw.utils import process
import glue.lal

import RIFT.lalsimutils as lalsimutils
import RIFT.LISA.lalsimutils_compat as lisa_lalsimutils_compat
import RIFT.likelihood.factored_likelihood as factored_likelihood
import RIFT.likelihood.factored_likelihood_LISA as factored_likelihood_LISA
import RIFT.integrators.mcsampler as mcsampler
from RIFT.integrators.rvs_record import (RvsRecord as _RvsRecord,   # see DESIGN_rvs_naming.md
                                         SamplerOutputMixin)
import RIFT.misc.sky_rotations as sky_rotations
try:
    import RIFT.integrators.mcsamplerEnsemble as mcsamplerEnsemble
    mcsampler_gmm_ok = True
except:
    print(" No mcsamplerEnsemble ")
    mcsampler_gmm_ok = False
try:
    import RIFT.integrators.mcsamplerGPU as mcsamplerGPU
    mcsampler_gpu_ok = True
except:
    print( " No mcsamplerGPU ")
    mcsampler_gpu_ok = False
try:
    import RIFT.integrators.mcsamplerAdaptiveVolume as mcsamplerAdaptiveVolume
    mcsampler_AV_ok = True
except:
    print(" No mcsamplerAV ")
    mcsampler_AV_ok = False
mcsampler_Portfolio_ok=False
try:
    import RIFT.integrators.mcsamplerPortfolio as mcsamplerPortfolio
    mcsampler_Portfolio_ok = True
except:
    print(" No mcsamplerPortolfio ")
lisa_lalsimutils_compat.install_choose_waveform_print_params_lisa()

import RIFT.likelihood.priors_utils as priors_utils
import RIFT.misc.xmlutils as xmlutils

###########################################################################################
class EvenBivariateLinearInterpolator:
    def __init__(self, x0, dx, y0, dy, f):
        self._x0 = x0
        self._dx = dx
        self._y0 = y0
        self._dy = dy
        self._fgrid = xpy_default.asarray(f)

        self._dx_inv = 1.0 / self._dx
        self._dy_inv = 1.0 / self._dy

        self._N, self._M = self._fgrid.shape

    def __call__(self, x, y):
        # Compute the fractional indices into the lookup table where the free
        # parameters lie.
        i_mid = self._dx_inv * (x - self._x0)
        j_mid = self._dy_inv * (y - self._y0)

        # Compute the floor and ceiling of the fractional indices to get the
        # indices of the boundaries of the bin `x` and `y` lie in.
        # NOTE: In the case where `x` or `y` lie directly on a boundary, the
        # floor and ceiling will be equal, but the output is written in such a
        # way that we'd still get the correct result.
        i_lo = xpy_default.floor(i_mid).astype(int)
        j_lo = xpy_default.floor(j_mid).astype(int)
        i_hi = xpy_default.ceil(i_mid).astype(int)
        j_hi = xpy_default.ceil(j_mid).astype(int)

        # Compute just the fractional part of each index from the low and high
        # points.
        p = i_mid - i_lo
        q = j_mid - j_lo
        p_ = 1-p
        q_ = 1-q

        # Compute the interpolated result.
        f_approx =  p_*q_ * self._fgrid[i_lo,j_lo]
        f_approx += p *q_ * self._fgrid[i_hi,j_lo]
        f_approx += p_*q  * self._fgrid[i_lo,j_hi]
        f_approx += p *q  * self._fgrid[i_hi,j_hi]

        return f_approx


__author__ = "Evan Ochsner <evano@gravity.phys.uwm.edu>, Chris Pankow <pankow@gravity.phys.uwm.edu>, R. O'Shaughnessy"



#
# Pinnable parameters -- for command line processing
#
LIKELIHOOD_PINNABLE_PARAMS = ["right_ascension", "declination", "psi", "distance", "phi_orb", "t_ref", "inclination"]

def get_pinned_params(opts):
    """
    Retrieve a dictionary of user pinned parameters and their pin values.
    """
    return dict([(p,v) for p, v in opts.__dict__.items() if p in LIKELIHOOD_PINNABLE_PARAMS and v is not None]) 

def get_unpinned_params(opts, params):
    """
    Retrieve a set of unpinned parameters.
    """
    return params - set([p for p, v in opts.__dict__.items() if p in LIKELIHOOD_PINNABLE_PARAMS and v is not None])

def zero_like(*args,**kwargs):
  if len(kwargs)>0:
    arg0 = kwargs['psi']
  else:
    arg0 = args[0]
  return xpy_default.zeros(len(arg0))

def unit_like(*args,**kwargs):
  if len(kwargs)>0:
    arg0 = kwargs['psi']
  else:
    arg0 = args[0]
  return xpy_default.ones(len(arg0))


def sampler_param_tuple(sampler,args):
  param_list = sampler.params_ordered
  return tuple( [param_list.index(x) for x in args])


####################### Option parsing ###################################
optp = OptionParser()
optp.add_option( "--zero-likelihood", action='store_true', help="Run with exactly zero likelihood.  Prior test (e.g., alternative distance priors)")
optp.add_option("-c", "--cache-file", default=None, help="LIGO cache file containing all data needed.")
optp.add_option("-C", "--channel-name", action="append", help="instrument=channel-name, e.g. H1=FAKE-STRAIN. Can be given multiple times for different instruments.")
optp.add_option("-p", "--psd-file", action="append", help="instrument=psd-file, e.g. H1=H1_PSD.xml.gz. Can be given multiple times for different instruments.")
optp.add_option("-k", "--skymap-file", help="Use skymap stored in given FITS file.")
optp.add_option("-x", "--coinc-xml", help="gstlal_inspiral XML file containing coincidence information.")
optp.add_option("-I", "--sim-xml", help="XML file containing mass grid to be evaluated")
optp.add_option("--sim-grid", help="Hyperpipeline ASCII grid file (RIFT_HYPERPIPELINE_V1) of intrinsic points to be evaluated. Row index selected by --event. LISA-fork analogue of the main ILE --sim-grid; populates mass1/mass2/spin1z/spin2z and ecliptic_longitude/latitude from the selected row.")
optp.add_option("-E", "--event", default=0,type=int, help="Event number used for this run")
optp.add_option("--n-events-to-analyze", default=1,type=int, help="Number of events to analyze from this XML")
optp.add_option("--soft-fail-event-range",action='store_true',help='Soft failure (exit 0) if event ID is out of range. This happens in pipelines, if we have pre-built a DAG attempting to analyze more points than we really have')
optp.add_option("-f", "--reference-freq", type=float, default=100.0, help="Waveform reference frequency. Required, default is 100 Hz.")
optp.add_option("--fmin-template", dest='fmin_template', type=float, default=40, help="Waveform starting frequency.  Default is 40 Hz. Also equal to starting frequency for integration") 
optp.add_option("--fmin-template-correct-for-lmax",action='store_true',help="Modify amount of data selected, waveform starting frequency to account for l-max, to better insure all requested modes start within the targeted band")
optp.add_option("--fmin-ifo", action='append' , help="Minimum frequency for each IFO. Implemented by setting the PSD=0 below this cutoff. Use with care.") 
#optp.add_option("--nr-params",default=None, help="List of specific NR parameters and groups (and masses?) to use for the grid.")
#optp.add_option("--nr-index",type=int,default=-1,help="Index of specific NR simulation to use [integer]. Mass used: mtot= m1+m2")
optp.add_option('--nr-group', default=None,help="If using a *ssingle specific simulation* specified on the command line, provide it here")
optp.add_option('--nr-param', default=None,help="If using a *ssingle specific simulation* specified on the command line, provide it here")
optp.add_option("--nr-lookup",action='store_true', help=" Look up parameters from an NR catalog, instead of using the approximant specified")
optp.add_option("--nr-lookup-group",action='append', help="Restriction on 'group' for NR lookup")
optp.add_option("--nr-hybrid-use",action='store_true',help="Enable use of NR (or ROM!) hybrid, using --approx as the default approximant and with a frequency fmin")
optp.add_option("--nr-hybrid-method",default="taper_add",help="Hybridization method for NR (or ROM!).  Passed through to LALHybrid. pseudo_aligned_from22 will provide ad-hoc higher modes, if the early-time hybridization model only includes the 22 mode")
optp.add_option("--rom-group",default=None)
optp.add_option("--rom-param",default=None)
optp.add_option("--rom-use-basis",default=False,action='store_true',help="Use the ROM basis for inner products.")
optp.add_option("--rom-limit-basis-size-to",default=None,type=int)
optp.add_option("--rom-integrate-intrinsic",default=False,action='store_true',help='Integrate over intrinsic variables. REQUIRES rom_use_basis at present. ONLY integrates in mass ratio as present')
optp.add_option("--nr-perturbative-extraction",default=False,action='store_true')
optp.add_option("--nr-perturbative-extraction-full",default=False,action='store_true')
optp.add_option("--nr-use-provided-strain",default=False,action='store_true')
optp.add_option("--no-memory",default=False,action='store_true', help="At present, turns off m=0 modes. Use with EXTREME caution only if requested by model developer")
optp.add_option("--restricted-mode-list-file",default=None,help="A list of ALL modes to use in likelihood. Incredibly dangerous. Only use when comparing with models which provide restricted mode sets, or otherwise to isolate the effect of subsets of modes on the whole")
optp.add_option("--use-gwsignal",default=False,action='store_true',help='Use gwsignal. In this case the approx name is passed as a string to the lalsimulation.gwsignal interface')
optp.add_option("--use-gwsignal-lmax-nyquist",default=None,type=int,help='Passes lmax_nyquist integer to the gwsignal waveform interface')
optp.add_option("--use-external-EOB",default=False,action='store_true')
optp.add_option("--maximize-only",default=False, action='store_true',help="After integrating, attempts to find the single best fitting point")
optp.add_option("--dump-lnL-time-series",default=False, action='store_true',help="(requires --sim-xml) Dump lnL(t) at the injected parameters")
optp.add_option("-a", "--approximant", default="TaylorT4", help="Waveform family to use for templates. Any approximant implemented in LALSimulation is valid.")
optp.add_option("-A", "--amp-order", type=int, default=0, help="Include amplitude corrections in template waveforms up to this e.g. (e.g. 5 <==> 2.5PN), default is Newtonian order.")
optp.add_option("--l-max", type=int, default=2, help="Include all (l,m) modes with l less than or equal to this value.")
optp.add_option("-s", "--data-start-time", type=float, default=None, help="GPS start time of data segment. If given, must also give --data-end-time. If not given, sane start and end time will automatically be chosen.")
optp.add_option("-e", "--data-end-time", type=float, default=None, help="GPS end time of data segment. If given, must also give --data-start-time. If not given, sane start and end time will automatically be chosen.")
optp.add_option("--data-integration-window-half",default=75*1e-3,type=float,help="Only change this window size if you are an expert. The window for time integration is -/+ this quantity around the event time")
optp.add_option("-F", "--fmax", type=float, help="Upper frequency of signal integration. Default is use PSD's maximum frequency.")
optp.add_option("--srate",default=16384,type=float,help="Sampling rate. Change ONLY IF YOU ARE ABSOLUTELY SURE YOU KNOW WHAT YOU ARE DOING.")
optp.add_option("-t", "--event-time", type=float, help="GPS time of the event --- probably the end time. Required if --coinc-xml not given.")
optp.add_option("-i", "--inv-spec-trunc-time", type=float, default=8., help="Timescale of inverse spectrum truncation in seconds (Default is 8 - give 0 for no truncation)")
optp.add_option("-w", "--window-shape", type=float, default=0, help="Shape of Tukey window to apply to data (default is no windowing)")
optp.add_option("--psd-window-shape", type=float, default=0, help="Shape of Tukey window that *was* applied to the PSD being passed. If nonzero, we will rescale the PSD by the ratio of window shape results")
optp.add_option("-m", "--time-marginalization", action="store_true", help="Perform marginalization over time via direct numerical integration. Default is false.")
#optp.add_option("--n-fairdraw-extrinsic-samples",default=None,type=int,help="Extracts a concrete number of fair draw extrinsic samples, bounded above by n_eff")
optp.add_option("--resample-time-marginalization",action='store_true', help="If time-marginalizaiton is true (and should almost always be true), at the end export step use resampling. REQUIRES using fairdraw-extrinsic-output")
optp.add_option("--srate-resample-time-marginalization", type=int, default=None,
                help="Requested lattice rate for time-posterior resampling. The LISA driver uses this as its interpolation/draw grid; default keeps its historical 1 ms grid.")
optp.add_option("--time-posterior-export", type="choice", choices=["auto", "continuous", "grid"], default="auto",
                help="How --resample-time-marginalization exports time. The LISA likelihood cannot yet evaluate the selected stencil at arbitrary times, so auto/grid use its interpolation lattice and continuous is accepted as the same best-available lattice export for drop-in CLI compatibility (with a warning).")
optp.add_option("-d", "--distance-marginalization", action="store_true", help="Perform marginalization over distance via a look-up table. Default is false.")
optp.add_option("-l", "--distance-marginalization-lookup-table", default=None, help="Look-up table for distance marginalization.")
optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures.  (Combine with --gpu to enable GPU use, where available)")
optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc")
optp.add_option("--force-gpu-only", action="store_true", help="Hard fail if no GPU present (assessed by cupy not loading)")
optp.add_option("--force-xpy", action="store_true", help="Use the xpy code path.  Use with --vectorized --gpu to use the fallback CPU-based code path. Useful for debugging.")
optp.add_option("-o", "--output-file", help="Save result to this file.")
optp.add_option("-O", "--output-format", default='xml', help="[xml|hdf5]")
optp.add_option("-S", "--save-samples", action="store_true", help="Save sample points to output-file. Requires --output-file to be defined.")
optp.add_option("-L", "--save-deltalnL", type=float, default=float("Inf"), help="Threshold on deltalnL for points preserved in output file.  Requires --output-file to be defined")
optp.add_option("-P", "--save-P", type=float,default=0.1, help="Threshold on cumulative probability for points preserved in output file.  Requires --output-file to be defined")
optp.add_option("--internal-hard-fail-on-error",action='store_true',help='If true, fails with exit code 1 if any point is unsuccessful')
optp.add_option("--internal-soft-fail-on-cuda-error",action='store_true',help='If true, returns with exit code 0 on any CUDA error. Use with care (e.g., if many jobs failing)')
optp.add_option("--internal-make-empty-file-on-error",action='store_true',help='If true, failed points generate empty output file. Protects against OSG workflow problems')
optp.add_option("--internal-waveform-fd-L-frame",action='store_true',help='If true, passes extra_waveform_kwargs = {fd_L_frame=True} to lalsimutils hlmoft. Impacts outputs of ChooseFDWaveform calls only.')
optp.add_option("--internal-waveform-fd-no-condition",action='store_true',help='If true, adds extra_waveform_kwargs = {no_condition=True} to lalsimutils hlmoft. Impacts outputs of ChooseFDWaveform calls only. Provided to enable controlled tests of conditioning impact on PE')
optp.add_option("--internal-waveform-extra-lalsuite-args",type=str,default=None)
optp.add_option("--verbose",action='store_true')
optp.add_option("--save-eccentricity", action="store_true")

# 
# Add LISA options
#
lisa_params = OptionGroup(optp, "LISA arguments")
lisa_params.add_option("--LISA", action="store_true", help = "Let the code know that you are analysing LISA data. For now makes the code use analyze_event_LISA instead of analyze_event.")
lisa_params.add_option("--h5-frame", action="store_true", help = "LISA injection frames are in h5 format. Data is assumed to be in time domain when this option is used.")
lisa_params.add_option("--h5-frame-FD", action="store_true", help = "LISA injection frames are in h5 format. Data is assumed to be in frequency domain when this option is used.")
lisa_params.add_option("--lisa-reference-time", default=0.0, type=float, help = "Time in seconds at reference frequency, if reference frequency is not provided then reference frequency is the frequency at max(f^2 * A_22(f)).")
lisa_params.add_option("--lisa-reference-frequency", default=None, type=float, help = "Reference frequency in Hz, if reference frequency is not provided then reference frequency is the frequency at max(f^2 * A_22(f)).")
lisa_params.add_option("--modes", default=None, help = "If you need specific modes, set modes to an array eg --modes '[(2,2),(3,3)]")
lisa_params.add_option("--lisa-fixed-sky", default=False, help="Not varying skylocation")
lisa_params.add_option("--ecliptic-latitude", default=0, help="Value of ecliptic latitude (beta) if sky location is fixed")
lisa_params.add_option("--ecliptic-longitude", default=0, help="Value of ecliptic longitude (lambda) if sky location is fixed")
optp.add_option_group(lisa_params)

#
# Add the integration options
#
integration_params = OptionGroup(optp, "Integration Parameters", "Control the integration with these options.")
# Default is actually None, but that tells the integrator to go forever or until n_eff is hit.
integration_params.add_option("--n-max", type=int, help="Total number of samples points to draw. If this number is hit before n_eff, then the integration will terminate. Default is 'infinite'.",default=1e7)
integration_params.add_option("--n-eff", type=int, default=100, help="Total number of effective samples points to calculate before the integration will terminate. Default is 100")
integration_params.add_option("--fairdraw-extrinsic-output", action='store_true' , help="Output is fair draw, rather than being comprehensive")
integration_params.add_option("--n-chunk", type=int, help="Chunk'.",default=10000)
integration_params.add_option("--convergence-tests-on",default=False,action='store_true')
integration_params.add_option("--seed", type=int, help="Random seed to use. Default is to not seed the RNG.  Seeds every backend the samplers draw through (numpy, cupy, torch), so a seeded run is reproducible on GPU as well as CPU.")
integration_params.add_option("--no-adapt", action="store_true", help="Turn off adaptive sampling. Adaptive sampling is on by default.")
integration_params.add_option("--force-adapt-all", action="store_true", help="Force adaptive sampling for all parameters.")
integration_params.add_option("--force-reset-all", action="store_true", help="Force reset of sampling every iteration. (Recommended if AC and not using no-adapt-after-first)")
integration_params.add_option("--no-adapt-distance", action="store_true", help="Turn off adaptive sampling, just for distance. Adaptive sampling is on by default.")
integration_params.add_option("--no-adapt-after-first",action='store_true',help="Disables adaptation after first iteration with significant lnL")
integration_params.add_option("--adapt-weight-exponent", type=float, default=1.0, help="Exponent to use with weights (likelihood integrand) when doing adaptive sampling. Used in tandem with --adapt-floor-level to prevent overconvergence. Default is 1.0.")
integration_params.add_option("--adapt-floor-level", type=float, default=0.1, help="Floor to use with weights (likelihood integrand) when doing adaptive sampling. This is necessary to ensure the *sampling* prior is non zero during adaptive sampling and to prevent overconvergence. Default is 0.1 (no floor)")
integration_params.add_option("--adapt-adapt",action='store_true',help="Adapt the tempering exponent")
integration_params.add_option("--adapt-log",action='store_true',help="Use a logarithmic tempering exponent")
integration_params.add_option("--interpolate-time", default=False,help="If using time marginalization, compute using a continuously-interpolated array. (Default=false)")
integration_params.add_option("--d-prior",default='Euclidean' ,type=str,help="Distance prior for dL.  Options are dL^2 (Euclidean), uniform and 'pseudo_cosmo'  .")
integration_params.add_option("--d-max", default=10000,type=float,help="Maximum distance in volume integral. Used to SET THE PRIOR; changing this value changes the numerical answer.")
integration_params.add_option("--d-min", default=1,type=float,help="Minimum distance in volume integral. Used to SET THE PRIOR; changing this value changes the numerical answer.")
integration_params.add_option("--declination-cosine-sampler",action='store_true',help="If specified, the parameter used for declination is cos(dec), not dec")
integration_params.add_option("--inclination-cosine-sampler",action='store_true',help="If specified, the parameter used for inclination is cos(dec), not dec")
integration_params.add_option("--internal-rotate-phase", action='store_true',help="If specified, the integration sampler uses phase_p ==phi+psi and phase_m == phi-psi as sampling coordinates, both ranging from 0 to 4 pi.  The prior is twice as large.")
integration_params.add_option("--internal-sky-network-coordinates",action='store_true',help="If specified, perform integration in sky coordinates aligned with the first two IFOs provided")
integration_params.add_option("--internal-sky-network-coordinates-raw",action='store_true',help="If specified, does not attempt to organize IFO network sensibly, uses them AS PROVIDED IN ORDER.")
integration_params.add_option("--manual-logarithm-offset",type=float,default=0,help="Target value of logarithm lnL. Integrand is reduced by exp(-manual_logarithm_offset).  Important for high-SNR sources!   Should be set dynamically")
integration_params.add_option("--auto-logarithm-offset",action='store_true',help="Use the 'guess_snr' field returned in the precompute stage to change --manual-logarithm-offset for each event.")
integration_params.add_option("--internal-use-lnL",action='store_true',help="likelihood returns lnL, and integrator integrates lnL")
integration_params.add_option("--sampler-method",default="adaptive_cartesian_gpu",help="adaptive_cartesian|GMM|adaptive_cartesian_gpu")
integration_params.add_option("--sampler-portfolio",default=None,action='append',type=str,help="comma-separated strings, matching sampler methods other than portfolio")
integration_params.add_option("--sampler-portfolio-args",default=None, action='append', type=str, help='eval-able dictionaryo to be passed to that sampler')
# Portfolio freeze/allocation policy.  Pure pass-through to the shared portfolio sampler.
# Definitions copied verbatim from bin/integrate_likelihood_extrinsic_batchmode; pinned by
# test_lisa_sampler_plumbing.py.
integration_params.add_option("--portfolio-adaptive-alloc",action='store_true',default=False,help="Portfolio: ENABLE (opt-in) adaptive-probe draw allocation -- concentrate draws on the best per-chunk-n_ess member. Good on strongly-correlated targets; NOT recommended for AV-favorable high-SNR events (it starves the slow-contracting AV workhorse). Off by default (legacy n_ess reweighting).")
integration_params.add_option("--portfolio-alloc-exponent",default=None,type=float,help="Portfolio: adaptive allocation ~ member_quality^exponent. Higher concentrates harder on the winner. Sampler default 1.0.")
integration_params.add_option("--portfolio-freeze-wt",default=None,type=float,help="Portfolio: a member whose balance weight is below this stops updating its proposal (subject to grace/revive/VARAHA-exemption). Sampler default 0.05.")
integration_params.add_option("--portfolio-grace-iters",default=None,type=int,help="Portfolio: never freeze ANY member during the first N integration chunks (let slow starters contract). Sampler default 25.")
integration_params.add_option("--portfolio-probe-period",default=None,type=int,help="Portfolio: round-robin probe one member at a raised draw share every N chunks (breaks the under-observation trap). 0 disables probing. Sampler default 4.")
integration_params.add_option("--portfolio-quality-signal",default=None,type=str,help="Portfolio adaptive allocation: which per-member quality signal to rank members by. 'global' (default) = marginal gain in POOLED n_eff per sample (credits weight mass, debits weight variance); 'credit' = q_mix-native MIS credit assignment, sum_i [frac_m q_m/q_mix]_i * w_i per drawn sample (credits a member for COVERING where the integrand is, even if it drew few samples there); 'ness' = legacy per-member Kish n_ess (scale-invariant, misranks a slow-contracting AV -- see DESIGN_portfolio_freeze_policy.md).")
integration_params.add_option("--portfolio-revive-period",default=None,type=int,help="Portfolio: every N chunks, update even a frozen member one step so it can recover. 0 disables. Sampler default 8.")
integration_params.add_option("--portfolio-varaha-can-freeze",action='store_true',default=False,help="Portfolio: DISABLE the VARAHA freeze-exemption, so VARAHA/AV members obey the grace/revive/weight freeze schedule like other members. Use only if a VARAHA member is a known-bad fit and you want to save its selfish-draw eval cycles.")
integration_params.add_option("--portfolio-varaha-max-frac",default=None,type=float,help="Portfolio: CAP the combined DRAW fraction of VARAHA/AV members (0/unset = no cap).  Use WITH --portfolio-varaha-min-frac to constrain the VARAHA share to a BAND.  Rationale: a floor alone stops the mixture degenerating to peaked-member-only (which strips q_mix of its broad backstop, so a missed mode goes uncovered and lnZ is silently low while n_eff looks GOOD), but the share can then run away the OTHER way to ~1 and the mixture degenerates to VARAHA-only instead.  A band (e.g. 0.25/0.75) keeps q_mix genuinely mixed by construction.  Unbiased either way (balance heuristic), so it costs at most draws, never correctness.")
integration_params.add_option("--portfolio-varaha-min-frac",default=None,type=float,help="Portfolio: reserve this combined DRAW fraction for VARAHA/AV members (0/unset = off). never-freeze keeps a VARAHA member UPDATING, but both allocation rules score by per-chunk n_ess, which sits at ~1 during VARAHA's slow cumulative contraction -- so a member that looks instantly good can take nearly the whole budget (measured on S250114ax post-#33: GMM took ~0.84 and the portfolio collapsed to n_eff ~2 vs ~100 for standalone AV). Unbiased for any allocation (q_mix); trades efficiency only.")
integration_params.add_option("--portfolio-varaha-never-freeze",action='store_true',default=False,help="Portfolio: VARAHA/AV members always update every chunk past their breakpoint (freeze-exempt). This is the sampler default; the flag is here for explicitness/pipe pass-through.")
integration_params.add_option("--portfolio-weight-clip",default=None,type=float,help="Portfolio: OPT-IN truncated importance sampling applied to the PROPOSAL-FIT INPUT ONLY. Caps the weights fed to member.update_sampling_prior (the GMM covariance fit) at tau = C*sqrt(n)*mean(w) (0/unset = off; C~1 is the standard Ionides choice), so one enormous weight cannot make that fit degenerate. The estimator (ln Z, n_eff), the n_ess report, and the allocation signal all use the TRUE unclipped weights, so they stay exactly unbiased and undistorted. Do NOT clip the estimator (measured on S250114ax: n_eff=100 2x faster than AV but ln Z biased -11.5 nats) or the n_ess report (clipping inflates the clipped member's n_ess and starves the AV workhorse). The withheld tail mass is tracked and reported as a diagnostic. NOTE: if huge weights come from q_mix UNDERFLOW (watch for the warning) they are a numerical artifact, not tail mass.")
integration_params.add_option("--sampler-xpy",default=None,help="numpy|cupy  if the adaptive_cartesian_gpu sampler is active, use that.")
# MC-error replicas.  Copied verbatim from bin/integrate_likelihood_extrinsic_batchmode;
# pinned by test_lisa_mc_error_replicas.py.
optp.add_option("--mc-error-replicas",default=0,type=int, help="MC-error stabilization: when the reported lnL error is untrustworthy (see the trigger options below), re-run the extrinsic integration this many EXTRA times as cold replicas (adaptation reset, sample cache dropped, fresh RNG draws) and report lnL from the LINEAR mean of the replica integrals with sigma from the max of the propagated error and the between-replica scatter (t-distributed, K-1 dof).  The naive per-run sigma is computed from the SAME weights as the integral, so it is small exactly when the run silently missed the peak; only independent replicas can see that.  NEVER combine replicas by inverse-variance weighting -- that overweights the worst replica.  The posterior/fairdraw export POOLS the replicas (weights renormalized so each contributes Z_k/K; fairdraw blocks contribute equal within-block weights, since those samples already carry their weights once), so the exported samples represent the same mixture as the reported evidence.  Default 0 = off (production behavior unchanged).")
optp.add_option("--mc-error-sigma-trigger",default=0.4,type=float, help="Replicate (see --mc-error-replicas) when the reported sigma_lnZ exceeds this value.")
optp.add_option("--mc-error-ess-trigger",default=30.,type=float, help="Replicate when the Kish effective sample size (sum w)^2/sum w^2 of the run's weights falls below this value.")
optp.add_option("--mc-error-khat-trigger",default=0.7,type=float, help="Replicate when the Pareto k-hat weight-tail diagnostic exceeds this value (0.7 = the PSIS reliability threshold: above it the weight variance is effectively unresolved and the naive sigma is a lower bound).")
# AV live-volume state, per-axis bin allocation, and the collapse gate.  Copied verbatim
# from bin/integrate_likelihood_extrinsic_batchmode; pinned by test_lisa_av_state.py.
integration_params.add_option("--sampler-save-state",default=None,help="AV only: after integration, write the adapted live-volume state (.npz) for reuse by later instances/iterations. Point --sampler-load-state at the same file across a grid to warm-start each point from the previous one.")
integration_params.add_option("--sampler-load-state",default=None,help="AV only: load a saved live-volume state (.npz from --sampler-save-state) to warm-start this integration. Overrides --sampler-warmstart-samples.")
integration_params.add_option("--sampler-anisotropic-bins",action="store_true",help="AV only: give each extrinsic axis a DIFFERENT number of bins during contraction -- fine where the live points cluster tightly (phase/polarization/sky), coarse where they are broad (distance/inclination) -- instead of the default equal split.  Keeps the same total bin budget, so the estimator is unchanged; helps AV wrap a correlated/degenerate posterior more tightly.")
optp.add_option("--reject-collapsed-live-volume",action='store_true',default=False, help="DROP an event whose adaptive-volume live volume degenerated (see the [AV COLLAPSE] report) instead of exporting it: the integration is treated as a failure, so no likelihood row, XML or posterior samples are written for it.  Such a run's lnZ and samples describe a single mode of the integrand and are NOT a fair posterior draw, and nothing downstream can distinguish them from a converged export.  Default off, because dropping the event silently THINS the posterior in an SNR-dependent way -- that was the pre-fix behaviour, when this case crashed.  Left off, the event is exported but announces itself loudly and (with --mc-error-replicas>0) triggers replication.  Turn it on when a contaminated point is worse than a missing one.")
# L0 auto-rescue.  Ported from bin/integrate_likelihood_extrinsic_batchmode; defaults and help
# text kept IDENTICAL there and here on purpose -- see test_lisa_l0_rescue.py, which pins them.
integration_params.add_option("--sampler-warmstart-retry-neff",type=float,default=None,help="AV or portfolio (L0 auto-rescue): if a pass finishes below this n_eff (i.e. it stalled on a very sharp / high-amplitude peak), automatically re-run a second pass warm-started from THIS point's own highest-likelihood samples.  Same-problem reuse in the sense that the seed provably contains the peak the cold pass found -- but NOT that every mode is represented, so the warm pass can be biased low if the seed missed one.  The rescue still runs as before; its result is rejected in favour of the cold pass only on positive evidence of lost mass (see --sampler-l0-rescue-reject-dlnZ).  A portfolio is unaffected: its GMM member carries a defensive component.  Directly targets the high-SNR n_eff LOTTERY (a large fraction of independent runs collapse to n_eff~1 by contracting onto the wrong spot); the rescue re-seeds a collapsed run from the peak it did find.  Recommended for high-SNR events; e.g. 5.")
integration_params.add_option("--sampler-l0-rescue-reject-dlnZ", type=float, default=3.0, help="Evidence threshold (nats) for rejecting the L0 rescue's warm pass: reject when the full-support cold pass reports lnZ this much HIGHER, which would indicate the seed missed mass.  Larger = more permissive.  DEFAULT RAISED 0.5 -> 3.0 ON MEASUREMENT (see test/expensive_before_merging/integrators/L0_REJECT_DLNZ_MEASUREMENT.md): across 160 known-lnZ passes the gate caught 0 of 55 genuinely truncated warm passes at EVERY threshold, while at 0.5 it binned 25% of GOOD portfolio warm passes.  0.5 was therefore strictly dominated -- it bought no detection and cost one good pass in four.  3.0 keeps a safety net for a genuinely large discrepancy at ~0% false-positive rate.  This gate is NOT a working truncation detector; do not rely on it as one.")
integration_params.add_option("--sampler-l0-rescue-accept-truncated", action='store_true', default=False, help="Report the L0 rescue's warm pass even when it lands well below the full-support cold pass (see --sampler-l0-rescue-reject-dlnZ).  Default OFF: on that evidence the cold result is kept instead, since the warm pass is confined to the seeded peak and may be missing a mode.  The rescue itself still runs either way.")
integration_params.add_option("--sampler-l0-rescue-puff-scale", type='choice', choices=['fixed','auto'], default='auto', help="How wide to puff the L0 rescue's seed when it is rank-deficient in the adaptive dimensions.  'auto' (default) measures the posterior scale AND correlations from every finite lnL the collapsed pass already drew; 'fixed' uses --sampler-l0-rescue-puff-width-frac of each parameter's prior range, which is the historical behaviour and knows nothing about the posterior (which narrows as 1/rho).  'auto' falls back to 'fixed' when there are too few finite points to estimate a covariance.")
integration_params.add_option("--sampler-l0-rescue-puff-width-frac", type=float, default=0.005, help="Isotropic puff width for the L0 rescue's rank-deficient seed, as a fraction of each parameter's prior range.  Used by --sampler-l0-rescue-puff-scale fixed, and as the 'auto' fallback.  Default 0.005 = the historical hardcoded 1/200.")
integration_params.add_option("--sampler-l0-rescue-puff-factor", type=float, default=2.0, help="Multiply the L0 rescue's puff width by this factor.  Default 2 is the measured optimum on a known-lnZ 6-D target (mean lnZ error +0.08 nats, ESS 52); BOTH tails are wrong, so do not treat wide as free -- x0.5 truncates (-8.5 nats), x6 biases high (+3.0) and costs efficiency, x12 is a cold start in all but name and re-collapses (-30).")
# Also consumed by the rescue (it is the lnL window build_warm_seed keeps), which is why it
# lands in this pass rather than with the sequential warm start it is named for.
integration_params.add_option("--sampler-sequential-warmstart-deltalnL",type=float,default=15.0,help="Keep previous-point samples within this lnL of the max as the warm seed for the next point.  Default 15.")
integration_params.add_option("--supplementary-likelihood-factor-code", default=None,type=str,help="Import a module (in your pythonpath!) containing a supplementary factor for the likelihood.  Used to impose supplementary external priors of arbitrary complexity and external dependence (e.g., EM observations). EXPERTS-ONLY")
integration_params.add_option("--supplementary-likelihood-factor-function", default=None,type=str,help="With above option, specifies the specific function used as an external prior. EXPERTS ONLY")
integration_params.add_option("--supplementary-likelihood-factor-ini", default=None,type=str,help="With above option, specifies an ini file that is parsed (here) and passed to the preparation code, called when the module is first loaded, to configure the module. EXPERTS ONLY")
optp.add_option_group(integration_params)

#
# Add the intrinsic parameters
#
intrinsic_params = OptionGroup(optp, "Intrinsic Parameters", "Intrinsic parameters (e.g component mass) to use.")
intrinsic_params.add_option("--pin-distance-to-sim",action='store_true', help="Pin *distance* value to sim entry. Used to enable source frame reconstruction with NR.")
intrinsic_params.add_option("--mass1", type=float, help="Value of first component mass, in solar masses. Required if not providing coinc tables.")
intrinsic_params.add_option("--mass2", type=float, help="Value of second component mass, in solar masses. Required if not providing coinc tables.")
intrinsic_params.add_option("--spin1z", type=float, help="Value of first component spin (aligned with angular momentum), dimensionless.")
intrinsic_params.add_option("--spin2z", type=float, help="Value of second  component spin (aligned with angular momentum), dimensionless.")
intrinsic_params.add_option("--eff-lambda", type=float, help="Value of effective tidal parameter. Optional, ignored if not given.")
intrinsic_params.add_option("--deff-lambda", type=float, help="Value of second effective tidal parameter. Optional, ignored if not given")
intrinsic_params.add_option("--export-eos-index",action='store_true')
optp.add_option_group(intrinsic_params)


#
# Add options to integrate over intrinsic parameters.  Same conventions as util_ManualOverlapGrid.py.  
# Parameters have special names, and we adopt priors that use those names.
# NOTE: Only 'q' implemented
#
intrinsic_int_params = OptionGroup(optp, "Intrinsic integrated parameters", "Intrinsic parameters to integrate over. ONLY currently used with ROM version")
intrinsic_int_params.add_option("--parameter",action='append')
intrinsic_int_params.add_option("--parameter-range",action='append',type=str)
intrinsic_int_params.add_option("--adapt-intrinsic",action='store_true')
optp.add_option_group(intrinsic_int_params)

#
# Add the pinnable parameters
pinnable = OptionGroup(optp, "Pinnable Parameters", "Specifying these command line options will pin the value of that parameter to the specified value with a probability of unity.")
for pin_param in LIKELIHOOD_PINNABLE_PARAMS:
    option = "--" + pin_param.replace("_", "-")
    pinnable.add_option(option, type=float, help="Pin the value of %s." % pin_param)
optp.add_option_group(pinnable)

opts, args = optp.parse_args()
from RIFT.likelihood.time_posterior import (
    legacy_time_interpolation_enabled, resolve_time_posterior_export_mode)
opts.interpolate_time = legacy_time_interpolation_enabled(opts.interpolate_time)
opts._time_posterior_export = resolve_time_posterior_export_mode(
    opts.time_posterior_export, "cubic" if opts.interpolate_time else "nearest",
    continuous_available=False)
if (opts.resample_time_marginalization and
        opts._time_posterior_export == "continuous"):
    print("WARNING: LISA ILE cannot evaluate the selected time stencil at arbitrary "
          "times; treating --time-posterior-export continuous as its best-available "
          "interpolation-lattice export for executable-swap compatibility.")
    opts._time_posterior_export = "grid"
if (opts.srate_resample_time_marginalization is not None and
        opts.srate_resample_time_marginalization <= 0):
    optp.error("--srate-resample-time-marginalization must be positive")

#
# Failure modes
#
ok_lnL_methods = ['GMM', 'adaptive_cartesian', 'adaptive_cartesian_gpu','AV','portfolio']
if opts.internal_use_lnL and not(opts.sampler_method  in ok_lnL_methods ):
  print(" OPTION MISMATCH : --internal-use-lnL not compatible with", opts.sampler_method, " can only use ", ok_lnL_methods)
  sys.exit(99)
# if we are on a GPU and using a GPU-accelerated likelihood, don't send the likelihood data back from the GPU needlessly
if cupy_success and opts.gpu and opts.sampler_method == 'adaptive_cartesian_gpu':
  identity_convert_lnL  = lambda x:x


if opts.resample_time_marginalization:
  import scipy.special
if opts.resample_time_marginalization and not(opts.fairdraw_extrinsic_output):
  raise Exception(" Resampled time output requires --fairdraw-extrinsic-output ")

# Fairdraw is NOT YET IMPLEMENTED for these other integrators!
#if opts.fairdraw_extrinsic_output:
#  if opts.sampler_method == 'adaptive_cartesian_gpu':
#    raise Exception(" Fairdraw not available  for this sampler")
#  if opts.sampler_method == 'GMM':
#    raise Exception(" Fairdraw not available  for this sampler")


supplemental_ln_likelihood= None
supplemental_ln_likelhood_prep=None
supplemental_ln_likelhood_parsed_ini=None
# Supplemental likelihood factor. Must have identical call sequence to 'likelihood_function'. Called with identical raw inputs (including cosines/etc)
if opts.supplementary_likelihood_factor_code and opts.supplementary_likelihood_factor_function:
  print(" EXTERNAL SUPPLEMENTARY LIKELIHOOD FACTOR : {}.{} ".format(opts.supplementary_likelihood_factor_code,opts.supplementary_likelihood_factor_function))
  __import__(opts.supplementary_likelihood_factor_code)
  external_likelihood_module = sys.modules[opts.supplementary_likelihood_factor_code]
  supplemental_ln_likelihood = getattr(external_likelihood_module,opts.supplementary_likelihood_factor_function)
  name_prep = "prepare_"+opts.supplementary_likelihood_factor_function
  if hasattr(external_likelihood_module,name_prep):
    supplemental_ln_likelhood_prep=getattr(external_likelihood_module,name_prep)
    # Check for and load in ini file associated with external library
    if opts.supplementary_likelihood_factor_ini:
      import configparser as ConfigParser
      config = ConfigParser.ConfigParser()
      config.optionxform=str # force preserve case! 
      config.read(opts.supplementary_likelihood_factor_ini)
      supplemental_ln_likelhood_parsed_ini=config

if opts.distance_marginalization:
  lookup_table = np.load(opts.distance_marginalization_lookup_table)
  opts.maximize_only  = False

if opts.gpu is None:
  opts.gpu = False  # None can be treated differently from false?

if opts.gpu and opts.force_gpu_only and not cupy_success:
  print("GPU requested but no GPU/cupy found: Hard fail by request")
  sys.exit(35)  # unique code to make resuming due to this error easier to identify

if opts.gpu and xpy_default is numpy:
    print( " Override --gpu  (not available);  use --force-xpy to require the identical code path is used (with xpy =[np|cupy]")
    opts.gpu=False
    if opts.force_xpy:
        opts.gpu=True

manual_avoid_overflow_logarithm=opts.manual_logarithm_offset
manual_avoid_overflow_logarithm_default =  manual_avoid_overflow_logarithm

# LISA check (Sampling rate)
deltaT = None
fSample= opts.srate # change sampling rate
if not(fSample is None):
  deltaT =1./fSample


# Load in restricted mode set, if available
restricted_mode_list=None
if not(opts.restricted_mode_list_file is None):
    modes =numpy.loadtxt(opts.restricted_mode_list_file,dtype=int) # columns are l m.  Must contain all. Only integers obviously
    restricted_mode_list = [ (l,m) for l,m in modes]
    print( " RESTRICTED MODE LIST target :", restricted_mode_list)

intrinsic_param_names = opts.parameter
valid_intrinsic_param_names = ['q']
if intrinsic_param_names:
 for param in intrinsic_param_names:
    
    # Check if in the valid list
    if not(param in valid_intrinsic_param_names):
            print( ' Invalid param ', param, ' not in ', valid_intrinsic_param_names)
            sys.exit(1)
    param_ranges = []
    if len(intrinsic_param_names) == len(opts.parameter_range):
        param_ranges = numpy.array(map(eval, opts.parameter_range))
        # Rescale mass-dependent ranges to SI units
        for p in ['mc', 'm1', 'm2', 'mtot']:
          if p in intrinsic_param_names:
            indx = intrinsic_param_names.index(p)
            param_ranges[indx]= numpy.array(param_ranges[indx])* lal.MSUN_SI


# LISA check (what if I am working in FD?)
# Check both or neither of --data-start/end-time given
if opts.data_start_time is None and opts.data_end_time is not None:
    raise ValueError("You must provide both or neither of --data-start-time and --data-end-time.")
if opts.data_end_time is None and opts.data_start_time is not None:
    raise ValueError("You must provide both or neither of --data-start-time and --data-end-time.")

#
# Import NR grid
#
NR_template_group=None
NR_template_param=None
if opts.nr_group and opts.nr_param:
    import NRWaveformCatalogManager3 as nrwf
    NR_template_group = opts.nr_group
    if nrwf.internal_ParametersAreExpressions[NR_template_group]:
        NR_template_param = eval(opts.nr_param)
    else:
        NR_template_param = opts.nr_param


#
# Hardcoded variables
#
template_min_freq = opts.fmin_template # minimum frequency of template
#t_ref_wind = 50e-3 # Interpolate in a window +/- this width about event time. 
t_ref_wind = opts.data_integration_window_half
T_safety = 2. # Safety buffer (in sec) for wraparound corruption

#
# Inverse spectrum truncation control
#
T_spec = opts.inv_spec_trunc_time
if T_spec == 0.: # Do not do inverse spectrum truncation
    inv_spec_trunc_Q = False
    T_safety += 8. # Add a bit more safety buffer in this case
else:
    inv_spec_trunc_Q = True

#
# Integrator options
#
n_max = opts.n_max # Max number of extrinsic points to evaluate at
n_eff = opts.n_eff # Effective number of points evaluated

#
# Initialize the RNG, if needed
#
# Seed EVERY backend a sampler can draw from, not just numpy: the samplers draw
# through self.xpy / xpy_default, which is cupy on GPU, and cupy has its own
# global generator.  See RIFT/integrators/seeding.py.
if opts.seed is not None:
    from RIFT.integrators.seeding import seed_everything
    seed_everything(opts.seed)

# LISA check, reference time instead of event time
if not(opts.LISA):
  if opts.event_time is not None:
      event_time = glue.lal.LIGOTimeGPS(opts.event_time)
      print( "Event time from command line: %s" % str(event_time))
  else:
      print( " Error! ")
      sys.exit(1)
  fiducial_epoch = lal.LIGOTimeGPS()
  fiducial_epoch = event_time.seconds + 1e-9*event_time.nanoseconds   # no more direct access to gpsSeconds
else:
    fiducial_epoch = opts.lisa_reference_time
#
# Template descriptors
#


print("###########################################################################################")
print("# Reading grid, data and psd")
print("###########################################################################################")
# Struct to hold template parameters, P.deltaF comes from data (data.deltaF)
P_list = None
P=None # force allocation so I can use the preferred event later
if opts.sim_grid:
    # Hyperpipeline intrinsic-grid handoff (the DAG/CEPP path passes one row per
    # ILE worker via --sim-grid + --event).  The LISA fork is single-point per
    # job, so resolve the selected row here into the mass1/mass2/spin*/ecliptic
    # opts and fall through to the standard single-point P construction below.
    from RIFT.misc import hyperpipeline_io as _hpio
    print("====Loading injection grid file:", opts.sim_grid, " event", opts.event, "=======")
    _grid_arr, _grid_cols = _hpio.read_table(opts.sim_grid)
    if opts.event is None or opts.event >= len(_grid_arr):
        print(" Event index out of range for grid; soft exit")
        sys.exit(0)
    _row = _grid_arr[opts.event]
    def _grid_get(name, default=0.0):
        return float(_row[name]) if name in _grid_cols else default
    opts.mass1 = _grid_get('m1')
    opts.mass2 = _grid_get('m2')
    opts.spin1z = _grid_get('a1z')
    opts.spin2z = _grid_get('a2z')
    if 'ecliptic_longitude' in _grid_cols:
        opts.ecliptic_longitude = _grid_get('ecliptic_longitude')
    if 'ecliptic_latitude' in _grid_cols:
        opts.ecliptic_latitude = _grid_get('ecliptic_latitude')
    print("   grid row -> m1={} m2={} s1z={} s2z={} ecl_long={} ecl_lat={}".format(
        opts.mass1, opts.mass2, opts.spin1z, opts.spin2z,
        opts.ecliptic_longitude, opts.ecliptic_latitude))
if opts.sim_xml:
    print(f"====Loading injection XML: {opts.sim_xml}, reading from event {opts.event} to event {opts.event+opts.n_events_to_analyze} =======")
    P_list = lalsimutils.xml_to_ChooseWaveformParams_array(str(opts.sim_xml))
    if  len(P_list) < opts.event: 
#+opts.n_events_to_analyze:
        print( " Event list of range; soft exit")
        sys.exit(0)
    n_event_max= np.min([len(P_list), opts.event+opts.n_events_to_analyze])
    P_list = P_list[opts.event:n_event_max]
    if len(P_list) ==0:
      print(" No events to analyze, terminating with success ")
      sys.exit(0)
    for P in P_list:
        P.radec =False  # do NOT propagate the epoch later
        # LISA 
        if opts.LISA:
           #P.fref = opts.lisa_reference_frequency
           P.fref = opts.reference_freq # we are distinguishing between tempalte fref and lisa_fref
           P.tref = 0.0 
           P.dist = factored_likelihood_LISA.distGpcRef * 1.e9 * lal.PC_SI 
        else:
          P.fref = opts.reference_freq
          P.tref = fiducial_epoch  # the XML table
          P.dist = factored_likelihood.distMpcRef * 1.e6 * lal.PC_SI   # use *nonstandard* distance
        P.fmin = template_min_freq
        m1 = P.m1/lal.MSUN_SI
        m2 =P.m2/lal.MSUN_SI
        lambda1, lambda2 = P.lambda1,P.lambda2
        P.phiref=0.0
        P.psi=0.0
        P.incl = 0.0       # only works for aligned spins. Be careful.
        #P.fref = opts.reference_freq  # twice?
        if opts.approximant != "TaylorT4": # not default setting
            P.approx = lalsimutils.lalsim.GetApproximantFromString(opts.approximant)  # allow user to override the approx setting. Important for NR followup, where no approx set in sim_xml!
        if opts.approximant == "EccentricTD":
            P.phaseO = 3

    P = P_list[0]  # Load in the physical parameters of the injection.  
else:
 if opts.mass1 is not None and opts.mass2 is not None:
    m1, m2 = opts.mass1, opts.mass2
 else:
    raise RuntimeError('Missing --mass1 or --mass2')
 s1z=s2z=0.
 if opts.spin1z is not None:
    s1z = opts.spin1z
 if opts.spin2z is not None:
    s2z = opts.spin2z

# LISA check (should I change things here? Seems like it is for Neutron stars, so maybe not)
 lambda1, lambda2 = 0, 0
 if opts.eff_lambda is not None:
        lambda1, lambda2 = lalsimutils.tidal_lambda_from_tilde(m1, m2, opts.eff_lambda, opts.deff_lambda or 0)
 P = lalsimutils.ChooseWaveformParams(
        approx = lalsimutils.lalsim.GetApproximantFromString(opts.approximant),
    fmin = template_min_freq,
    radec = False,   # do NOT propagate the epoch later
    incl = 0.0,       # only works for aligned spins. Be careful.
    phiref = 0.0,
    theta = 0.0,
    phi = 0.0,
    psi = 0.0,
    m1 = m1 * lal.MSUN_SI,
    m2 = m2 * lal.MSUN_SI,
    lambda1 = lambda1,
    lambda2 = lambda2,
    s1z = s1z,
    s2z = s2z,
    ampO = opts.amp_order,
    fref = opts.reference_freq,
    tref = fiducial_epoch,
    dist = factored_likelihood.distMpcRef * 1.e6 * lal.PC_SI
    )
 P_list= [P]

# Time information not needed for FD data
if not(opts.h5_frame_FD):
  # User requested bounds for data segment
  if not (opts.data_start_time == None) and  not (opts.data_end_time == None):
      start_time =  opts.data_start_time
      end_time =  opts.data_end_time
      print( "Fetching data segment with start=", start_time)
      print( "                             end=", end_time)

  # Automatically choose data segment bounds so region of interest isn't corrupted
  # FIXME: Use estimate, instead of painful full waveform generation call here.
  else:
      approxTmp = P.approx
      LmaxEff = 2
      if opts.fmin_template_correct_for_lmax:
        LmaxEff=opts.l_max
      T_tmplt = lalsimutils.estimateWaveformDuration(P,LmaxEff) + 4  # Much more robust than the previous (slow, prone-to-crash approach)
  #     if (m1+m2)<30:
  #         P.approx = lalsimutils.lalsim.TaylorT4  # should not impact length much.  IMPORTANT because EOB calls will fail at the default sampling rate
  #         htmplt = lalsimutils.hoft(P)   # Horribly wasteful waveform gneeration solely to estimate duration.  Will also crash if spins are used.
  #         P.approx = approxTmp
  #         T_tmplt = - float(htmplt.epoch)    # ASSUMES time returned by hlm is a RELATIVE time that does not add tref. P.radec=False !
  #         print fiducial_epoch, event_time, htmplt.epoch
  #     else:
  #         print " Using estimate for waveform length in a high mass regime: beware!"
  #         T_tmplt = lalsimutils.estimateWaveformDuration(P) + 4
      T_seg = T_tmplt + T_spec + T_safety # Amount before and after event time
  #    if opts.use_external_EOB:
  #        T_seg *=2   # extra safety factor
      start_time = float(event_time) - T_seg
      end_time = float(event_time) + T_seg
      print( "Fetching data segment with start=", start_time)
      print( "                             end=", end_time)
      print( "\t\tEvent time is: ", float(event_time))
      print( "\t\tT_seg is: ", T_seg)

#
# Load in data and PSDs
#
data_dict, psd_dict = {}, {}

# this if statement will handle LISA FD data.
if opts.h5_frame_FD:
    for inst, chan in map(lambda c: c.split("="), opts.channel_name):
      print("Reading channel %s from cache %s" % (inst+":"+chan, opts.cache_file))
      data_dict[inst] = lisa_lalsimutils_compat.frame_h5_to_hoff(opts.cache_file,inst)
      print("Frequency binning: %e (1/deltaF = %f), length %d" % (data_dict[inst].deltaF, 1/data_dict[inst].deltaF,  data_dict[inst].data.length) )

# this can handle both LIGO and LISA TD data
else:  
  for inst, chan in map(lambda c: c.split("="), opts.channel_name):
      print( "Reading channel %s from cache %s" % (inst+":"+chan, opts.cache_file))
      data_dict[inst] = lisa_lalsimutils_compat.frame_data_to_non_herm_hoff(opts.cache_file,
              inst+":"+chan, start=start_time, stop=end_time,
              window_shape=opts.window_shape, deltaT=deltaT, h5_frame=opts.h5_frame)
      print( "Frequency binning: %f, length %d" % (data_dict[inst].deltaF,
              data_dict[inst].data.length) )

flow_ifo_dict = {}
if opts.fmin_ifo:
 for inst, freq_str in map(lambda c: c.split("="), opts.fmin_ifo):
    freq_low_here = float(freq_str)
    print( "Reading low frequency cutoff for instrument %s from %s" % (inst, freq_str), freq_low_here)
    flow_ifo_dict[inst] = freq_low_here

for inst, psdf in map(lambda c: c.split("="), opts.psd_file):
    print( "Reading PSD for instrument %s from %s" % (inst, psdf))
    psd_dict[inst] = lalsimutils.get_psd_series_from_xmldoc(psdf, inst)

    deltaF = data_dict[inst].deltaF
    psd_dict[inst] = lalsimutils.resample_psd_series(psd_dict[inst], deltaF)
    print( "PSD deltaF after interpolation %e" % psd_dict[inst].deltaF)

    # Implement PSD window rescaling: see T1900249
    if opts.psd_window_shape > 0 or opts.window_shape > 0:
      # Note this is just a constant scale factor in the likelihood, but it *can* have some effect on parameters if one windowing size is small (as it stretches/compresses the likleihood)
      window_fac_psd = lalsimutils.psd_windowing_factor(opts.psd_window_shape, len(psd_dict[inst].data.data))  # assume the windowing factor IS accounted for, and we have to undo it
      window_fac_data = lalsimutils.psd_windowing_factor(opts.window_shape, len(data_dict[inst].data.data))
      psd_dict[inst].data.data *= window_fac_data/window_fac_psd   # scale the windowing factor to be appropriate for our actual data window: the PSD was measured using a windowing different than the one we use

    # implement cutoff.  
    if inst in flow_ifo_dict.keys():
        if isinstance(psd_dict[inst], lal.REAL8FrequencySeries):
            psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(psd_dict[inst].data.length)
            psd_dict[inst].data.data[ psd_fvals < flow_ifo_dict[inst]] = 0 # 
        else:
            print('FAIL on PSD import')
            sys.exit(1)
#        elif isinstance(psd_dict[inst], pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries):  # for backward compatibility
#            psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(len(psd_dict[inst].data))
#            psd_dict[inst].data[psd_fvals < ifo_dict[inst]] =0


    assert psd_dict[inst].deltaF == deltaF

    # Highest freq. at which PSD is defined
    # if isinstance(psd_dict[inst],
    #         pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries):
    #     fmax = psd_dict[inst].f0 + deltaF * (len(psd_dict[inst].data) - 1)
    if isinstance(psd_dict[inst], lal.REAL8FrequencySeries):
        fmax = psd_dict[inst].f0 + deltaF * (psd_dict[inst].data.length - 1)

    # Assert upper limit of IP integral does not go past where PSD defined
    assert opts.fmax is None or opts.fmax<= fmax
    # Allow us to target a smaller upper limit than provided by the PSD. Important for numerical PSDs that turn over at high frequency
    if opts.fmax and opts.fmax < fmax:
        fmax = opts.fmax # fmax is now the upper freq. of IP integral

# Ensure data and PSDs keyed to same detectors
# For LISA we may only provde one PSD, so this check is not necessary.
if sorted(psd_dict.keys()) != sorted(data_dict.keys()) and not(opts.LISA):
    print >>sys.stderr, "Got a different set of instruments based on data and PSDs provided."

# Ensure waveform has same sample rate, padded length as data
#
# N.B. This assumes all detector data has same sample rate, length
#
# data_dict holds 2-sided FrequencySeries, so their length is the same as
# that of the original TimeSeries that was FFT'd = Nsamples
# Also, deltaF = 1/T, with T = the duration (in sec) of the original TimeSeries
# Therefore 1/(data.length*deltaF) = T/Nsamples = deltaT
P.deltaT = 1./ (data_dict[list(data_dict.keys())[0]].data.length * deltaF)
P.deltaF = deltaF
for Psig in P_list:
  Psig.deltaT = P.deltaT
  Psig.deltaf = P.deltaF


#
# Set up parameters and bounds
#

# PROBLEM: if too large, you can MISS a source. Does NOT need to be fixed for all masses *IF* the problem really has strong support
dmin = opts.d_min    # min distance
dmax = opts.d_max  # max distance FOR ANY SOURCE EVER. EUCLIDEAN



dmax_sampling_guess = dmax
distBoundGuess = dmax

print( "Recommended distance for sampling ", dmax_sampling_guess, " and probably near ", distBoundGuess, " smaller than  ", dmax)
print( "    (recommendation not yet used) ")
param_limits = { "psi": (0, 2*numpy.pi),
    "phi_orb": (0, 2*numpy.pi),
    "distance": (dmin, dmax),   # CAN LEAD TO CATASTROPHIC FAILURE if dmax is too large (adaptive fails - too few bins)
    "right_ascension": (0, 2*numpy.pi),
    "declination": (-numpy.pi/2, numpy.pi/2),
    "t_ref": (-t_ref_wind, t_ref_wind),
    "inclination": (0, numpy.pi)
}
if opts.internal_rotate_phase:
  param_limits['psi'] = (0, 4*numpy.pi)
  param_limits['phi_orb'] = (0, 4*numpy.pi)

#
# Parameter integral sampling strategy
#
print("\n###########################################################################################\nInitiating integrator\n###########################################################################################\n")
# Oracle for portfolio if needed


# Portfolio
use_portfolio=False
params = {}
sampler = mcsampler.MCSampler()
xpy_asarray_already = functools.partial(xpy_default.asarray,dtype=np.float64)

use_gmm_member=False   # set when a portfolio carries a GMM member (see the portfolio setup loop)
if opts.sampler_method == "adaptive_cartesian_gpu":
    print(" ILE: {}".format(opts.sampler_method))
    sampler = mcsamplerGPU.MCSampler()
    sampler.xpy = xpy_default
    sampler.identity_convert=identity_convert
    mcsampler  = mcsamplerGPU  # force use of routines in that file, for properly configured GPU-accelerated code as needed

    xpy_asarray_already = lambda x: x  # do nothing because we are already on the board for GPU-generated 

    if opts.sampler_xpy == "numpy":
      mcsampler.set_xpy_to_numpy()
      sampler.xpy= numpy
      sampler.identity_convert= lambda x: x
    
elif opts.sampler_method == "GMM":
    print(" ILE: {}".format(opts.sampler_method))
    sampler = mcsamplerEnsemble.MCSampler()
  
elif opts.sampler_method == 'AV':
    print(" ILE: {}".format(opts.sampler_method))
    opts.internal_use_lnL=True

    sampler = mcsamplerAdaptiveVolume.MCSampler(n_chunk=opts.n_chunk) # note larger is better, but keep GPU mem limits in mind
    sampler.xpy = xpy_default
    sampler.identity_convert=identity_convert
#    mcsampler  = mcsampler  # force use of routines in that file, for properly configured GPU-accelerated code as needed

    xpy_asarray_already = lambda x: x  # do nothing because we are already on the board for GPU-generated 

    if opts.sampler_xpy == "numpy":
      mcsampler.set_xpy_to_numpy()
      sampler.xpy= numpy
      sampler.identity_convert= lambda x: x

elif opts.sampler_method == "portfolio":
    use_portfolio=True
    opts.internal_use_lnL=True  # required, we only implement those scenarios right now
    sampler_list = []
    sampler_types = opts.sampler_portfolio

    # prep xpy, etc
    my_xpy = xpy_default
    my_identity_convert=identity_convert
    my_identity_convert_togpu=identity_convert_togpu
    print(" PORTFOLIO ", opts.sampler_portfolio)
    xpy_asarray_already = lambda x: x  # do nothing because we are already on the board for GPU-generated 
    if opts.sampler_xpy == "numpy":
      mcsampler.set_xpy_to_numpy()
      my_xpy= numpy
      my_identity_convert= lambda x: x
      my_identity_convert_togpu= lambda x: x
    for name in sampler_types:
        if name =='AV':
            sampler = mcsamplerAdaptiveVolume.MCSampler(n_chunk=opts.n_chunk) # enforce now, so provided for setup phase
        elif name =='GMM':
            sampler = mcsamplerEnsemble.MCSampler()
            # A GMM member needs the GMM-specific argument blocks below to run so its config is
            # forwarded.  This used to CLOBBER opts.sampler_method='GMM', which silently broke
            # every downstream `sampler_method == "portfolio"` test -- most importantly the L0
            # auto-rescue gate, which then NEVER FIRED for a portfolio carrying a GMM member --
            # and made a portfolio take GMM-only branches (e.g. return_lnI).  Flag it
            # non-destructively instead: sampler_method stays 'portfolio', and the GMM blocks
            # below key off `use_gmm_args` = standalone GMM OR a portfolio with a GMM member.
            # Ported from bin/integrate_likelihood_extrinsic_batchmode, which fixed this.
            use_gmm_member = True
        elif name == "adaptive_cartesian_gpu" or name == 'AC':
            sampler = mcsamplerGPU.MCSampler()
            mcsampler  = mcsamplerGPU  # force use of routines in that file, for properly configured GPU-accelerated code as needed
        elif name in mcsamplerPortfolio.known_pipelines:  # everything else, including nflow
            sampler =  mcsamplerPortfolio.known_pipelines[name]()
        else:
            # No else clause here meant an unrecognized name left `sampler` bound to its
            # previous value -- the plain MCSampler built before this chain, or, on the second
            # and later iterations, the PREVIOUS member -- and appended it silently.  The
            # portfolio then ran with a member the user never asked for, and a typo in
            # --sampler-portfolio produced a duplicate rather than an error.  Ported from
            # bin/integrate_likelihood_extrinsic_batchmode.  (The chain above is now elif for
            # the same reason: with plain `if`, a name matching no branch fell through every
            # test and reused whatever `sampler` still held.)
            raise Exception(" --sampler-portfolio: unknown member '{}'.  Known: AV, GMM, AC/adaptive_cartesian_gpu, {}".format(name, sorted(mcsamplerPortfolio.known_pipelines)))
        print('PORTFOLIO: adding {} '.format(name))
        # enable xpy for low level sampler as needed 
        if hasattr(sampler, 'xpy'):
          sampler.xpy = my_xpy
          sampler.identity_convert= my_identity_convert
          sampler.identity_convert_togpu=identity_convert_togpu
        sampler_list.append(sampler)
    sampler = mcsamplerPortfolio.MCSampler(portfolio=sampler_list)
    sampler.xpy = my_xpy
    sampler.identity_convert= my_identity_convert
    sampler.identity_convert_togpu= my_identity_convert_togpu
    # sampler weights will be CPU-typed, so don't change them

else:
    print(" ILE: original sampler")
    print(" ILE: {}".format(opts.sampler_method))

#
# Psi -- polarization angle
# sampler: uniform in [0, pi)
#
psi_sampler = mcsampler.ret_uniform_samp_vector_alt( 
    param_limits["psi"][0], param_limits["psi"][1])
psi_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, 
    param_limits["psi"][0], param_limits["psi"][1])
sampler.add_parameter("psi", 
    pdf = psi_sampler, 
    cdf_inv = psi_sampler_cdf_inv, 
    left_limit = param_limits["psi"][0], 
    right_limit = param_limits["psi"][1],
    prior_pdf = mcsampler.uniform_samp_psi,
    adaptive_sampling=opts.internal_rotate_phase or opts.force_adapt_all)

#
# Phi - orbital phase
# sampler: uniform in [0, 2*pi)
#
if not (opts.distance_marginalization and lookup_table["phase_marginalization"]):
    phi_sampler = mcsampler.ret_uniform_samp_vector_alt( 
        param_limits["phi_orb"][0], param_limits["phi_orb"][1])
    phi_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, 
        param_limits["phi_orb"][0], param_limits["phi_orb"][1])
    sampler.add_parameter("phi_orb",
        pdf = phi_sampler,
        cdf_inv = phi_sampler_cdf_inv,
        left_limit = param_limits["phi_orb"][0], 
        right_limit = param_limits["phi_orb"][1],
        prior_pdf = mcsampler.uniform_samp_phase,
        adaptive_sampling=opts.internal_rotate_phase or opts.force_adapt_all)


#
# inclination - angle of system angular momentum with line of sight
# sampler: cos(incl) uniform in [-1, 1)
#

adapt_extra_extrinsic=False
if (opts.sampler_method == "adaptive_cartesian_gpu"  or opts.sampler_method == 'AV' ) or use_portfolio or opts.force_adapt_all:  # this is a better/more stable/faster adaptive code, trust it to adapt in more extrinsic dimensions
  adapt_extra_extrinsic=True

if not opts.inclination_cosine_sampler:
 incl_sampler = mcsampler.cos_samp_vector # this is NOT dec_samp_vector, because the angular zero point is different!
 incl_sampler_cdf_inv = mcsampler.cos_samp_cdf_inv_vector
 sampler.add_parameter("inclination", 
    pdf = incl_sampler, 
    cdf_inv = incl_sampler_cdf_inv, 
    left_limit = param_limits["inclination"][0], 
    right_limit = param_limits["inclination"][1],
    prior_pdf = mcsampler.uniform_samp_theta)  # do not adapt in parameter going to zero at edge
else:
 incl_sampler = mcsampler.ret_uniform_samp_vector_alt(-1.0,1.0)
 incl_sampler_cdf_inv = lambda x: x*2.0-1.  #functools.partial(mcsampler.uniform_samp_cdf_inv_vector,-1,1) 
 sampler.add_parameter("inclination", 
    pdf = incl_sampler, 
    cdf_inv = incl_sampler_cdf_inv, 
    left_limit = -1, 
    right_limit = 1,
    prior_pdf = incl_sampler,
    adaptive_sampling=adapt_extra_extrinsic)

#
# Distance - luminosity distance to source in parsecs
# sampler: uniform distance over [dmin, dmax), adaptive sampling
#
print(f" Using {opts.d_prior} prior in distance.")
if not opts.distance_marginalization:
  dist_sampler = mcsampler.ret_uniform_samp_vector_alt( param_limits["distance"][0], param_limits["distance"][1])
  dist_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector,     param_limits["distance"][0], param_limits["distance"][1])
  #dist_sampler=functools.partial( mcsampler.uniform_samp_withfloor_vector, numpy.min([distBoundGuess,param_limits["distance"][1]]), param_limits["distance"][1], 0.001)
  dist_prior_pdf =   lambda x: x**2/(param_limits["distance"][1]**3/3.                   - param_limits["distance"][0]**3/3.) 
  if opts.d_prior == 'pseudo_cosmo':
    nm = priors_utils.dist_prior_pseudo_cosmo_eval_norm(param_limits["distance"][0],param_limits["distance"][1])
    dist_prior_pdf =functools.partial( priors_utils.dist_prior_pseudo_cosmo, nm=nm,xpy=xpy_default)
  elif opts.d_prior == 'uniform':
      dist_prior_pdf = dist_sampler
  elif opts.d_prior != 'Euclidean':
    print(" ==== WARNING UNKNOWN DISTANCE PRIOR === ")
    raise Exception('distance prior')
  #dist_sampler_cdf_inv=None
  sampler.add_parameter("distance", 
                        pdf = dist_sampler, 
                        cdf_inv = dist_sampler_cdf_inv,
                        left_limit = param_limits["distance"][0], 
                        right_limit = param_limits["distance"][1],
                        prior_pdf = dist_prior_pdf,
                        adaptive_sampling = (adapt_extra_extrinsic and not (opts.no_adapt or opts.no_adapt_distance)) or opts.force_adapt_all)

# 
# Rotate sky coordinates
#
if opts.internal_sky_network_coordinates:
  ifo_list = list(psd_dict)
  if not(opts.internal_sky_network_coordinates_raw):
    # remove V and K : the sky ring is almost always only HL
    if 'K1' in ifo_list:
      ifo_list.remove('K1')
    if 'V1' in ifo_list:
      ifo_list.remove('V1')
  # problem: ordering of psd_dict is rarely rational; almost always we want an HL network, rarely V
  if len(ifo_list) <2 or opts.LISA:
    opts.internal_sky_network_coordinates = False # stop using this code
  else:
    sky_rotations.assign_sky_frame(ifo_list[0], ifo_list[1], fiducial_epoch)
    frm = identity_convert_togpu(sky_rotations.frm)
    my_rotation = functools.partial(lalsimutils.polar_angles_in_frame_alt,frm,xpy=xpy_default) 
    my_rotation_cpu = functools.partial(lalsimutils.polar_angles_in_frame_alt,sky_rotations.frm,xpy=np) 

#
# Intrinsic parameters
#
sampler_lookup = {}
sampler_inv_lookup = {}
sampler_lookup['q'] = mcsampler.q_samp_vector   # only one intrinsic parameter possible
sampler_lookup['M'] = mcsampler.M_samp_vector
sampler_inv_lookup['q'] = mcsampler.q_cdf_inv_vector
sampler_inv_lookup['M'] = None # mcsampler.M_cdf_inv_vector

if opts.rom_use_basis and opts.rom_integrate_intrinsic:
 for p in intrinsic_param_names:
    indx = intrinsic_param_names.index(p)
    qmin,qmax = param_ranges[indx]
    q_pdf =mcsampler.ret_uniform_samp_vector_alt( qmin, qmax)   # sample uniformly by default
#    q_pdf =functools.partial(sampler_lookup[p], qmin, qmax)   # sample uniformly by default
    q_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, qmin,qmax) # sample uniformly by default
#    q_cdf_inv= None
    q_pdf_prior = functools.partial(sampler_lookup[param], qmin, qmax)  # true prior, from lookup
    sampler.add_parameter(p, 
        pdf=q_pdf, 
        cdf_inv=q_cdf_inv, 
        left_limit=qmin, right_limit=qmax,prior_pdf=q_pdf_prior, 
        adaptive_sampling = opts.adapt_intrinsic)



if False: #opts.skymap_file is not None:
    from ligo.skymap.io import fits as bfits
    #
    # Right ascension and declination -- use a provided skymap
    #
    smap, _ = bfits.read_sky_map(opts.skymap_file)
    # FIXME: Uncomment for 'mixed' map
    #smap = 0.9*smap + 0.1*numpy.ones(len(smap))/len(smap)
    ss_sampler = mcsampler.HealPixSampler(smap)
    #isotropic_2d_sampler = numpy.vectorize(lambda dec, ra: mcsampler.dec_samp_vector(dec)/2/numpy.pi)
    isotropic_bstar_sampler = numpy.vectorize(lambda dec, ra: 1.0/len(smap))

    # FIXME: Should the left and right limits be modified?
    sampler.add_parameter(("declination", "right_ascension"), 
        pdf = ss_sampler.pseudo_pdf,
        cdf_inv = ss_sampler.pseudo_cdf_inverse, 
        left_limit = (param_limits["declination"][0], param_limits["right_ascension"][0]),
        right_limit = (param_limits["declination"][1], param_limits["right_ascension"][1]),
        prior_pdf = isotropic_bstar_sampler)

else:
    #
    # Right ascension - angle in radians from prime meridian plus hour angle
    # sampler: uniform in [0, 2pi), adaptive sampling
    #
    ra_sampler = mcsampler.ret_uniform_samp_vector_alt(
        param_limits["right_ascension"][0], param_limits["right_ascension"][1])
    ra_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector,
        param_limits["right_ascension"][0], param_limits["right_ascension"][1])
    sampler.add_parameter("right_ascension", 
        pdf = ra_sampler, 
        cdf_inv = ra_sampler_cdf_inv, 
        left_limit = param_limits["right_ascension"][0],
        right_limit =  param_limits["right_ascension"][1],
        prior_pdf = mcsampler.uniform_samp_phase,
        adaptive_sampling = opts.force_adapt_all or ((not opts.no_adapt) and (not opts.internal_sky_network_coordinates)))  # TOO DANGEROUS to double-adapt in sky, ends up overconverging too easily. Just leave the whole sky ring in place if we want to do this, it's usually there and we don't lose that much.  If we have full localization, we should just remove the sky network coordinates argument!

    #
    # declination - angle in radians from the north pole piercing the celestial
    # sky sampler: cos(dec) uniform in [-1, 1), adaptive sampling
    #
    if not opts.declination_cosine_sampler:
     dec_sampler = mcsampler.dec_samp_vector
     dec_sampler_cdf_inv = mcsampler.dec_samp_cdf_inv_vector
     sampler.add_parameter("declination", 
        pdf = dec_sampler, 
        cdf_inv = dec_sampler_cdf_inv, 
        left_limit = param_limits["declination"][0], 
        right_limit = param_limits["declination"][1],
        prior_pdf = mcsampler.uniform_samp_dec,
        adaptive_sampling = opts.force_adapt_all or (not opts.no_adapt))
    else:
     # Sample uniformly in cos(polar_theta), =1 for north pole, -1 for south pole.
     # Propagate carefully in conversions: time of flight libraries use RA,DEC
     dec_sampler = mcsampler.ret_uniform_samp_vector_alt(-1.0,1.0)
     dec_sampler_cdf_inv = lambda x: x*2.0-1. # functools.partial(mcsampler.uniform_samp_cdf_inv_vector,-1,1)
     sampler.add_parameter("declination", 
        pdf = dec_sampler, 
        cdf_inv = dec_sampler_cdf_inv, 
        left_limit = -1, 
        right_limit = 1,
        prior_pdf = dec_sampler,
        adaptive_sampling = opts.force_adapt_all or (not opts.no_adapt))

# LISA: What will change here if we don't marginalize in time for LISA?
if not opts.time_marginalization:
    #
    # tref - GPS time of geocentric end time
    # sampler: uniform in +/-2 ms window around estimated end time 
    #
    tref_sampler = mcsampler.ret_uniform_samp_vector_alt(
        param_limits["t_ref"][0], param_limits["t_ref"][1])

    tref_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, 
                                            param_limits["t_ref"][0], param_limits["t_ref"][1])
    sampler.add_parameter("t_ref", 
                          pdf = tref_sampler, 
                          cdf_inv = tref_sampler_cdf_inv, 
                          left_limit = param_limits["t_ref"][0], 
                          right_limit = param_limits["t_ref"][1],
                          prior_pdf = functools.partial(mcsampler.uniform_samp_vector, param_limits["t_ref"][0], param_limits["t_ref"][1]))


# skymap oracle
oracleRS=None
if opts.skymap_file:
  # Read file
  from RIFT.misc.reference_samples import ReferenceSamples
  rs_object = ReferenceSamples()
  if opts.skymap_file.endswith('fits'):
    rs_object.from_skymap_fits(fname=opts.skymap_file,cos_dec=opts.declination_cosine_sampler)
  else:
      rs_object.from_ascii(fname=opts.skymap_file, reference_params=['ra', 'dec'])
      if opts.declination_cosine_sampler:
        rs_object.reference_params[:,1] = np.arccos(rs.reference_params[:,1])
      # relabel
      rs_object.reference_params  = ['right_ascension', 'declination']

  # Pass to resampling oracle
  from RIFT.integrators.unreliable_oracle.resampling import ResamplingOracle
  oracleRS = ResamplingOracle()
  for p in sampler.params_ordered:
        oracleRS.add_parameter(p,pdf=None, left_limit =sampler.llim[p], right_limit = sampler.rlim[p])
  oracleRS.setup(reference_samples=rs_object.reference_samples, reference_params=rs_object.reference_params)

#
# Determine pinned and non-pinned parameters
#

# LISA
# I am calling ecliptic latitude (beta) as declination and ecliptic longitude (lambda) as right-ascension.

pinned_params = get_pinned_params(opts)
unpinned_params = get_unpinned_params(opts, sampler.params)

if opts.LISA and opts.lisa_fixed_sky:
    lisa_sky_lamda = float(opts.ecliptic_longitude)
    lisa_sky_beta = float(opts.ecliptic_latitude)
    pinned_params["right_ascension"] = 0.0
    pinned_params["declination"] = 0.0
if opts.LISA and not(opts.lisa_fixed_sky): # read the grid for sky location co-ordinate, but fixing it to prevent sampler from wasting time sampling skylocation.
    pinned_params["right_ascension"] = 0.0
    pinned_params["declination"] = 0.0


print( "{0:<25s} {1:>5s} {2:>5s} {3:>20s} {4:<10s}".format("parameter", "lower limit", "upper limit", "pinned?", "pin value"))
plen = len(sorted(sampler.params, key=lambda p: len(p))[-1])
for p in sampler.params:
    if p in pinned_params:
        pinned, value = True, "%1.3g" % pinned_params[p]
    else:
        pinned, value = False, ""

    if isinstance(p, tuple):
        for subp, subl, subr in zip(p, sampler.llim[p], sampler.rlim[p]):
            subp = subp + " "*min(0, plen-len(subp))
            print( "|{0:<25s} {1:>1.3g}   {2:>1.3g} {3:>20s} {4:<10s}".format(subp, subl, subr, str(False), ""))
    else:
        p = p + " "*min(0, plen-len(p))
        print("{0:<25s} {1:>5g} {2:>5g} {3:>20s} {4:<10s}".format(p, sampler.llim[p], sampler.rlim[p], str(pinned), value))

# Special case: t_ref is assumed to be relative to the epoch
if "t_ref" in pinned_params:
    pinned_params["t_ref"] -= float(fiducial_epoch)

#
# Provide convergence tests
# FIXME: Currently using hardcoded thresholds, poorly hand-tuned
#
test_converged = {}

#
# Merge options into one big ol' kwargs dict
#

pinned_params.update({ 
    # Iteration settings and termination conditions
    "n": min(opts.n_chunk, n_max), # Number of samples in a chunk
    "nmax": n_max, # Total number of samples to draw before termination
    "neff": n_eff, # Total number of effective samples to collect before termination

    "convergence_tests" : test_converged,    # Dictionary of convergence tests

    # Adaptive sampling settings
    "tempering_exp": opts.adapt_weight_exponent if not opts.no_adapt else 0.0, # Weights will be raised to this power to prevent overconvergence
    "tempering_log": opts.adapt_log,
    "tempering_adapt": opts.adapt_adapt,

    "floor_level": opts.adapt_floor_level if not opts.no_adapt else 0.0, # The new sampling distribution at the end of each chunk will be floor_level-weighted average of a uniform distribution and the (L^tempering_exp p/p_s)-weighted histogram of sampled points.
    "history_mult": 10, # Multiplier on 'n' - number of samples to estimate marginalized 1-D histograms
    "n_adapt": 100 if not opts.no_adapt else 0, # Number of chunks to allow adaption over

    # Verbosity settings
    "verbose": True, #not opts.rom_integrate_intrinsic, 
    "extremely_verbose": False, 

    # Sample caching
    "save_intg": opts.save_samples, # Cache the samples (and integrand values)?
    "igrand_threshold_deltalnL": opts.save_deltalnL, # Threshold on distance from max L to save sample
    "igrand_threshold_p": opts.save_P, # Threshold on cumulative probability contribution to cache sample
    "igrand_fairdraw_samples": opts.fairdraw_extrinsic_output,
    "igrand_fairdraw_samples_max": opts.n_eff
})
if opts.sampler_method == "adaptive_cartesian_gpu":
  pinned_params.update({"save_no_samples":True})   # do not exhaust GPU memory with MC samples!  
return_lnL=False
if opts.sampler_method=="GMM"  and opts.internal_use_lnL:
  return_lnL=True
  pinned_params.update({"use_lnL":True,"return_lnI":True})
if opts.sampler_method =="adaptive_cartesian_gpu" and opts.internal_use_lnL:
  return_lnL=True
  pinned_params.update({"use_lnL":True})
if opts.sampler_method =="AV" and opts.internal_use_lnL:
  # AV integrates in log space natively (integrate() is a thin wrapper over integrate_log);
  # without this, --internal-use-lnL --sampler-method AV passed the ok_lnL_methods check but
  # silently did nothing, so exp(lnL) overflowed at high SNR when no logarithm offset was set.
  #
  # PORTED FROM THE MAIN DRIVER, where this branch already exists.  It was missing here, and
  # the drift audit could not see it: a missing `if` branch is not a FUNC/OPTION/CONST/ATTR,
  # so it produces no gap item.  High-SNR is the LISA MBHB regime, which is exactly the case
  # the main driver's comment describes.
  return_lnL=True
  pinned_params.update({"use_lnL":True})
if opts.sampler_method =="portfolio":
  return_lnL=True
  pinned_params.update({"use_lnL":True})

# What the sampler will actually STORE in _rvs['integrand'], derived from the pinned params
# above rather than from the CLI.  This is not the same predicate as opts.internal_use_lnL:
# that option is accepted for adaptive_cartesian_gpu and portfolio too (see the branches
# directly above), which set use_lnL WITHOUT return_lnI and therefore still store linear L.
# Keying the weight helpers off the option would compute L + ln p - ln p_s for those, which
# is the failure the main driver documents at ln_weights_from_rvs.
rvs_integrand_is_lnL = bool(pinned_params.get("return_lnI", False))


# ---------------------------------------------------------------------------------------
# Fair-draw weighting helpers.  Ported from bin/integrate_likelihood_extrinsic_batchmode
# (PR #87); see test/expensive_before_merging/integrators/RVS_FAIRDRAW_AUDIT.md.
#
# WHY THESE ARE HERE, given this driver has no .dgrid/.dslice/proposal-breadcrumb exports
# (the three consumers whose double-weighting PR #87 actually fixed): this driver DOES set
# igrand_fairdraw_samples from --fairdraw-extrinsic-output, so its _rvs can be a fair draw,
# and every shared sampler already sets the provenance marker at its rebind.  The marker was
# arriving here and nothing was reading it.  The helpers are the correct thing for the next
# person to reach for, which is the whole argument of the audit's Recommendation 1.
#
# KEEP IN STEP WITH THE MAIN DRIVER.  These are deliberate copies, not an import, because the
# two drivers are a deliberate fork; audit_lisa_driver_drift.py is what makes the copy visible.
# ---------------------------------------------------------------------------------------
def _rvs_lnL_convention(use_lnL=None):
    """Resolve the stored-'integrand' convention for a helper call.

    Returns the explicit argument when given, else the run's `rvs_integrand_is_lnL`.  Falls
    back to False (the historical linear reading) when that global is absent, which is what
    happens when these helpers are lifted out of the driver by the unit tests.  Read through
    globals() rather than by name so a missing global cannot become a NameError swallowed by
    a caller's bare `except Exception`.
    """
    if use_lnL is not None:
        return bool(use_lnL)
    return bool(globals().get('rvs_integrand_is_lnL', False))


def ln_weights_from_rvs(rvs, convert=None, use_lnL=False):
    """THE importance log-weight of an _rvs record: lnL + ln(prior) - ln(sampling_prior).

    ONE definition, because the alternative has already cost us.  A stored 'log_weights'
    column does not mean the same thing in every sampler: mcsamplerPortfolio stores the true
    importance weight, but mcsamplerGPU stores tempering_exp*lnL + ln p - ln p_s -- the
    ADAPTATION weight, with --adapt-weight-exponent baked in.  That exponent is not 1 in
    production and --no-adapt drives it to 0, removing the likelihood from the column
    entirely.  A consumer preferring that cache silently reweights its output by L^(e-1).

    So the cache is never read here: the weight is DERIVED from the canonical components --
    log form first, then the linear (mcsamplerEnsemble) form, out-of-support rows -inf.
    Raises when neither set is present: an explicit failure beats a plausible wrong number.

    `use_lnL` is REQUIRED to read the linear form correctly, because mcsamplerEnsemble reuses
    'integrand' for BOTH conventions (it stores lnL when given return_lnI).  Taking log() of
    lnL compresses tens of nats into log(tens), leaving an almost flat weight vector, and the
    positivity cut is wrong in that mode too: non-positive means a low-likelihood point, not
    a rejected one, so `ig > 0` would discard every sample with lnL <= 0.

    PASS THE STORED CONVENTION, NOT THE CLI OPTION -- `rvs_integrand_is_lnL`, not
    `opts.internal_use_lnL`.  In this driver the two genuinely differ: --internal-use-lnL is
    also accepted for adaptive_cartesian_gpu and portfolio, which set use_lnL without
    return_lnI and still store linear L.
    """
    conv = convert if convert is not None else (lambda x: x)
    if all(k in rvs for k in ('log_integrand', 'log_joint_prior', 'log_joint_s_prior')):
        return (numpy.asarray(conv(rvs['log_integrand']), dtype=float)
                + numpy.asarray(conv(rvs['log_joint_prior']), dtype=float)
                - numpy.asarray(conv(rvs['log_joint_s_prior']), dtype=float))
    if all(k in rvs for k in ('integrand', 'joint_prior', 'joint_s_prior')):
        ig = numpy.asarray(conv(rvs['integrand']), dtype=float)
        jp = numpy.asarray(conv(rvs['joint_prior']), dtype=float)
        js = numpy.asarray(conv(rvs['joint_s_prior']), dtype=float)
        out = numpy.full(len(ig), -numpy.inf)
        if use_lnL:
            # 'integrand' already holds lnL: do not log it again, do not cut on its sign.
            keep = numpy.isfinite(ig) & (jp > 0) & (js > 0)
            out[keep] = ig[keep] + numpy.log(jp[keep]) - numpy.log(js[keep])
        else:
            keep = (ig > 0) & (jp > 0) & (js > 0)
            out[keep] = numpy.log(ig[keep]) + numpy.log(jp[keep]) - numpy.log(js[keep])
        return out
    raise Exception("cannot build importance weights from sampler._rvs (keys={})".format(
        sorted(rvs.keys())))


def _rvs_len(rvs):
    """Rows in a raw `_rvs` column dict -> int.

    ONE row-count rule, and it lives with the record (`rvs_record.n_rows`).  Flattening
    whichever column came first was wrong for the ORDINARY case, not a corner: `_rvs` is
    seeded parameters-first, and a combined parameter is stored (ndim, N) under a TUPLE key,
    so any run registering one reported ndim*N.  Here that number is the length the pooled
    export's weight vector is checked against, so the check failed and the pooled record went
    out weight-mixed -- the exact degradation `_export_rvs_equal_weight` exists to prevent.

    Imported INSIDE the function deliberately: the test harnesses exec these helpers out of
    the driver into a bare namespace, so a module-level name here would have to be threaded
    through every one of them -- and this staying a one-line delegation is the point.
    """
    from RIFT.integrators.rvs_record import n_rows as _n_rows_of_columns
    return _n_rows_of_columns(rvs)


def _rvs_is_export_resample(sampler):
    """True when the ROWS of _rvs were drawn in proportion to weight.

    Set by the samplers at the rebind itself, so it means "the draw FIRED", which is NOT the
    same predicate as `opts.fairdraw_extrinsic_output`: the draw is skipped when it would not
    shrink the record (n_extr >= len(_rvs)), and then the rows are still the retained set
    carrying real importance weights.  Keying off the CLI flag would flatten those -- the same
    class of error in the other direction.

    SURVIVES POOLING by design, and this driver now does pool (--mc-error-replicas): a pooled
    record built from fair-drawn replicas still has posterior-resampled rows, so anything that
    must not re-weight them keeps seeing True here.  Whether the record is GLOBALLY
    equal-weight is a different question; see _rvs_is_equal_weight.
    """
    return bool(getattr(sampler, '_rvs_is_fairdraw', False))


def _rvs_is_equal_weight(sampler):
    """True when EVERY row of _rvs carries the same posterior weight.

    Two properties, deliberately not one flag:

      rows resampled  -- each row drawn proportional to w   (per-BLOCK property)
      equal weight    -- the record as a whole is uniform   (property of the WHOLE record)

    A single fair draw has both.  A POOLED record has the first and not the second, because
    pooling weights block k by the replica evidence Z_k/K.  Conflating them broke two things
    in opposite directions in the main driver (audit Finding 6), which is why the split is
    carried over here even though this driver has no pooling yet.
    """
    return (bool(getattr(sampler, '_rvs_is_fairdraw', False))
            and not bool(getattr(sampler, '_rvs_is_pooled', False)))


def ln_weights_for_posterior(rvs, sampler, convert=None, use_lnL=None):
    """The weights to use when treating an _rvs record as a POSTERIOR SAMPLE SET.

    NOT the same question as `ln_weights_from_rvs`, which answers "what is the importance
    weight of this record" and is always right about that.  The question here is "how should
    these rows be weighted to represent the posterior", and the answer depends on whether the
    fair draw already did it.

    A fair-drawn record was resampled WITH REPLACEMENT proportional to w, so its rows are
    already an equal-weight draw from the posterior.  Weighting them by w again applies w^2
    and over-concentrates the result -- measured at a 13% shift in the posterior mean of a
    weight-correlated coordinate (verify_skew.py).

    So: uniform (zero log-weight) for a fair-drawn record, the derived importance weight
    otherwise.  Returns a float array the length of the record.

    `use_lnL` is passed THROUGH UNRESOLVED, exactly as in the main driver: a caller that
    omits it gets the linear reading, not the run's convention.  That is a trap in both
    drivers, and it is deliberately reproduced rather than fixed here -- a helper of the
    same name behaving differently in the two forked drivers would be a worse defect than
    the one it fixes.  Callers must pass `use_lnL=rvs_integrand_is_lnL`, or route through
    `_rvs_lnL_convention` first, the way the main driver's call sites do.
    """
    # MIGRATION (DESIGN_rvs_naming.md), the first consumer to move.  This is the exact
    # site where the one-flag-two-questions defect lived, so it is the one worth converting
    # first: `is_equal_weight()` is a named question rather than two booleans a caller has to
    # combine, and it cannot be answered with the wrong one.
    #
    # The flags stay as the fallback while the other six samplers are unconverted -- and while
    # both exist they MUST agree, which is asserted directly in test_rvs_record.py rather than
    # left as a comment, because "two sources of truth" is the risk this migration runs.
    _rec = _rvs_record_for(sampler, rvs)
    if _rec is not None:
        if _rec.is_equal_weight():
            return numpy.zeros(_rvs_len(rvs), dtype=float)
        # THE WEIGHT ITSELF now comes from the record, not from ln_weights_from_rvs -- which is
        # the point of the record: it knows its own convention, so there is no `use_lnL` to
        # thread through and no way for a caller to pass the wrong one.
        #
        # Verified equivalent before switching, not after: the two implementations were fuzzed
        # against each other over 1200 randomized records spanning all three column families
        # with NaN / -inf / 0 sprinkled through every column.  That found a REAL divergence
        # first -- log_weights() had been computing lnL + ln(pi) - ln(q) term by term, which
        # yields NaN where the canonical form's conjunctive keep-mask yields -inf -- and it is
        # fixed there rather than papered over here.
        return numpy.asarray(_rec.log_weights(convert=convert), dtype=float)
    if _rvs_is_equal_weight(sampler):
        return numpy.zeros(_rvs_len(rvs), dtype=float)
    return numpy.asarray(ln_weights_from_rvs(rvs, convert=convert, use_lnL=use_lnL),
                         dtype=float)
use_gmm_args = (opts.sampler_method == "GMM") or use_gmm_member
if use_gmm_args:   # standalone GMM, or a portfolio carrying a GMM member
    n_step =pinned_params["n"]
    n_max_blocks = ((1.0*int(opts.n_max))/n_step)
    # pairing coordinates for adaptive integration: see definition of order below
    #    (distance,inclination)
    #    (ra, dec)
    #    (psi) (orb_phase)   # only do 1d adaptivity there, 
#    gmm_dict = {tuple([2]):None,tuple([4]):None,(3,5):None,(0,1):None} 
#    comp_dict = {tuple([2]):1,tuple([4]):1,(3,5):3,(0,1):4} 
    default_phipsi = mcsamplerEnsemble.create_wide_single_component_prior( [ param_limits['psi'], param_limits['phi_orb']])
#    pair_sky = sampler_param_tuple(sampler, ['right_ascension','declination'])
    pair_ra_dec =   sampler_param_tuple(sampler, ["right_ascension", "declination"])
    n_d =  2  # 2 distance-inclination lobes by default
    n_sky = 4
    n_phase =4
    adapt_phase = False
    if 'distance' in sampler.params:
         pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination'])
    else:
         pair_d_incl = sampler_param_tuple(sampler, ['inclination'])
         n_d = 1 # one peak in distance by default
    if 'phi_orb' in sampler.params:
         pair_phi_psi = sampler_param_tuple(sampler, ['psi', "phi_orb"])
    else:
         pair_phi_psi = sampler_param_tuple(sampler, ['psi'])
         adapt_phase = True
         n_phase = 2
    gmm_dict = {pair_ra_dec:None,pair_d_incl:None,pair_phi_psi:default_phipsi} 
    gmm_adapt = {pair_ra_dec: True, pair_d_incl: True, pair_phi_psi: False}
    if opts.internal_rotate_phase:
      gmm_adapt[pair_phi_psi] = True
#    gmm_dict = {tuple([2]):None,tuple([4]):None,(3,5):None,tuple([0]):None,tuple([1]):None} 
#    if len(psd_dict.keys()) > 2:
#        n_sky=2  # presumably we have localized better, avoid failure modes
#    comp_dict = {tuple([2]):1,tuple([4]):1,(3,5):3,tuple([0]):n_sky,tuple([1]):n_sky} 
    comp_dict = {pair_ra_dec:n_sky,pair_d_incl:n_d,pair_phi_psi:n_phase} 
    extra_args = {'n_comp':comp_dict,'max_iter':n_max_blocks,'gmm_dict':gmm_dict, 'gmm_adapt':gmm_adapt}  # made up for now, should adjust
    print("GMM:",extra_args)
    print("GMM:",sampler.params_ordered)
    # if opts.distance_marginalization:
    #   if lookup_table["phase_marginalization"]:
    #     gmm_dict = {pair_ra_dec:None,(3,4):None,(0,):None}
    #     gmm_adapt = {pair_ra_dec: True, (3,4): True, (0,): True}
    #     comp_dict = {pair_ra_dec:n_sky,tuple([4]):1,(0,):2}
    #   else:
    #     gmm_dict = {pair_ra_dec:None,tuple([4]):None,(0,1):None}
    #     gmm_adapt = {pair_ra_dec:True,tuple([4]):True,(0,1):True}
    #     comp_dict = {pair_ra_dec:n_sky,tuple([4]):1,(0,1):4}
#      extra_args = {'n_comp':comp_dict,'max_iter':n_max_blocks,'gmm_dict':gmm_dict}  # made up for now, should adjust

    # force adapt in all parameters, if we request it. Requires pass-by-reference in above
    if opts.force_adapt_all:
        for param in gmm_adapt:
            gmm_adapt[param] = True
    pinned_params.update(extra_args)

    # cannot pin params at present
    # unpinned params MUST be full call signature of likelihood, IN ORDER
    if not(opts.time_marginalization):
      unpinned_params =['right_ascension','declination', 't_ref','phi_orb','inclination','psi','distance']
    else:
      unpinned_params =['right_ascension','declination', 'phi_orb','inclination','psi','distance']
    if opts.distance_marginalization:
      unpinned_params.remove('distance')
      if lookup_table["phase_marginalization"]:
        unpinned_params.remove('phi_orb')
if opts.sampler_method == "AV":
  return_lnL=True
  pinned_params.update( { 'enforce_bounds':True})  # don't go out of range : choose integer bin sizes


# set up sampler, as needed.  Mainly for portfolio integrator
if use_portfolio:
      print(" PORTFOLIO : setup")
      if opts.sampler_portfolio_args:
          print(" PRE_EVAL", opts.sampler_portfolio_args)
          #opts.sampler_portfolio_args = list(map(lambda x: eval(' "{}" '.format(x)), opts.sampler_portfolio_args))
          opts.sampler_portfolio_args = list(map(eval, opts.sampler_portfolio_args))
          # confirm all are dict
          for indx in range(len(opts.sampler_portfolio_args)):
            if not(isinstance(opts.sampler_portfolio_args[indx], dict)):
                print(indx,opts.sampler_portfolio_args[indx]) 
          print(" ARGS ", opts.sampler_portfolio_args)
      # Assemble freeze-policy overrides from the CLI.  Only include options the user actually
      # set (None = unset) so the sampler keeps its built-in defaults otherwise.  The two VARAHA
      # flags are mutually exclusive; --portfolio-varaha-can-freeze wins if both are given.
      _freeze_policy_kwargs = {}
      if opts.portfolio_grace_iters is not None:
          _freeze_policy_kwargs['portfolio_grace_iters'] = opts.portfolio_grace_iters
      if opts.portfolio_revive_period is not None:
          _freeze_policy_kwargs['portfolio_revive_period'] = opts.portfolio_revive_period
      if opts.portfolio_freeze_wt is not None:
          _freeze_policy_kwargs['portfolio_freeze_wt'] = opts.portfolio_freeze_wt
      if opts.portfolio_varaha_can_freeze:
          _freeze_policy_kwargs['portfolio_varaha_never_freeze'] = False
      elif opts.portfolio_varaha_never_freeze:
          _freeze_policy_kwargs['portfolio_varaha_never_freeze'] = True
      # adaptive-probe draw allocation (OPT-IN; off by default in the sampler)
      if opts.portfolio_adaptive_alloc:
          _freeze_policy_kwargs['portfolio_adaptive_alloc'] = True
      if opts.portfolio_varaha_min_frac is not None:
          _freeze_policy_kwargs['portfolio_varaha_min_frac'] = opts.portfolio_varaha_min_frac
      if opts.portfolio_varaha_max_frac is not None:
          _freeze_policy_kwargs['portfolio_varaha_max_frac'] = opts.portfolio_varaha_max_frac
      if opts.portfolio_weight_clip is not None:
          _freeze_policy_kwargs['portfolio_weight_clip'] = opts.portfolio_weight_clip
      if opts.portfolio_quality_signal is not None:
          _freeze_policy_kwargs['portfolio_quality_signal'] = opts.portfolio_quality_signal
      if opts.portfolio_alloc_exponent is not None:
          _freeze_policy_kwargs['portfolio_alloc_exponent'] = opts.portfolio_alloc_exponent
      if opts.portfolio_probe_period is not None:
          _freeze_policy_kwargs['portfolio_probe_period'] = opts.portfolio_probe_period
      print(" PORTFOLIO freeze-policy overrides: ", _freeze_policy_kwargs)
      sampler.setup(portfolio_args=opts.sampler_portfolio_args, **_freeze_policy_kwargs, **pinned_params) # directly pass all parameters set above to low-level portfolios.  In particular, GMM setup

# initialize sampler, before we call integrate, so we can seed it
if opts.sampler_method == 'adaptive_cartesian_gpu' and opts.skymap_file:
  sampler.setup()


# ---------------------------------------------------------------------------------------
# L0 auto-rescue.  Ported from bin/integrate_likelihood_extrinsic_batchmode (PR #79/#84/#87);
# see test/expensive_before_merging/integrators/RVS_FAIRDRAW_AUDIT.md Findings 1 and 5.
#
# ONE DELIBERATE STRUCTURAL DIVERGENCE FROM THE MAIN DRIVER.  There the rescue is inline in
# the single analyze_event.  This driver has TWO -- analyze_event_LISA (used with --LISA) and
# analyze_event (the non-LISA fallback) -- each with its own integrate call and export block,
# already ~50% duplicated.  Inlining the rescue twice would create a third copy to keep in
# step, which is the failure mode this whole exercise exists to prevent.  So the block lives
# in _maybe_l0_rescue below and both call it.  The helpers are byte-identical to main's and
# are pinned that way by test_lisa_l0_rescue.py.
# ---------------------------------------------------------------------------------------
def _rvs_record_for(sampler, rvs):
    """The record describing THESE columns, or None.  See DESIGN_rvs_naming.md.

    THE IDENTITY CHECK IS THE POINT.  A record holds a reference to a column dict that other
    code replaces in place, so "the sampler has a record" and "the record describes the rows I
    am holding" are different questions -- the same shape as everything else in this file's
    history. A record that has fallen out of step is not consulted; the caller falls back to
    the provenance flags, which are maintained separately and are still correct.

    One lookup rather than the check repeated per consumer, for the reason the reserve lookup
    was centralised in #87: two copies of a guard drift.
    """
    # `samples()` is the public accessor; the getattr guard is for an object that predates the
    # mixin (an old pickle, a test double), not for the six samplers, all of which have it.
    _get = getattr(sampler, 'samples', None)
    rec = _get() if callable(_get) else None
    if rec is None or getattr(rec, 'columns', None) is not rvs:
        return None
    return rec


def _sampler_keeps_records(sampler):
    """Does this sampler populate `_rvs_record` at all?  See DESIGN_rvs_naming.md.

    A PRODUCER's question, not a consumer's, and deliberately a different function from
    `_rvs_record_for`.  The pooling step is about to REPLACE `sampler._rvs`, so asking "does a
    record describe the rows I hold" is the wrong question there -- it would be answered `None`
    and the pooled record would silently not be built.  What it needs to know is whether this
    sampler participates in the record scheme at all.

    Two questions, two names.  That is the entire lesson of this file's last four review rounds.
    """
    # PARTICIPATION, not "is one present right now".  Every sampler clears _rvs_record at the
    # top of integrate(), so a replica that raised leaves None behind while the sampler is still
    # a full participant -- and keying on presence would silently skip building the pooled
    # record for it.  Ask whether the sampler implements the scheme at all.
    return isinstance(sampler, SamplerOutputMixin) or (
        callable(getattr(sampler, 'samples', None))
        and callable(getattr(sampler, 'set_samples', None)))


def _internal_record_of(sampler):
    """This pass's record, marked INTERNAL for threading -> RvsRecord or None.

    Replica pooling needs each block's record to derive that block's weights with the right
    convention.  Marking them internal is the difference between "we had to hand the structure
    back" and "this is now something consumers may use": set_samples() refuses an internal
    record, so nothing on this list can reappear from samples().
    """
    _get = getattr(sampler, 'samples', None)
    rec = _get() if callable(_get) else None
    return rec.as_internal() if rec is not None else None


def _rebound_record(sampler, columns):
    """A copy of the sampler's record whose `.columns` is `columns` -> RvsRecord or None.

    Snapshot/restore installs a COPY of the column dict, so a record still pointing at the
    original would fail every identity check and silently do nothing.
    """
    _get = getattr(sampler, 'samples', None)
    rec = _get() if callable(_get) else None
    if rec is None:
        return None
    out = rec.snapshot()
    out.columns = columns
    return out


def _lw_of(rvs, record, use_lnL):
    """Importance log-weights for `rvs`, preferring a record that describes it.

    ONE resolver, so the two estimators below cannot drift in which source they trust.  A
    record is used only when its `.columns` IS this dict: `_rvs` is copied and replaced all
    over this file, and a record describing different columns must not be believed.  Otherwise
    fall back to the canonical derivation with the stored convention -- the two are verified
    equivalent by a randomized comparison in test_rvs_record.py, so this is a source choice,
    not a semantics choice.
    """
    if record is not None and getattr(record, 'columns', None) is rvs:
        return numpy.asarray(record.log_weights(), dtype=float)
    return numpy.asarray(ln_weights_from_rvs(rvs, use_lnL=_rvs_lnL_convention(use_lnL)),
                         dtype=float)


def _lnZ_of_rvs(rvs, already_pooled=True, use_lnL=None, record=None):
    """log of the evidence implied by an _rvs record.

    For a POOLED record the weights already carry their 1/(K n_k) factor, so the estimate is the
    plain sum; for a single run it is the mean.  Returns None when the weights cannot be rebuilt.
    """
    try:
        try:
            lw = _lw_of(rvs, record, use_lnL)
        except Exception:
            return None
        lw = lw[numpy.isfinite(lw)]
        if lw.size == 0:
            return None
        m = numpy.max(lw)
        tot = m + numpy.log(numpy.sum(numpy.exp(lw - m)))
        return float(tot if already_pooled else tot - numpy.log(lw.size))
    except Exception:
        return None


def _kish_neff_of_rvs(rvs, use_lnL=None, record=None):
    """Kish effective sample size of an _rvs record, or None if the weights are not reconstructible."""
    try:
        try:
            lw = _lw_of(rvs, record, use_lnL)
        except Exception:
            return None
        lw = lw[numpy.isfinite(lw)]
        if lw.size == 0:
            return None
        lw = lw - numpy.max(lw)
        w = numpy.exp(lw)
        return float(numpy.sum(w) ** 2 / numpy.sum(w ** 2))
    except Exception:
        return None


def _lnZ_of_reserve_or_rvs(sampler, rvs, reserve=None):
    """lnZ of a completed pass, from the points it RETAINED where that is available.

    The rescue's reject gate compares the warm pass's lnZ against the cold pass's, and both
    were read out of _rvs -- which the fair draw has already replaced with
    min(n_extr, 1.5*eff_samp, 1.5*neff) rows resampled WITH REPLACEMENT, proportional to
    weight.  That is not a smaller unbiased sample of the same estimator, it is a DIFFERENT
    and biased one: _lnZ_of_rvs forms logsumexp(w)/n, so drawing n rows proportional to w
    returns something near max(w) rather than mean(w), high by roughly

        log(n_retained / eff_samp)

    and the two passes are drawn at wildly different n and eff_samp.  The gate was reading a
    multi-nat artifact of its own two subsample sizes as evidence that the warm seed had
    missed mass.

    Falls back to the old _rvs reading when no reserve was kept, so the comparison degrades to
    the previous behaviour rather than to no gate at all.

    Returns (lnZ, source) -- the caller MUST check that both sides came from the same source,
    because the two readings are not interchangeable.
    """
    _res = reserve if reserve is not None else getattr(sampler, '_warm_seed_reserve', None)
    if isinstance(_res, dict) and 'log_joint_prior' in _res and 'log_joint_s_prior' in _res:
        try:
            # NOT _lnZ_of_rvs: it averages over the rows it is handed, and the reserve is
            # neither the draw set nor a uniform sample of it -- non-finite rows were dropped
            # and the remainder may have been capped.  lnZ_from_reserve restores the original
            # proposal-draw normalization from n_finite/n_retained.  Without it a PORTFOLIO
            # reading is high by ~log(n_retained/n_finite), ~11 nats on a collapsed pass, and
            # the error does NOT cancel in the gate: the cold and warm passes have different
            # finite fractions, so it is the difference of two different-sized errors.
            _v = mcsamplerAdaptiveVolume.lnZ_from_reserve(_res)
            if _v is not None and numpy.isfinite(_v):
                return _v, 'retained'
        except Exception:
            pass
    return _lnZ_of_rvs(rvs, already_pooled=False,
                       record=_rvs_record_for(sampler, rvs)), 'fairdraw'


def _snapshot_pass_state(sampler, res, var, neff, dict_return, rvs=None):
    """Everything that must move TOGETHER when a completed pass is put back -> dict.

    THE POINT IS THE WORD "everything".  A pass is described by more than its samples, and the
    reject path used to restore only some of it: `_rvs`, the estimate and `dict_return` went
    back to the cold pass while `_warm_seed_reserve` was left holding the REJECTED warm cloud.
    In the main driver that stayed latent until --sampler-sequential-warmstart began seeding
    the next intrinsic point from the reserve, at which point a rejected, truncated warm pass
    became the seed for the next point -- the exact failure the reject gate exists to prevent,
    reintroduced one attribute over.  THAT OPTION DOES NOT EXIST IN THIS DRIVER YET, so the
    reserve restore is pre-emptive here; it is also what makes porting the capture safe, which
    is why it lands first.

    So the snapshot carries the reserve and the fair-draw marker as well, including the
    per-member reserves: `_warm_seed_reserve_for` falls through to `portfolio_realizations`,
    so restoring only the aggregate would leave that fallback pointing at the warm pass.
    """
    return dict(
        rvs=(dict(sampler._rvs) if rvs is None else rvs),
        res=res, var=var, neff=neff, dict_return=dict_return,
        warm_seed_reserve=getattr(sampler, '_warm_seed_reserve', None),
        rvs_is_fairdraw=bool(getattr(sampler, '_rvs_is_fairdraw', False)),
        rvs_is_pooled=bool(getattr(sampler, '_rvs_is_pooled', False)),
        # The record too.  A stale one is already declined by _rvs_record_for's identity check,
        # so this is belt-and-braces -- but "everything describing the pass moves together" is
        # the invariant, and carving an exception into it is how round 1 happened.
        # REBOUND to the snapshot's columns.  The record held a reference to the LIVE dict, and
        # the restore installs a COPY -- so storing it as-is produced a record whose identity
        # check could never match, i.e. inert rather than belt-and-braces.  Rebinding makes it
        # describe what is actually put back.
        rvs_record=_rebound_record(sampler, dict(sampler._rvs) if rvs is None else rvs),
        member_reserves=[getattr(_m, '_warm_seed_reserve', None)
                         for _m in list(getattr(sampler, 'portfolio_realizations', []) or [])],
    )


def _restore_pass_state(sampler, state):
    """Undo of _snapshot_pass_state -> (res, var, neff, dict_return).

    Both callers (the reject path and the exception handler) go through here, so the set of
    attributes that travels with a restored pass cannot drift between them.
    """
    sampler._rvs = state['rvs']
    sampler._warm_seed_reserve = state['warm_seed_reserve']
    sampler._rvs_is_fairdraw = state['rvs_is_fairdraw']
    sampler._rvs_is_pooled = state['rvs_is_pooled']
    if callable(getattr(sampler, 'set_samples', None)):
        sampler.set_samples(state.get('rvs_record'))
    _members = list(getattr(sampler, 'portfolio_realizations', []) or [])
    for _m, _r in zip(_members, state.get('member_reserves', [])):
        _m._warm_seed_reserve = _r
    return state['res'], state['var'], state['neff'], state['dict_return']


def _warm_seed_reserve_for(sampler):
    """The retained-sample reserve a completed pass left behind, or None.

    THE ONE LOOKUP FOR SEED CONSUMERS.  In the main driver there are two -- the L0 auto-rescue
    (which re-seeds a collapsed pass from its own peak) and the --sampler-sequential-warmstart
    capture (which seeds the NEXT intrinsic point).  ONLY THE RESCUE EXISTS IN THIS DRIVER; the
    shared lookup is kept so the pair cannot drift once the capture is ported.  Both otherwise
    fall back to sampler._rvs, which by then has been rebound to a fair-draw subset
    taken WITH REPLACEMENT -- and the whole point of the reserve is that on the collapsed pass
    a warm start exists for, that subset is a handful of rows several of which are the same
    point twice.

    A PORTFOLIO keeps the reserve on the aggregate, not on its members, but a bare AV member
    can be the one that has it; check the sampler first, then its realizations.

    COLUMN ORDER MUST MATCH or the seed is scrambled: the reserve stores X in the column order
    of the sampler that built it, and a seed handed to bootstrap_from_samples is read
    positionally against params_ordered.  A mismatch is silent and produces a seed in the
    wrong coordinates, so decline the reserve rather than use it.

    NOT SHARED WITH `_lnZ_of_reserve_or_rvs` above, deliberately.  That one reads only lnL and
    the two prior columns, never X, so a column-order mismatch is harmless to it and declining
    would throw away a good lnZ reading and silently downgrade the gate to its fallback.
    """
    _res = getattr(sampler, '_warm_seed_reserve', None)
    if _res is None:
        for _m in list(getattr(sampler, 'portfolio_realizations', []) or []):
            _res = getattr(_m, '_warm_seed_reserve', None)
            if _res is not None:
                break
    if _res is not None and list(_res.get('params_ordered', [])) != list(sampler.params_ordered):
        return None
    return _res


def _warm_seed_geometry(sampler):
    """Which columns a warm seed must span, and the box it must lie in -> (axes, lo, hi).

    The seed is judged on the ADAPTIVE axes, because those are the only ones the live-volume
    grid resolves (the rest get a single bin), and that is the set the [AV COLLAPSE] report
    counts against.  Ask the sampler that will consume the seed rather than assuming all
    dimensions: with --force-adapt-all they coincide, without it a rank test over every column
    would demand a seed span directions the grid cannot resolve and puff for nothing.

    A PORTFOLIO has no adaptive axes of its own -- they live on its AV-style members -- so fall
    through to the first member that can answer.  If nobody can, every column it is.
    """
    _lo = np.array([sampler.llim[p] for p in sampler.params_ordered], dtype=float)
    _hi = np.array([sampler.rlim[p] for p in sampler.params_ordered], dtype=float)
    for _s in [sampler] + list(getattr(sampler, 'portfolio_realizations', []) or []):
        if hasattr(_s, 'warm_seed_axes'):
            try:
                return list(_s.warm_seed_axes()), _lo, _hi
            except Exception:
                pass
    return list(range(len(sampler.params_ordered))), _lo, _hi


def _clear_warm_state(sampler):
    """Clear a warm-start seed AND any grid it installed, reaching PORTFOLIO MEMBERS too.

    `sampler._warm = None` alone is not enough for mcsamplerPortfolio: `_warm` and the
    contracted AV grid live on each MEMBER, and portfolio.integrate_log() does not rerun each
    member's setup(), so the next point would silently draw from the PREVIOUS point's
    contracted live volume.  If the new point's support falls outside it, lnZ is biased low
    with a healthy-looking n_eff and no error.  Portfolio exposes clear_warm_state();
    everything else keeps the old behaviour.
    """
    # Deliberately NOT wrapped in try/except.  A reset that quietly did not happen leaves the
    # next point drawing from the previous point's contracted grid -- the exact silent bias
    # this guards against -- so a failure must abort the point rather than degrade to a log
    # line nobody reads.
    if hasattr(sampler, 'clear_warm_state'):
        sampler.clear_warm_state()
    else:
        sampler._warm = None
        sampler._warm_applied = False


def _maybe_load_av_state(sampler):
    """Warm-start this integration from a saved AV live-volume state (--sampler-load-state)."""
    try:
        if opts.sampler_load_state and hasattr(sampler, 'load_state'):
            print("  warm-start: loading saved sampler state from", opts.sampler_load_state)
            sampler.load_state(opts.sampler_load_state)
    except Exception as _e_ls:
        print("  AV state load skipped (", _e_ls, ")")


def _maybe_save_av_state(sampler):
    """Persist the adapted live-volume state for reuse by later instances/iterations."""
    if opts.sampler_method == 'AV' and opts.sampler_save_state and hasattr(sampler, 'save_state'):
        if not getattr(sampler, '_av_state_reuse_safe', True):
            print("  AV: not saving live-volume state from a rejected/failed rescue")
            return
        try:
            sampler.save_state(opts.sampler_save_state)
            print("  AV: saved live-volume state to", opts.sampler_save_state)
        except Exception as _e_ss:
            print("  AV: could not save state (", _e_ss, ")")


def _maybe_enable_anisotropic_bins(sampler):
    """Opt-in per-axis bin allocation, on the AV sampler AND any AV portfolio members."""
    if getattr(opts, 'sampler_anisotropic_bins', False):
        _aniso_targets = [sampler] + list(getattr(sampler, 'portfolio_realizations', []))
        for _t in _aniso_targets:
            if hasattr(_t, 'anisotropic_bins'):
                _t.anisotropic_bins = True
        print("  AV: anisotropic per-axis bin allocation ENABLED")


def _extract_mc_diag(dd):
    dd = dd if isinstance(dd, dict) else {}
    return dd.get('pareto_khat', None), dd.get('sigma_lnZ_block', None), dd.get('n_ESS', None), dd.get('lnZ_ci90', None)


def _pool_replica_rvs(rep_rvs, sampler, rep_lnZ=None, already_resampled=False, use_lnL=None,
                      records=None):
    """Concatenate the replicas' samples into one correctly-weighted set.

    Each replica k is an independent importance-sampling estimate with weights w_ki and its own
    sample count n_k, and the reported evidence is the linear mean (1/K) sum_k Z_k.  The posterior
    that matches THAT estimator is the concatenation with weights w_ki/(K n_k) -- equivalently the
    importance weight against the pooled proposal q'_ki = q_ki * K * n_k, which is the actual
    density of "pick a replica uniformly, then one of its n_k draws".  So the K*n_k factor goes
    into the sampling prior, where every downstream weight computation already accounts for it.

    Falls back to the first replica if the record shape is unexpected: a degraded export is
    recoverable, a silently mis-weighted one is not.

    `use_lnL` is the stored convention of the RAW ('integrand') columns -- see
    `_rvs_lnL_convention`.  It matters here because this function REWRITES joint_s_prior to force a
    block's weights, and the equation to solve is convention-dependent (see below).
    """
    _lnL_here = _rvs_lnL_convention(use_lnL)
    # `already_resampled` may be a single bool or a PER-REPLICA sequence.  It has to be the
    # latter in general: each pass decides independently whether to fair-draw (the draw is
    # skipped when it would not shrink that pass's record), so one global flag either flattens
    # a replica whose weights are genuine, or leaves a resampled replica double-weighted.  A
    # mixture of raw and resampled replicas is the normal case near the n_extr boundary.
    _ar_list = (list(already_resampled)
                if isinstance(already_resampled, (list, tuple, numpy.ndarray))
                else None)
    # PER-REPLICA RECORDS, threaded in so each block's lnZ is derived with ITS OWN convention
    # instead of one `use_lnL` asserted over the whole set.  These are INTERNAL: they are
    # plumbing for this function, marked as such, and refused by set_samples() so they cannot
    # escape through the public samples() accessor.  Having had to pass the structure around is
    # not a reason for anyone else to reach for it.
    _rec_list = list(records) if records is not None else None
    # Drop empty records in LOCKSTEP with their metadata.  The filter used to run on rep_rvs
    # alone, so a single empty replica shifted every later block against its own lnZ -- and
    # would now shift it against its own resampled flag too.
    _keep = [i for i, r in enumerate(rep_rvs) if r]
    rep_rvs = [rep_rvs[i] for i in _keep]
    if rep_lnZ is not None:
        rep_lnZ = [rep_lnZ[i] for i in _keep if i < len(rep_lnZ)]
    if _ar_list is not None:
        _ar_list = [_ar_list[i] for i in _keep if i < len(_ar_list)]
    if _rec_list is not None:
        _rec_list = [_rec_list[i] for i in _keep if i < len(_rec_list)]

    def _block_record(i, r):
        """The record for block i, but only if it describes THAT block's columns."""
        if _rec_list is None or i >= len(_rec_list):
            return None
        rec = _rec_list[i]
        return rec if getattr(rec, 'columns', None) is r else None

    def _block_resampled(i):
        if _ar_list is not None:
            return bool(_ar_list[i]) if i < len(_ar_list) else False
        return bool(already_resampled)

    def _block_column(k, v):
        """One block's column for key `k`, in the layout the KEY implies.
        THE KEY SAYS WHERE THE ROW AXIS IS -- the same rule `_rvs_len` delegates to
        (rvs_record._column_n_rows), applied to the rows themselves.  A combined parameter is
        stored (ndim, N) under a TUPLE key, so ravelling it and concatenating on axis 0 turns
        it into ONE 1-D column of length ndim*sum(N) while the scalar columns have sum(N) rows.
        Consumers still require (ndim, N) -- the sample exporter unpacks the combined sky column
        as `samples["latitude"], samples["longitude"] = samples[("declination",
        "right_ascension")]` -- so --mc-error-replicas produced a malformed record and could
        abort the export.  Per-row columns keep the flatten they always had.
        """
        v = numpy.asarray(v)
        return numpy.atleast_2d(v) if isinstance(k, tuple) else numpy.atleast_1d(v).ravel()

    def _empty_column(k):
        """No block contributed any rows -- an empty column that still has the key's LAYOUT.
        Handing back a bare `array([])` for a tuple key would fail to unpack in the exporter
        for the shape reason above rather than for the real one (there are no samples).
        """
        if not isinstance(k, tuple):
            return numpy.array([])
        ndim = _block_column(k, sampler.identity_convert(rep_rvs[0][k])).shape[0]
        return numpy.empty((ndim, 0))

    if len(rep_rvs) <= 1:
        return rep_rvs[0] if rep_rvs else {}
    # `already_resampled` -- the records are FAIRDRAW output.  Those samples were already drawn in
    # proportion to their own importance weights, so reusing those weights applies them a second
    # time and the pooled block follows w^2 instead of w.  Renormalizing to Z_k/K fixes the block's
    # SCALE but not its SHAPE, so it does not help here.  A fairdraw block is an equal-weight draw
    # from its own posterior, so that is what it must contribute: constant weights within the
    # block, summing to Z_k/K.
    #
    # DO NOT assume the records are raw importance samples.  integrate() may have thresholded or
    # fairdraw-resampled _rvs before we see it, in which case sum_i w_ki over the RETAINED rows is
    # no longer Z_k * n_k and a 1/n_k rescale would mis-weight the replica (a fairdraw record is
    # already posterior-resampled, so scaling it by its retained length weights it twice).  When
    # the reported per-replica lnZ is available, renormalize each block so it contributes exactly
    # Z_k/K -- correct whether the rows are raw, pruned or resampled, since only their RELATIVE
    # weights need be right.
    keys = set(rep_rvs[0])
    for r in rep_rvs[1:]:
        keys &= set(r)
    log_key = 'log_joint_s_prior' if 'log_joint_s_prior' in keys else None
    lin_key = 'joint_s_prior' if (log_key is None and 'joint_s_prior' in keys) else None
    if log_key is None and lin_key is None:
        print(" [mc error] pooling skipped: no sampling-prior column in the replica records; "
              "exporting the FIRST replica (consistent weights, fewer samples)")
        return rep_rvs[0]
    K = len(rep_rvs)
    out = {}
    try:
        cols = {k: [] for k in keys}
        for _i, r in enumerate(rep_rvs):
            n_k = _rvs_len(r)
            if n_k <= 0:
                continue
            _flat_block = False
            if _block_resampled(_i) and rep_lnZ is not None and _i < len(rep_lnZ) \
                    and numpy.isfinite(rep_lnZ[_i]):
                # equal weights within the block, summing to Z_k/K
                _flat_block = True
                _target_lw = float(rep_lnZ[_i]) - numpy.log(float(K)) - numpy.log(float(n_k))
                scale = 0.0
            elif rep_lnZ is not None and _i < len(rep_lnZ) and numpy.isfinite(rep_lnZ[_i]):
                # target: this block's weights sum to Z_k/K
                _cur = _lnZ_of_rvs(r, already_pooled=True, use_lnL=_lnL_here,
                                   record=_block_record(_i, r))
                if _cur is None or not numpy.isfinite(_cur):
                    scale = numpy.log(float(K) * float(n_k))
                else:
                    scale = _cur - (float(rep_lnZ[_i]) - numpy.log(float(K)))
            else:
                scale = numpy.log(float(K) * float(n_k))
            if _flat_block and log_key is not None:
                # force lw_i = log_integrand + log_joint_prior - log_joint_s_prior == _target_lw
                _li = numpy.atleast_1d(numpy.asarray(
                    sampler.identity_convert(r['log_integrand']), dtype=float)).ravel()
                _lp = numpy.atleast_1d(numpy.asarray(
                    sampler.identity_convert(r['log_joint_prior']), dtype=float)).ravel()
                _forced = _li + _lp - _target_lw
            for k in keys:
                v = _block_column(k, sampler.identity_convert(r[k]))
                if _flat_block and log_key is not None and k == log_key:
                    v = _forced
                elif _flat_block and lin_key is not None and k == lin_key:
                    # Same forcing for a RAW-field record: choose joint_s_prior so the
                    # reconstructed weight is exactly _target_lw.  WHICH equation that is depends
                    # on the convention 'integrand' is stored in -- the same ambiguity
                    # ln_weights_from_rvs handles:
                    #   linear:  lw = log(ig) + log(jp) - log(js)  ->  js = ig*jp/exp(target)
                    #   log:     lw = ig      + log(jp) - log(js)  ->  js = exp(ig + log(jp) - target)
                    # Applying the linear form to an lnL record gives js < 0 for every row with
                    # lnL < 0 -- a NEGATIVE proposal density -- and the block weights it produces
                    # are not constant at all, which is the entire point of the flat block.  Worse,
                    # it corrupts the canonical columns BEFORE ln_weights_from_rvs ever reads them,
                    # so fixing the helper alone does not rescue this path.
                    _ig = numpy.atleast_1d(numpy.asarray(
                        sampler.identity_convert(r['integrand']), dtype=float)).ravel()
                    _jp = numpy.atleast_1d(numpy.asarray(
                        sampler.identity_convert(r['joint_prior']), dtype=float)).ravel()
                    if _lnL_here:
                        v = numpy.exp(_ig + numpy.log(_jp) - _target_lw)
                    else:
                        v = _ig * _jp / numpy.exp(_target_lw)
                elif k == log_key:
                    v = v + scale
                elif k == lin_key:
                    # The linear counterpart of 'log_joint_s_prior += scale'.  This used to be a
                    # hardcoded K*n_k, which is only the FALLBACK value of `scale` -- so whenever a
                    # reported per-replica lnZ was available the raw-field path silently skipped
                    # the renormalization the log path applied, and a pruned or thresholded replica
                    # was mis-weighted.  exp(scale) reduces to K*n_k in the fallback case.
                    v = v * numpy.exp(scale)
                cols[k].append(v)
        for k in keys:
            # ...and concatenate along THAT row axis: axis 1 for a combined (ndim, N) parameter,
            # axis 0 for everything else.  One rule, stated in _block_column, applied twice.
            out[k] = (numpy.concatenate(cols[k], axis=1 if isinstance(k, tuple) else 0)
                      if cols[k] else _empty_column(k))
        # CACHED WEIGHTS MUST FOLLOW THE COMPONENTS.  _rvs may carry a precomputed 'log_weights'
        # (mcsamplerPortfolio writes one), and the .dgrid and calibration-posterior exporters
        # PREFER it -- they only fall back to log_integrand + log_joint_prior - log_joint_s_prior
        # when it is absent.  Concatenating the per-replica caches unchanged would hand those
        # scientific outputs the ORIGINAL weights while the estimate used the corrected ones:
        # replica rebalancing ignored, and fairdraw blocks double-weighted again in exactly the
        # products this pooling exists to make consistent.  Recompute from the canonical columns.
        # Rebuild through the ONE canonical definition rather than a second inline copy of it:
        # the copy that used to live here carried the same linear-only assumption as the helper's
        # old second branch, so under the log convention it re-logged lnL and cut on its sign --
        # writing exactly the flattened weights the exporters prefer.
        try:
            _lw_pooled = ln_weights_from_rvs(out, use_lnL=_lnL_here)
        except Exception:
            _lw_pooled = None
        if _lw_pooled is not None:
            if 'log_weights' in out:
                out['log_weights'] = _lw_pooled
            if 'weights' in out:
                out['weights'] = numpy.exp(_lw_pooled - numpy.max(_lw_pooled[numpy.isfinite(_lw_pooled)]))
        elif 'log_weights' in out or 'weights' in out:
            # cannot rebuild them -> DROP, so consumers fall through to whatever components exist
            # rather than silently trusting a stale cache.
            out.pop('log_weights', None)
            out.pop('weights', None)
            print(" [mc error] pooled record: dropped stale cached weights (components unavailable"
                  " to rebuild them); consumers will reconstruct from what remains")
    except Exception as e:
        print(" [mc error] pooling failed ({}); exporting the FIRST replica".format(e))
        return rep_rvs[0]
    return out


def _export_rvs_equal_weight(rvs, sampler, use_lnL=None):
    """The version of a POOLED record that may be written to the SimInspiral XML.

    THE XML KEEPS NO WEIGHT THIS DRIVER WRITES.  xmlutils maps 'joint_prior'/'joint_s_prior'
    onto alpha2/alpha3 -- which the ILE export below overwrites with zeros -- and the
    log_joint_* columns the pool actually carries have no mapping at all.  So every exported
    row is read downstream with the same weight, whatever the record says.

    For an ordinary run that is the long-standing convention (the export is a fair draw, or is
    treated as one).  For a POOLED record it is wrong in a NEW way: _pool_replica_rvs gives
    block k weights summing to Z_k/K, deliberately unequal BETWEEN blocks, so equal-weight rows
    mix the replicas by ROW COUNT instead of by evidence -- silently discarding exactly the
    disagreement the replicas were run to measure, in the one output a human looks at.

    So convert here rather than hope: draw rows in proportion to the pool's reconstructed
    posterior weights, turning weights the format drops into row multiplicities it keeps.

    Returns its argument UNCHANGED for any record that is not a pooled mixture (every
    non-replica run is untouched), and on any failure to rebuild or apply the weights -- with a
    message, because a weighted export is a real defect, not a silent degradation.
    """
    if not bool(getattr(sampler, '_rvs_is_pooled', False)):
        return rvs
    _conv = getattr(sampler, 'identity_convert', None)
    n = _rvs_len(rvs)
    if n <= 1:
        return rvs

    def _bail(why):
        print(" [mc error] pooled export left AS-IS ({}); the XML preserves no weight column,"
              " so its consumers will mix the replicas by row count".format(why))
        return rvs

    try:
        # ln_weights_for_posterior, not ln_weights_from_rvs: it is the "how should these rows be
        # weighted as a posterior" question.  The pooled marker is what makes it answer with the
        # reconstructed per-row weights instead of zeros -- a flat block contributes constant
        # weights summing to Z_k/K (already an equal-weight draw, correctly scaled), a raw block
        # its genuine importance weights, likewise scaled.  Resolve the stored convention here:
        # the helper deliberately passes use_lnL through unresolved, as the main driver's does.
        lw = numpy.asarray(ln_weights_for_posterior(rvs, sampler, convert=_conv,
                                                    use_lnL=_rvs_lnL_convention(use_lnL)),
                           dtype=float).ravel()
    except Exception as e:
        return _bail(e)
    if lw.size != n:
        return _bail("weights are {} long for {} rows".format(lw.size, n))
    _ok = numpy.isfinite(lw)
    if not numpy.any(_ok):
        return _bail("no finite weights")
    w = numpy.zeros(n, dtype=float)
    w[_ok] = numpy.exp(lw[_ok] - numpy.max(lw[_ok]))
    _tot = float(numpy.sum(w))
    if not numpy.isfinite(_tot) or _tot <= 0:
        return _bail("weights do not sum to anything usable")
    w = w / _tot
    # HOW MANY ROWS.  The Kish n_eff of the pooled weights, capped at the rows available: the
    # honest count, and the same quantity the samplers' own fair draw caps on.  Drawing the full
    # K*n_k rows instead would report K times the independent information whenever one replica
    # dominates -- the case pooling exists to expose.
    n_out = int(min(n, max(1, int(round(1.0 / float(numpy.sum(w ** 2)))))))
    # SYSTEMATIC resampling, not multinomial: one uniform offset, then n_out equally spaced
    # positions through the cumulative weight.  Unbiased in the same way, but each block gets its
    # evidence share of the rows deterministically rather than with O(sqrt(n)) draw noise on top
    # of the replica scatter being measured -- and when the replicas agree (equal weights,
    # n_out == n) it returns every row exactly once, where a bootstrap would duplicate ~37% of
    # them for nothing.
    cdf = numpy.cumsum(w)
    cdf[-1] = 1.0
    pos = (numpy.random.uniform() + numpy.arange(n_out)) / float(n_out)
    idx = numpy.clip(numpy.searchsorted(cdf, pos, side='left'), 0, n - 1)
    out = {}
    for k, v in rvs.items():
        try:
            arr = numpy.asarray(_conv(v) if _conv is not None else v)
        except Exception as e:
            return _bail("column {!r}: {}".format(k, e))
        # Index the LAST axis: _rvs may hold tuple-keyed pairs stored as (2, n).
        if arr.ndim < 1 or arr.shape[-1] != n:
            return _bail("column {!r} is not row-shaped".format(k))
        out[k] = arr[..., idx]
    print(" [mc error] pooled export re-drawn to equal weight: {} rows -> {} (weights the XML"
          " cannot carry are now row multiplicities)".format(n, n_out))
    return out


def _reject_if_collapsed(dd, stage):
    """Apply --reject-collapsed-live-volume to whatever the CURRENT verdict is.

    Called TWICE, in both drivers: once on the first run and again on the replica POOL,
    because replication can turn a healthy first run into a collapsed pool and gating only
    the first would bypass the flag for exactly the case pooling introduces.  Both calls live
    in _maybe_replicate_for_mc_error, which is why analyze_event must not gate directly.
    """
    if not opts.reject_collapsed_live_volume:
        return
    if not (isinstance(dd, dict) and dd.get('live_volume_collapsed', False)):
        return
    # Route through the ordinary failure path, so the caller skips this binary and writes no
    # result row -- the pre-fix outcome, but for a stated reason.
    _exc = mcsamplerAdaptiveVolume.LiveVolumeCollapse if mcsampler_AV_ok else RuntimeError
    raise _exc(
        "extrinsic integration collapsed ({}): live volume degenerated ({}); "
        "--reject-collapsed-live-volume is set, so this event is being dropped "
        "rather than exported".format(stage, dd.get('collapse_reason', '')))


def _report_and_gate_collapse(dict_return, stage="first run"):
    """Announce a collapsed live volume, then apply the rejection gate."""
    _collapsed = bool(dict_return.get('live_volume_collapsed', False)) if isinstance(dict_return, dict) else False
    _collapse_reason = (dict_return or {}).get('collapse_reason', '') if isinstance(dict_return, dict) else ''
    if _collapsed:
        print(" [mc error] *** LIVE VOLUME COLLAPSED *** {}".format(_collapse_reason))
        print(" [mc error] this event's lnZ and exported samples are NOT a fair draw from the posterior.")
    _reject_if_collapsed(dict_return, stage)


def _maybe_replicate_for_mc_error(sampler, res, var, neff, dict_return,
                                  log_res, sqrt_var_over_res,
                                  like_to_integrate, unpinned_params, pinned_params,
                                  lnL_offset=0.0):
    """MC-error replication + replica pooling.

    Returns (res, var, neff, log_res, sqrt_var_over_res, dict_return); returns them unchanged
    when nothing triggers, so the call site is one unconditional assignment.

    A module-level helper rather than inline (as in the main driver) because THIS DRIVER HAS
    TWO analyze_event variants -- inlining ~200 lines twice would be a third copy to keep in
    step.  Same reason as _maybe_l0_rescue.

    IT OWNS BOTH COLLAPSE-GATE CALLS.  The main driver gates once on the first run and again
    on the POOLED verdict, because replication can turn a healthy first run into a collapsed
    pool; gating only the first silently bypasses --reject-collapsed-live-volume for exactly
    the case pooling introduces.  Callers must therefore NOT call _report_and_gate_collapse
    themselves -- this does it.

    `lnL_offset` is the event's lnL_offset, used only for printing
    absolute lnZ.
    """
    _khat, _sig_block, _n_ess, _ci90 = _extract_mc_diag(dict_return)
    if _sig_block is not None and numpy.isfinite(_sig_block) and _sig_block > sqrt_var_over_res:
        print(" [mc error] sigma_lnZ raised to the between-chunk scatter: {:.4f} -> {:.4f}".format(float(sqrt_var_over_res), float(_sig_block)))
        sqrt_var_over_res = float(_sig_block)
    if _khat is not None:
        print(" [mc error] Pareto k-hat = {:.3f}{}".format(float(_khat), "  (> {:.2f}: weight tail unresolved; the reported sigma is a LOWER BOUND)".format(opts.mc_error_khat_trigger) if _khat > opts.mc_error_khat_trigger else ""))
    if _ci90 is not None:
        print(" [mc error] bootstrap lnZ 5/50/95 quantiles: {}".format(numpy.array2string(numpy.asarray(_ci90) + lnL_offset, precision=4)))

    # First-run collapse report AND gate (see the docstring: the pooled gate is below).
    _report_and_gate_collapse(dict_return, "first run")

    # AV live-volume state is persisted HERE, and the position is doubly constrained:
    #   * AFTER the first-run gate, so a collapsed grid that --reject-collapsed-live-volume
    #     rejects is never written.  Otherwise the event is correctly dropped while the NEXT
    #     intrinsic point warm-starts from the degenerate volume via --sampler-load-state,
    #     biased toward the surviving mode with nothing flagged.
    #   * BEFORE the replica loop, because afterwards the sampler holds the LAST replica's
    #     adapted grid rather than the run being reported.
    # The main driver satisfies only the second: it saves ~80 lines above its own first-run
    # gate, so a collapsed grid CAN be persisted there.  Deliberate divergence, and the main
    # driver should take the same reordering -- flagged rather than changed here.
    _maybe_save_av_state(sampler)
    _collapsed = bool(dict_return.get('live_volume_collapsed', False)) if isinstance(dict_return, dict) else False
    _collapse_reason = (dict_return or {}).get('collapse_reason', '') if isinstance(dict_return, dict) else ''

    _trigger_reasons = []
    if _collapsed and opts.mc_error_replicas > 0:
        _trigger_reasons.append('live volume collapsed ({})'.format(_collapse_reason))
    if opts.mc_error_replicas > 0:
        _neff_target = pinned_params.get('neff', None)
        if sqrt_var_over_res > opts.mc_error_sigma_trigger:
            _trigger_reasons.append('sigma={:.3f}>{:.2f}'.format(float(sqrt_var_over_res), opts.mc_error_sigma_trigger))
        if _khat is not None and _khat > opts.mc_error_khat_trigger:
            _trigger_reasons.append('khat={:.2f}>{:.2f}'.format(float(_khat), opts.mc_error_khat_trigger))
        if _n_ess is not None and _n_ess < opts.mc_error_ess_trigger:
            _trigger_reasons.append('ESS={:.1f}<{:g}'.format(float(_n_ess), opts.mc_error_ess_trigger))
        if _neff_target is not None and float(neff) < float(_neff_target):
            _trigger_reasons.append('neff={:.1f}<target {:g} (ran out of nmax)'.format(float(neff), float(_neff_target)))
    if _trigger_reasons:
        print(" [mc error] REPLICATING ({} extra cold runs): ".format(int(opts.mc_error_replicas)) + ", ".join(_trigger_reasons))
        _rep_lnZ = [float(log_res)]; _rep_sig = [float(sqrt_var_over_res)]; _rep_neff = [float(neff)]
        # Keep EVERY replica's samples.  Exporting only the highest-n_eff replica made the
        # posterior disagree with the evidence it was reported alongside: lnZ is the linear mean
        # over K replicas, so the samples must represent that same mixture.  Worse, n_eff is the
        # wrong selector -- it measures weight CONCENTRATION, not coverage, so a mode-collapsed
        # replica scores HIGHEST and would be the one exported.  (Measured elsewhere in this work:
        # the copy with the highest n_eff in its arm was the most biased, 11 nats low.)
        _rep_rvs = [sampler._rvs]
        # Provenance PER REPLICA, captured beside the record it describes.  The CLI flag is
        # not this: each pass decides independently whether to fair-draw (skipped when it
        # would not shrink that pass's record), so passing opts.fairdraw_extrinsic_output to
        # the pooler either flattens a replica whose importance weights are genuine, or leaves
        # a resampled replica double-weighted.  Near the n_extr boundary a run can produce a
        # MIXTURE of raw and resampled replicas, which one global boolean cannot describe.
        _rep_fairdraw = [bool(getattr(sampler, '_rvs_is_fairdraw', False))]
        # Collapse status must be aggregated over EVERY replica that ends up in the pool.
        # The exported posterior is the pooled mixture, so one collapsed replica taints it
        # even if the first run was healthy -- and the status sidecar is written from
        # dict_return, which only ever held the FIRST run's verdict.
        _rep_collapsed = [bool(dict_return.get('live_volume_collapsed', False))
                          if isinstance(dict_return, dict) else False]
        _rep_collapse_why = [("run 1: " + dict_return['collapse_reason'])
                             if isinstance(dict_return, dict) and dict_return.get('collapse_reason')
                             else None]
        for _irep in range(int(opts.mc_error_replicas)):
            # cold restart: drop the sample cache and reset per-parameter adaptation so
            # this replica is independent of the runs before it (AV cold-starts by
            # construction; the AC/GPU sampler needs the explicit reset; the portfolio and
            # the standalone GMM each need their own, below.  Any sampler matching none of
            # these reruns warm -- still an independent realization of the draws, just not
            # of the adaptation).
            sampler._rvs = {}
            # PORTFOLIO first.  Neither mcsamplerPortfolio, mcsamplerAdaptiveVolume nor
            # mcsamplerEnsemble defines reset_sampling (only the AC/GPU sampler does), so the
            # loop below was a no-op for them -- and because portfolio.integrate_log does NOT
            # call self.setup(), a portfolio replica reran with the PREVIOUS replica's adapted
            # grid and fitted GMM.  Such replicas share the very adaptation whose failure they
            # exist to detect, so their scatter understates the true MC error.  clear_warm_state
            # rebuilds every member from its original setup arguments.  (Standalone AV is
            # already cold: its integrate_log calls setup() itself.)
            # reset_adaptation() is the FULL reset: member proposals via clear_warm_state PLUS the
            # portfolio's own learned state (draw allocation, quality EMAs and their counts, probe
            # pointer, iteration counter, n_ess histories).  Clearing only the members leaves each
            # replica scheduling itself from what the previous replicas learned, so they are not
            # adaptation-independent and the between-replica scatter -- the whole quantity being
            # measured -- still understates the error.
            if hasattr(sampler, 'reset_adaptation'):
                sampler.reset_adaptation()
            elif hasattr(sampler, 'clear_warm_state'):
                sampler.clear_warm_state()
            elif hasattr(sampler, 'integrator'):
                # STANDALONE GMM (mcsamplerEnsemble) -- the one supported sampler with NONE of
                # the resets above, and the one that warm-starts hardest.  It carries the fit
                # forward by TWO routes, so both have to be cut or the replicas share the very
                # adaptation whose failure they exist to detect:
                #   1. its integrate() builds a fresh integrator and then deliberately
                #      TRANSFERS every fitted model from the previous self.integrator into it
                #      (the warm-start-survival path).  Dropping self.integrator is the whole
                #      reset -- integrate() rebuilds it from its own arguments, and with no
                #      previous integrator the transfer is skipped.
                #   2. the gmm_dict passed in pinned_params is ALIASED, not copied: the
                #      MonteCarloEnsemble integrator holds that very object and _train writes
                #      each refitted model back into it, so the driver's dict accumulates the
                #      previous run's proposals and would rewarm the next replica through
                #      kwargs even with self.integrator gone.
                # Only the ADAPTING groups are blanked, and for them None is the exact state
                # this dict had before the first run.  A non-adapting group (by default
                # (psi,phi_orb), seeded with the wide phase prior) is skipped by _train and so
                # is still pristine; blanking it would leave that group with no model at all
                # for the rest of the run -- a silent downgrade to uniform sampling, not a
                # cold start.  With --internal-rotate-phase that group does adapt, and its
                # seed was updated IN PLACE by the first run, so no pristine copy survives
                # anywhere to restore: it cold-starts from uniform, which is what the wide
                # single-component prior approximates anyway.
                sampler.integrator = None
                _gmm_dict_rep = pinned_params.get('gmm_dict', None)
                _gmm_adapt_rep = pinned_params.get('gmm_adapt', None)
                if not isinstance(_gmm_adapt_rep, dict):
                    _gmm_adapt_rep = {}
                if isinstance(_gmm_dict_rep, dict):
                    for _grp in list(_gmm_dict_rep.keys()):
                        if _gmm_adapt_rep.get(_grp, True):
                            _gmm_dict_rep[_grp] = None
            if hasattr(sampler, 'reset_sampling'):
                for _p in list(getattr(sampler, 'params_ordered', [])):
                    try:
                        sampler.reset_sampling(_p)
                    except Exception:
                        pass
            try:
                _r2 = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params)
            except Exception as _e_rep:
                print(" [mc error] replica {} failed ({}); continuing with {} results".format(_irep + 1, _e_rep, len(_rep_lnZ)))
                continue
            if (_r2 is None) or (_r2[0] is None):
                continue
            _res2, _var2, _neff2, _dd2 = _r2
            if not(opts.internal_use_lnL):
                _lr2 = numpy.log(_res2); _sig2 = numpy.sqrt(_var2)/_res2
            else:
                _lr2 = _res2; _sig2 = numpy.exp(_var2/2 - _lr2)
            _kh2, _sb2, _ne2, _unused = _extract_mc_diag(_dd2)
            if _sb2 is not None and numpy.isfinite(_sb2):
                _sig2 = max(float(_sig2), float(_sb2))
            _rep_lnZ.append(float(_lr2)); _rep_sig.append(float(_sig2)); _rep_neff.append(float(_neff2))
            _rep_rvs.append(sampler._rvs)
            _rep_fairdraw.append(bool(getattr(sampler, '_rvs_is_fairdraw', False)))
            _rep_collapsed.append(bool(_dd2.get('live_volume_collapsed', False))
                                  if isinstance(_dd2, dict) else False)
            if isinstance(_dd2, dict) and _dd2.get('collapse_reason'):
                _rep_collapse_why.append("replica {}: {}".format(_irep + 1, _dd2['collapse_reason']))
        # POOL the replicas rather than picking one, so the exported posterior is a draw from the
        # same mixture the reported evidence describes.
        #   Zhat = (1/K) sum_k (1/n_k) sum_i w_ki    ->    pooled weight  w_ki / (K n_k)
        # which is exactly the importance weight against the POOLED proposal density
        # q'_ki = q_ki * K * n_k  (pick a replica uniformly, then draw one of its n_k samples).
        # Folding the factor into log_joint_s_prior is therefore a statement of the real pooled
        # sampling density, not a fudge -- and it leaves every downstream weight computation
        # (which all form log_integrand + log_joint_prior - log_joint_s_prior) correct untouched.
        _pooled_rvs = _pool_replica_rvs(_rep_rvs, sampler, rep_lnZ=_rep_lnZ,
                                        already_resampled=_rep_fairdraw,
                                        use_lnL=rvs_integrand_is_lnL)
        # A POOLED RECORD IS NOT A FAIR DRAW, even when every block that went into it was.
        # _pool_replica_rvs gives block k weights summing to Z_k/K: equal WITHIN a block (each
        # block really is an equal-weight draw from its own posterior) but differing BETWEEN
        # blocks by exactly the replica evidences.  Leaving the marker set would make
        # ln_weights_for_posterior return zeros, and .dgrid and the proposal breadcrumb would
        # then mix the replicas by exported ROW COUNT instead of by evidence -- silently
        # discarding the disagreement the replicas were run to measure.  The reconstructed
        # per-row weights already encode it, so clear the marker and let them be read.
        #
        # Only when it actually pooled: every fallback path in _pool_replica_rvs returns one of
        # its INPUT records unchanged (too few replicas, no sampling-prior column, an exception),
        # and such a record is still the fair draw it arrived as.  Identity, not length, is the
        # reliable test for that.
        _did_pool = not any(_pooled_rvs is _r for _r in _rep_rvs)
        if _did_pool:
            # POOLED, not equal-weight.  The rows are still posterior-resampled wherever their
            # block was (so _rvs_is_fairdraw stays, and the .dslice safeguard keeps firing),
            # but the record as a whole is a mixture weighted by the replica evidences, so
            # ln_weights_for_posterior must read the reconstructed per-row weights.
            sampler._rvs_is_pooled = True
            sampler._rvs_is_fairdraw = any(_rep_fairdraw)
            # DELIBERATELY record-less, and this line is why it is deliberate.  The main
            # driver builds an _RvsRecord.pooled() here from per-replica records; the LISA
            # replica path does not collect them, so there is nothing honest to publish and
            # the weight route falls back to the flags above.  That fallback is correct --
            # but WITHOUT this line it would be correct only by accident: the record left on
            # the sampler describes the pre-pool columns, and it is declined solely because
            # `sampler._rvs` is about to become a different dict and _rvs_record_for compares
            # by IDENTITY.  Anything that later made the pooled dict reuse an input dict, or
            # added a samples() consumer here, would silently start reading a per-pass record
            # as if it described the mixture.  Clear it, so the absence is a statement.
            if _sampler_keeps_records(sampler):
                sampler.set_samples(None)
        # Did pooling FLATTEN any block?  That, not "is the record resampled", is what makes
        # the pooled Kish n_eff meaningless below -- a flattened block's rows carry its export
        # size rather than its integration quality.
        _blocks_flattened = bool(_did_pool and any(_rep_fairdraw))
        sampler._rvs = _pooled_rvs
        # The pooled export is a mixture over every replica in _rep_rvs, so its collapse
        # status is the OR over them: one collapsed member taints the pool.  Fold that back
        # into dict_return, which is what the status sidecar and the downstream reporting
        # read -- otherwise a healthy first run followed by a collapsed replica would export
        # the pooled posterior while recording "collapsed": false.
        if isinstance(dict_return, dict):
            _any_collapsed = any(_rep_collapsed)
            _why = [w for w in _rep_collapse_why if w]
            dict_return['live_volume_collapsed'] = bool(_any_collapsed)
            dict_return['n_replicas_pooled'] = int(len(_rep_lnZ))
            dict_return['n_replicas_collapsed'] = int(sum(1 for c in _rep_collapsed if c))
            if _any_collapsed:
                dict_return['collapse_reason'] = "; ".join(_why) if _why else "a pooled replica collapsed"
                print(" [mc error] *** LIVE VOLUME COLLAPSED in {} of {} pooled replicas ***".format(
                    dict_return['n_replicas_collapsed'], dict_return['n_replicas_pooled']))
                print(" [mc error] {}".format(dict_return['collapse_reason']))
                print(" [mc error] the POOLED posterior therefore contains degenerate samples.")
        # Re-apply the rejection gate to the POOLED verdict.  The early call above saw only
        # the first run, so without this a healthy first run followed by a collapsed replica
        # would export the pooled, collapsed result with --reject-collapsed-live-volume set.
        _reject_if_collapsed(dict_return, "pooled over {} replicas".format(len(_rep_lnZ)))
        if len(_rep_lnZ) > 1:
            _K = len(_rep_lnZ)
            _l = numpy.array(_rep_lnZ); _s = numpy.array(_rep_sig)
            _lref = numpy.max(_l)
            _Z = numpy.exp(_l - _lref)
            _Zbar = numpy.mean(_Z)
            _lnZ_comb = numpy.log(_Zbar) + _lref          # linear mean over replicas: unbiased in Z
            _sig_prop = float(numpy.sqrt(numpy.sum((_s*_Z)**2))/(_K*_Zbar))
            _sig_scatter = float(numpy.std(_l, ddof=1)/numpy.sqrt(_K))   # t_{K-1}: small-K quantiles are wider than Gaussian, hence the max() below
            _sig_comb = max(_sig_prop, _sig_scatter)
            print(" [mc error] combined {} replicas: lnZ {} -> {:.4f} (shift {:+.3f} vs first); sigma propagated {:.3f} / scatter {:.3f} -> {:.3f}; neff {} -> {:.1f}".format(
                _K, numpy.array2string(_l + lnL_offset, precision=3), float(_lnZ_comb + lnL_offset), float(_lnZ_comb - _rep_lnZ[0]),
                _sig_prop, _sig_scatter, _sig_comb, numpy.array2string(numpy.asarray(_rep_neff), precision=1), float(numpy.sum(_rep_neff))))
            log_res = float(_lnZ_comb)
            sqrt_var_over_res = _sig_comb
            # Report the POOLED n_eff, not the sum.  The sum claims the posterior carries the
            # combined effective sample size of K independent runs, which is only true if they
            # agree; when they disagree -- the case these replicas exist to detect -- the pooled
            # Kish n_eff is smaller, and that disagreement is exactly what should show up here.
            #
            # ...but NOT the Kish n_eff OF THE POOLED RECORD when that record is the fair-draw
            # export.  _pool_replica_rvs deliberately FLATTENS each block in that case (equal
            # weights within a block, summing to Z_k/K), and the Kish n_eff of piecewise-constant
            # weights is just the row count -- i.e. K*min(n_max, 1.5*eff_samp, 1.5*neff), the size
            # of the EXPORT, which says nothing about how well the integral converged.  NOTE this
            # driver has no --fairdraw-extrinsic-output-n-max: it caps the export at opts.n_eff
            # (igrand_fairdraw_samples_max), so the bogus figure would be K*n_eff, not 5K.
            #
            # Do the same computation one level up, where the quantities are still meaningful:
            # Kish over the BLOCKS, each carrying its own Z_k and its own n_eff,
            #
            #     neff_pooled = (sum_k Z_k)^2 / sum_k (Z_k^2 / neff_k)
            #
            # which has exactly the property the paragraph above asks for: it reduces to
            # sum_k neff_k when the replicas agree, and falls below it when they disagree --
            # the disagreement these replicas exist to detect.
            if _blocks_flattened:
                _l_rel = numpy.asarray(_rep_lnZ, dtype=float) - float(numpy.max(_rep_lnZ))
                _Zk = numpy.exp(_l_rel)
                _nk = numpy.asarray(_rep_neff, dtype=float)
                _ok = numpy.isfinite(_Zk) & numpy.isfinite(_nk) & (_nk > 0)
                _neff_pooled = (float(numpy.sum(_Zk[_ok]) ** 2 / numpy.sum(_Zk[_ok] ** 2 / _nk[_ok]))
                                if numpy.any(_ok) else None)
                _neff_how = 'block Kish over replicas (the export is fair-drawn)'
            else:
                _neff_pooled = _kish_neff_of_rvs(sampler._rvs)
                _neff_how = 'Kish over the pooled samples'
            neff = float(_neff_pooled) if _neff_pooled is not None else float(numpy.sum(_rep_neff))
            if _neff_pooled is not None:
                print(" [mc error] pooled posterior: {} samples, n_eff {:.1f} via {} (sum over replicas was {:.1f})".format(
                    len(numpy.atleast_1d(list(sampler._rvs.values())[0])) if sampler._rvs else 0,
                    float(_neff_pooled), _neff_how, float(numpy.sum(_rep_neff))))
            # keep the (res, var) pair consistent for any downstream reader
            if not(opts.internal_use_lnL):
                res = numpy.exp(log_res); var = (sqrt_var_over_res*res)**2
            else:
                res = log_res; var = 2*numpy.log(sqrt_var_over_res) + 2*log_res
    return res, var, neff, log_res, sqrt_var_over_res, dict_return


def _maybe_l0_rescue(sampler, res, var, neff, dict_return,
                     like_to_integrate, unpinned_params, pinned_params,
                     lnL_offset=0.0):
    """Run the L0 auto-rescue if this pass stalled -> (res, var, neff, dict_return).

    Returns its arguments unchanged when the rescue does not apply, so the call site is a
    single unconditional assignment.

    MUST BE CALLED BEFORE the `if not(res): raise` guard.  A degenerate early termination
    (mcsamplerPortfolio/AV returning (None,None,None,None) when the live volume never found
    finite in-volume samples) is the STRONGEST rescue trigger, not a reason to skip -- such a
    pass still populated _rvs, so the peak seed is available.  In the main driver that guard
    sits ~200 lines further down and the ordering is implicit; here it is immediately after
    integrate, so the ordering is stated and pinned by a test.

    `lnL_offset` is this event's manual_avoid_overflow_logarithm, used only to print absolute
    lnZ values.  It is a local of the caller in both analyze_event variants, hence a parameter.
    """
    # Reset per event before ANY early return.  The sampler is reused: a rejected rescue on
    # one event must not suppress saving a later healthy event that needs no rescue at all.
    sampler._av_state_reuse_safe = True
    # APPLICABILITY FIRST, then n_eff.  The main driver evaluates
    #     _neff_val = None if neff is None else float(sampler.identity_convert(neff))
    # BEFORE its guard, which is safe there only by luck: identity_convert comes from
    # MCSamplerGeneric, and RIFT.integrators.mcsampler.MCSampler -- the object this driver
    # keeps for --sampler-method adaptive_cartesian -- does NOT inherit it.  Evaluating it
    # unconditionally therefore raises AttributeError on EVERY adaptive_cartesian event, at
    # the end of a completed integration and before --output-file is written, losing the
    # whole point's compute.  The rescue is AV/portfolio-only regardless, so nothing is lost
    # by asking whether it applies before touching the sampler's conversion helpers.
    #
    # DELIBERATE DIVERGENCE from the main driver, which has the same latent defect on the
    # line above its own guard and should take the same reordering.
    if not (opts.sampler_method in ('AV', 'portfolio') and opts.sampler_warmstart_retry_neff
            and hasattr(sampler, 'bootstrap_from_samples')):
        return res, var, neff, dict_return
    # A DEGENERATE EARLY TERMINATION (neff None) counts as below threshold, not as "skip".
    _neff_val = None if neff is None else float(sampler.identity_convert(neff))
    _needs_l0_rescue = (_neff_val is None) or (_neff_val < float(opts.sampler_warmstart_retry_neff or 0))
    if not _needs_l0_rescue:
        return res, var, neff, dict_return

    # Cold state to fall back on, captured only once the warm pass is actually about to run.
    # `None` means nothing has been disturbed yet, so the handler must not "restore".
    _cold_state_l0 = None
    try:
        # SEED FROM THE POINTS THE PASS RETAINED, not from what survived the fair draw.
        # sampler._rvs has by now been REBOUND to a fair-draw subset taken WITH REPLACEMENT --
        # a resample built for EXPORT.  On the collapsed pass this rescue exists for the
        # effective sample size is ~1, so _rvs can be a single row, or a handful several of
        # which are the same point twice.  The live set held a thousand.
        _res_l0 = _warm_seed_reserve_for(sampler)
        if _res_l0 is not None:
            _cols = np.asarray(_res_l0['X'], dtype=float)
            _lnv = np.asarray(_res_l0['lnL'], dtype=float).ravel()
            print("  [L0 auto-rescue] seeding from {} retained sample(s) of {} (fair draw left {} in _rvs)".format(
                len(_lnv), _res_l0.get('n_retained', '?'),
                len(np.asarray(sampler.identity_convert(sampler._rvs['log_integrand'])).ravel())
                if 'log_integrand' in sampler._rvs else '?'))
        else:
            _lnkey = 'log_integrand' if 'log_integrand' in sampler._rvs else ('integrand' if 'integrand' in sampler._rvs else None)
            _lnv = np.asarray(sampler.identity_convert(sampler._rvs[_lnkey]), dtype=float).ravel() if _lnkey else np.array([])
            _cols = (np.vstack([np.asarray(sampler.identity_convert(sampler._rvs[p]), dtype=float).ravel()
                                for p in sampler.params_ordered]).T if _lnv.size else np.zeros((0, len(sampler.params_ordered))))
        if _lnv.size >= 1 and np.any(np.isfinite(_lnv)):
            # RANK, not count, decides whether this seed can define a live volume.  A 2-to-5
            # point seed passes a count test and is still rank-deficient in 6 adaptive
            # dimensions, so the warm start contracts onto a degenerate subspace and reports a
            # healthy n_eff over a sliver of the support.  build_warm_seed applies the rank
            # test through the SAME seed_affine_rank the grid builder uses, and puffs to full
            # rank when it is short.
            _ax_l0, _lo_l0, _hi_l0 = _warm_seed_geometry(sampler)
            _seed, _seed_info = mcsamplerAdaptiveVolume.build_warm_seed(
                _cols, _lnv, _lo_l0, _hi_l0, _ax_l0,
                deltalnL=opts.sampler_sequential_warmstart_deltalnL,
                puff_scale=opts.sampler_l0_rescue_puff_scale,
                puff_width_frac=opts.sampler_l0_rescue_puff_width_frac,
                puff_factor=opts.sampler_l0_rescue_puff_factor)
            print("  [L0 auto-rescue] cold n_eff {} < {}; re-running warm from this point's peak ({} pts)".format(
                "DEGENERATE (early termination)" if _neff_val is None else "{:.1f}".format(_neff_val),
                opts.sampler_warmstart_retry_neff, len(_seed)))
            if _seed_info['puffed']:
                print("  [L0 auto-rescue] seed of {} point(s) had affine rank {}/{}: PUFFED to rank"
                      " {}/{} with {} points ({} scale, x{:g}), keeping the original point(s)".format(
                          _seed_info['n_core'], _seed_info['rank_core'], _seed_info['dim'],
                          _seed_info['rank_final'], _seed_info['dim'], _seed_info['n_puff'],
                          _seed_info['puff_scale'], opts.sampler_l0_rescue_puff_factor))
                if _seed_info['rank_final'] < _seed_info['dim']:
                    print("  [L0 auto-rescue] *** the puffed seed is STILL rank-deficient"
                          " ({}/{}); the warm pass will be reported as collapsed.".format(
                              _seed_info['rank_final'], _seed_info['dim']))
            # The warm pass is an estimate over TRUNCATED support: the seeded box provably
            # contains the peak the cold pass found, and says nothing about what that pass did
            # not reach, so it is biased low by any missed mode.  The rescue still runs, because
            # it exists to fix the high-SNR n_eff lottery and removing it by default would be a
            # certain production regression traded against a possible bias.  What the gate below
            # changes is only the case where there is POSITIVE EVIDENCE of lost mass.
            #
            # SNAPSHOT, not an alias: integrate_log repopulates sampler._rvs IN PLACE, so
            # `_cold_rvs = sampler._rvs` would be holding the warm samples by the time the
            # restore ran -- i.e. the reject path would report the cold lnZ while exporting the
            # warm cloud, exactly what it exists to prevent.
            _cold_rvs = dict(sampler._rvs)
            # Snapshot the RESERVE for the same reason and at the same moment: the warm pass's
            # integrate_log clears and rewrites it, so reading it after the fact would compare
            # the warm pass against itself.
            _cold_reserve_l0 = getattr(sampler, '_warm_seed_reserve', None)
            _cold_lnZ, _cold_src = _lnZ_of_reserve_or_rvs(sampler, _cold_rvs,
                                                          reserve=_cold_reserve_l0)
            # dict_return too: khat, block scatter, ESS and the confidence interval all read it,
            # so keeping the warm pass's diagnostics beside a restored cold result would describe
            # a run we did not report.  And the reserve and the fair-draw marker, one level out.
            _cold_state_l0 = _snapshot_pass_state(sampler, res, var, neff, dict_return,
                                                  rvs=_cold_rvs)
            sampler.bootstrap_from_samples(_seed, cover_frac=0.0)
            res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params)
            _warm_lnZ, _warm_src = _lnZ_of_reserve_or_rvs(sampler, sampler._rvs)
            # BOTH SIDES FROM THE SAME READING, or the difference is not a difference.  A
            # fair-drawn lnZ sits ~log(n_retained/eff_samp) above a retained-set one, so a mixed
            # comparison manufactures a gap of several nats in whichever direction the mismatch
            # happens to fall.  If the two passes did not produce the same kind of estimate, read
            # BOTH from _rvs -- the old behaviour, at least self-consistent.
            if _cold_src != _warm_src:
                print("  [L0 auto-rescue] lnZ provenance differs (cold={}, warm={});"
                      " re-reading both from the fair-draw record so the comparison is"
                      " like-for-like.".format(_cold_src, _warm_src))
                _cold_lnZ = _lnZ_of_rvs(_cold_rvs, already_pooled=False)
                _warm_lnZ = _lnZ_of_rvs(sampler._rvs, already_pooled=False)
                _cold_src = _warm_src = 'fairdraw'
            _evidence_of_loss = (
                (_cold_lnZ is not None) and (_warm_lnZ is not None)
                and numpy.isfinite(_cold_lnZ) and numpy.isfinite(_warm_lnZ)
                and (_cold_lnZ - _warm_lnZ) > float(opts.sampler_l0_rescue_reject_dlnZ))
            if _evidence_of_loss:
                print("  [L0 auto-rescue] *** REJECTING the warm pass *** its lnZ {:.3f} is"
                      " {:.3f} nats BELOW the full-support cold pass ({:.3f}), which is evidence"
                      " the seed missed mass the cold pass reached.".format(
                          _warm_lnZ + lnL_offset, _cold_lnZ - _warm_lnZ, _cold_lnZ + lnL_offset))
                if opts.sampler_l0_rescue_accept_truncated:
                    print("  [L0 auto-rescue] --sampler-l0-rescue-accept-truncated set:"
                          " reporting the warm pass anyway (may be biased LOW).")
                else:
                    print("  [L0 auto-rescue] keeping the COLD (full-support) result; its n_eff"
                          " is lower but it is not missing mass.  A portfolio avoids this"
                          " trade entirely -- its GMM member carries a defensive component.")
                    # The RESERVE goes back too.  Once --sampler-sequential-warmstart is
                    # ported here, omitting this would seed the next intrinsic point from the
                    # warm cloud this gate just rejected: _warm_seed_reserve_for would return
                    # the warm pass's record while _rvs, the estimate and the diagnostics all
                    # describe the cold one.  Nothing reads it in this driver today.
                    res, var, neff, dict_return = _restore_pass_state(sampler, _cold_state_l0)
                    sampler._av_state_reuse_safe = False
            _clear_warm_state(sampler)
    except Exception as _e_l0:
        # "skipped" is only true if the warm pass never started.  If it raised PARTWAY THROUGH
        # sampler.integrate(), the assignment never completed, so res/var/neff/dict_return still
        # hold the COLD pass -- while sampler._rvs was repopulated in place and now holds the
        # WARM samples.  Reporting cold k-hat / ESS / lnZ beside a warm export describes a run
        # that was never made, and it did so silently for a whole campaign.
        print("  [L0 auto-rescue] *** FAILED *** (", _e_l0, ")")
        import traceback as _tb_l0
        _tb_l0.print_exc()
        if _cold_state_l0 is not None:
            print("  [L0 auto-rescue] the warm pass may already have replaced the stored"
                  " samples; restoring the COLD pass so the reported diagnostics and the"
                  " exported samples describe the same integral.")
            res, var, neff, dict_return = _restore_pass_state(sampler, _cold_state_l0)
            sampler._av_state_reuse_safe = False
        _clear_warm_state(sampler)
    return res, var, neff, dict_return


def resample_samples_LISA(my_samples, rholms, cross_terms, right_ascension, declination, P, modes, reference_distance):
  """This function takes in extrinsic samples and for each sample samples a time shift. This is done by generating a likelihood time series at an extrinsic sample and then weighted sampling in time."""
  # How many time samples? Same as the extrinsic samples being passed
  n_samples = len(my_samples['longitude'])
  print(" Time resampling size : {} ".format(n_samples))

  # Hardcoded time sampling limits, should change it in future
  low_t_lim, high_t_lim = -5, 5
  delta_t = (1.0 / opts.srate_resample_time_marginalization
             if opts.srate_resample_time_marginalization else 0.001)
  
  # t_ref_wind is defined as data_window_integration_half

  # This is the time series over which the lnL series is defined.
  tvals = xpy_default.linspace(-t_ref_wind, t_ref_wind - P.deltaT,int(t_ref_wind*2/P.deltaT))  # the array should have P.deltaT as deltaT between entries and for that you need to subtract P.deltaT from the final point.
  
  # Generating likelihood time series for these samples
  P.phi =  identity_convert_togpu(my_samples['right_ascension'])  # cast to float
  P.theta =   identity_convert_togpu(my_samples['declination'])
  P.tref = float(fiducial_epoch) # defined earlier as lisa-reference-time
  P.phiref = identity_convert_togpu(my_samples['coa_phase'])
  P.incl = identity_convert_togpu(my_samples['inclination'])
  P.psi = identity_convert_togpu(my_samples['psi'])
  P.dist = identity_convert_togpu(my_samples['distance']* 1.e6 * lalsimutils.lsu_PC) # luminosity distance
  
  if not(cupy_success):
    # convert objects
    for name in ['phi','theta', 'phiref','incl', 'psi','dist']:
      setattr(P,name, getattr(P,name).astype(float) )
  
  # plugging the selected samples
  lnLt = factored_likelihood_LISA.FactoredLogLikelihoodAlignedSpinLISA(rholms, cross_terms, declination, right_ascension, P.psi, P.incl, P.phiref, P.dist, modes, reference_distance, return_lnLt = True)
  
  # axis 0  is time axis 1 is each point, transposing it
  lnLt = lnLt.T 
  lnLt = identity_convert(lnLt) # back to CPU.  Note we have removed offsets 
  if opts.zero_likelihood:
    lnLt =np.zeros(lnLt.shape)  # zero likelihood
  lnLt_norm = scipy.special.logsumexp(lnLt,axis=-1)
  tvals = identity_convert(tvals) # back to CPU
  
  # Draw continuously when sub-sample interpolation requested it; retain the
  # historical fixed 1 ms LISA export grid as the explicit compatibility path.
  t_out = np.zeros(n_samples)
  lnL_out = np.zeros(n_samples)
  # identifiy max lnL(t) point, and then interpolate around that point to avoid Nan weights.
  index_max_lnLt = np.argmax(lnLt[0]).flatten()
  tval_at_max = tvals[index_max_lnLt]
  print(f"Identified max lnL(t) as {tval_at_max}s. Interpolating around that point.  Sampling for time from {tval_at_max + low_t_lim}s to {tval_at_max + high_t_lim}s with deltaT {delta_t}s.")
  
  # tvals over which we will resample
  interp_tvals  = xpy_default.arange(low_t_lim + tval_at_max, high_t_lim + tval_at_max, delta_t)
  indx_list = np.arange(len(interp_tvals))
  # Loop over each extrinsic sample
  for indx in np.arange(n_samples):
    interp = scipy.interpolate.interp1d(tvals, lnLt[indx] - lnLt_norm[indx], kind='cubic')
    # interpolating
    interp_lnLt = interp(interp_tvals)
    # weighted resampling
    indx_choose = np.random.choice(indx_list, p=np.exp(interp_lnLt)/np.sum(np.exp(interp_lnLt)))
    t_out[indx] = interp_tvals[indx_choose]
    lnL_out[indx] = interp_lnLt[indx_choose]
  print(t_out)
  # saving 
  my_samples['t_ref'] = fiducial_epoch + t_out # add sample time jitter from reweighting to samples
  my_samples["lnL_raw"] = lnL_out   # export likelihoods (equivalent to SNR). Note needs downstream code filters to catch this and put it somewhere.
  return my_samples


def resample_samples(my_samples,
                     lookupNKDict=None, rholmArrayDict=None, ctUArrayDict=None, ctVArrayDict=None,epochDict=None): # uses LOTS of global variables, don't pass them all
  # will look a LOT like the likelihood function definitions, unfortunately
  if not opts.vectorized:
    raise Exception( ' resample_samples currently only for vectorized ')
  if opts.distance_marginalization:
    raise Exception( ' resample_samples currently only final extrinsic samples; you should NOT have distance marginalization on ')

  n_samples = len(my_samples['longitude'])
  print(" Time resampling size : {} ".format(n_samples))


  # if we are using GPU-based generation this is ok; if itis mcsampler, we will have a lot of 'object' casts to fix, arg
  
  tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally
  P.phi =  identity_convert_togpu(my_samples['right_ascension'])  # cast to float
  P.theta =   identity_convert_togpu(my_samples['declination'])
  P.tref = float(fiducial_epoch)
  P.phiref = identity_convert_togpu(my_samples['coa_phase'])
  P.incl = identity_convert_togpu(my_samples['inclination'])
  P.psi = identity_convert_togpu(my_samples['psi'])
  P.dist = identity_convert_togpu(my_samples['distance']* 1.e6 * lalsimutils.lsu_PC) # luminosity distance
  if not(cupy_success):
    # convert objects
    for name in ['phi','theta', 'phiref','incl', 'psi','dist']:
      setattr(P,name, getattr(P,name).astype(float) )

  lnLt = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,
                        P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default,return_lnLt=True)

  lnLt = identity_convert(lnLt) # back to CPU.  Note we have removed offsets 
  if opts.zero_likelihood:
    lnLt =np.zeros(lnLt.shape)  # zero likelihood
  lnLt_norm = scipy.special.logsumexp(lnLt,axis=-1)
  tvals = identity_convert(tvals) # back to CPU
#  print(lnLt.shape, lnLt_norm.shape,tvals.shape)
  # Match the main driver's sub-sample export contract on this legacy
  # ground-based path retained in the LISA executable.
  t_out = np.zeros(n_samples)
  lnL_out = np.zeros(n_samples)
  indx_list =np.arange(len(tvals))
  for indx in np.arange(n_samples):
    indx_choose = np.random.choice(indx_list, p=np.exp(lnLt[indx] - lnLt_norm[indx]))
    t_out[indx] = tvals[indx_choose]
    lnL_out[indx] = lnLt[indx][indx_choose]
#    print(' Resampled time offset {} '.format(t_out[indx])) #, lnLt[indx]-lnLt_norm[indx])
    
  my_samples['t_ref'] = fiducial_epoch+t_out # add sample time jitter from reweighting to samples
  my_samples["lnL_raw"] = lnL_out   # export likelihoods (equivalent to SNR). Note needs downstream code filters to catch this and put it somewhere.
  
  return my_samples

#
# Main analysis functions for LIGO (analyze_event) and LISA (analyze_event_LISA). This precomputes terms (eq 23 10.1103/PhysRevD.92.023002) and then based on your settings\
# creates a likelihood function which is then passed to a sampler.
#

def analyze_event_LISA(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec):
    print("\n###########################################################################################\nPrecomputing\n###########################################################################################")
    nEvals=0
    # PROVENANCE RESET, ON ENTRY, BEFORE ANYTHING CAN FAIL.
    #
    # _rvs_is_pooled describes the record this call is about to build, and it is set by THIS
    # function (the replica block) rather than by a sampler, so no sampler-side per-pass reset
    # can clear it.  Clearing it only on the normal return is not enough: _reject_if_collapsed
    # RAISES after pooling, the caller's `except Exception` swallows that and moves to the next
    # event, and the marker survives.  The next ordinary fair draw is then read as "pooled",
    # _rvs_is_equal_weight goes False, and any consumer that weights rows applies importance
    # weights to rows that already carry them -- the w^2 defect, resurrected on the event after
    # any failure.
    #
    # On ENTRY rather than in a `finally`: entry is reached on every call by construction, and
    # it leaves the state correct even for a caller that never returns normally at all.
    sampler._rvs_is_pooled = False
    P = P_list[indx_event]
    # if pin-distance-to-sim, change the distance prior accordingly
    if opts.pin_distance_to_sim:
      pinned_params['distance']=P.dist/(1.e6 * lal.PC_SI)

    # Call external likleihood preparation if any
    if supplemental_ln_likelhood_prep:
      supplemental_ln_likelhood_prep(P=P,config=supplemental_ln_likelhood_parsed_ini)
    
    print("Generating modes for parameters:")
    P.print_params_lisa()
    # Generate modes
    if opts.modes:
      modes=np.array(eval(opts.modes))
      print(f"modes = \n{list(modes)}")
    P.deltaF=deltaF # why do we need this? why is it becoming None? 
    hlms_FD = lisa_lalsimutils_compat.hlmoff_for_LISA(P, opts.l_max, modes)
    fNyq = 0.5/P.deltaT
    modes = np.array(list(hlms_FD.keys()))
    reference_distance = P.dist
    print(f"Reference distance = {reference_distance/lal.PC_SI/1e9} Gpc, 1/deltaF = {1/P.deltaF} s, deltaT = {P.deltaT} s")

    # Are we varying skylocation? If we are, read it from the grid and pass it to precompute. If not, take in the value as is passed as an argument to ILE.
    if not(opts.lisa_fixed_sky):
        lisa_sky_beta = P.theta
        lisa_sky_lamda = P.phi
    if opts.lisa_fixed_sky:
        lisa_sky_lamda = float(opts.ecliptic_longitude)
        lisa_sky_beta = float(opts.ecliptic_latitude)

    # Precompute
    # fNyq is the max resolvable frequency for a waveform. It is 0.5/deltaT. RIFT needs deltaT, deltaF for waveform generation (information present in P) and for integration it needs fmax (fmax <= fNyq)
    print(f"Sky location lambda = {lisa_sky_lamda}, sky location beta = {lisa_sky_beta}")

    rholms_intp, cross_terms, cross_terms_V,  rholms,  guess_snr, rest = factored_likelihood_LISA.PrecomputeAlignedSpinLISA(opts.lisa_reference_time, opts.lisa_reference_frequency, opts.data_integration_window_half, hlms_FD, None, data_dict, psd_dict, flow_ifo_dict["A"], fNyq, fmax, P.deltaT, lisa_sky_beta, lisa_sky_lamda, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0)
    
    # reset to default.  Should not be needed, but weird python scoping error
    manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default 

    def likelihood_function(right_ascension, declination, phi_orb, inclination, psi, distance):
      P.phi = xpy_asarray_already(right_ascension) 
      P.theta = xpy_asarray_already(declination)
      P.phiref = xpy_asarray_already(phi_orb)
      if opts.inclination_cosine_sampler:
        P.incl = xpy_default.arccos(xpy_asarray_already(inclination))
      else:
        P.incl = xpy_asarray_already(inclination)
      P.psi = xpy_asarray_already(psi)
      P.dist = xpy_asarray_already(distance* 1.e6 * lalsimutils.lsu_PC) # luminosity distance

      # rotate phase if needed
      # make copies of arrays
      # Sampling assumes P.phiref == phi+ \psi ,   P.psi == phi - psi
      if opts.internal_rotate_phase:
        phi_orb_true = (P.phiref + P.psi)/2.
        psi_true = (P.phiref - P.psi)/2.
        P.psi= psi_true
        P.phiref = phi_orb_true
      #print(f"Sky location lambda = {lisa_sky_lamda}, sky location beta = {lisa_sky_beta}")


      lnL = factored_likelihood_LISA.FactoredLogLikelihoodAlignedSpinLISA(rholms, cross_terms, lisa_sky_beta, lisa_sky_lamda, P.psi, P.incl, P.phiref, P.dist, modes, reference_distance)
      return lnL
      #return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm))
    
    if opts.sampler_method == "adaptive_cartesian_gpu":
      # reset sampling parameter as needed (distance, inclination but not sky location)
      # distance (and inclination) can conceivably correlate with mass, and we don't want to truncate in distance/inclination prematurely from history effects
      # method ONLY implemented for adaptcart!
      if 'distance' in sampler.params:
        sampler.reset_sampling('distance')
      sampler.reset_sampling('inclination')
    elif use_gmm_args:   # standalone GMM or a portfolio with a GMM member (gmm_dict exists in both)
      if 'distance' in sampler.params:
         pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination'])
         if  pair_d_incl in gmm_dict:
            gmm_dict[pair_d_incl] = None # reset inclination/distance sampling
      single_incl = sampler_param_tuple(sampler, ['inclination'])
      if (opts.distance_marginalization) and single_incl in gmm_dict:
        gmm_dict[single_incl] = None # reset inclination sampling
    print("\n###########################################################################################\nIntegrating\n###########################################################################################")
    # Integrate
    args = likelihood_function.__code__.co_varnames[:likelihood_function.__code__.co_argcount]
    print( " --------> Arguments ", args)
    # if exactly zero likelihood function
    like_to_integrate = likelihood_function
    if opts.zero_likelihood:
              if opts.internal_use_lnL:
                   like_to_integrate = zero_like
              else:
                   like_to_integrate = unit_like
    if oracleRS:
         print("  ORACLE  - seeding sampling  with" , oracleRS.params_ordered)
         if hasattr(sampler, 'update_sampling_prior'):
             rvs_train = {}
             _, _, rv_oracle = oracleRS.draw_simplified(opts.n_chunk)
             for indx,p  in enumerate(sampler.params_ordered):
               rvs_train[p] = rv_oracle[:,indx]
             # train with equal weight - no likelihood information
             lnL_oracles  = np.zeros(opts.n_chunk)
             sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level)

    _maybe_load_av_state(sampler)
    _maybe_enable_anisotropic_bins(sampler)
    res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params)

    # L0 auto-rescue: on a very sharply-peaked (high-amplitude) point a cold AV can stall at
    # n_eff ~ 1 because it never draws near the tiny peak.  If so, seed a SECOND pass from this
    # same point's own highest-likelihood samples and re-run.  Opt-in via
    # --sampler-warmstart-retry-neff.  MUST run BEFORE the not(res) guard below: a degenerate
    # early termination returns (None,None,None,None) and is the strongest rescue trigger, so
    # raising on it first would skip exactly the case the rescue exists for.
    res, var, neff, dict_return = _maybe_l0_rescue(
        sampler, res, var, neff, dict_return,
        like_to_integrate, unpinned_params, pinned_params,
        lnL_offset=manual_avoid_overflow_logarithm)

    if not(res): # no resut
      raise ValueError(" No integral result returned")

    if not(opts.internal_use_lnL):
      log_res = numpy.log(res)
      sqrt_var_over_res =  numpy.sqrt(var)/res
    else:
      log_res = res
      sqrt_var_over_res =  numpy.exp(var/2 - log_res)

    # MC-error diagnostics, the collapse report/gate, and replica replication+pooling.  This
    # OWNS both collapse-gate calls (first run and pooled verdict) -- do not add a separate
    # _report_and_gate_collapse call here, or the pooled one gets bypassed.  It needs
    # log_res/sqrt_var_over_res, hence its position after they are computed; the main driver
    # has the same ordering inline.
    res, var, neff, log_res, sqrt_var_over_res, dict_return = _maybe_replicate_for_mc_error(
        sampler, res, var, neff, dict_return, log_res, sqrt_var_over_res,
        like_to_integrate, unpinned_params, pinned_params,
        lnL_offset=manual_avoid_overflow_logarithm)

    # Report results
    if opts.output_file:
        fname_output_txt = opts.output_file +"_"+str(indx_event)+"_" + ".dat"
        m1 =P.m1/lal.MSUN_SI
        m2 =P.m2/lal.MSUN_SI
        if opts.sim_xml: 
            event_id = opts.event
        else:
            event_id = -1
        if opts.event == None:
            event_id = -1
        # Hyperpipeline ASCII output (opt-in via RIFT_HYPERPIPELINE_FORMAT).  The
        # DAG/CEPP join (util_CleanILE_hyperpipeline.py) consumes self-describing
        # header-bearing shards with lnL/sigma_lnL as columns 0/1.  The legacy
        # headerless savetxt below cannot be read by that path.
        from RIFT.misc import hyperpipeline_io as _hpio
        if _hpio.is_active():
            _cols = _hpio.build_column_list(use_sky=True)
            _vals = {
                "lnL": log_res + manual_avoid_overflow_logarithm,
                "sigma_lnL": sqrt_var_over_res,
                "m1": m1, "m2": m2,
                "a1x": P.s1x, "a1y": P.s1y, "a1z": P.s1z,
                "a2x": P.s2x, "a2y": P.s2y, "a2z": P.s2z,
                "ecliptic_longitude": lisa_sky_lamda,
                "ecliptic_latitude": lisa_sky_beta,
            }
            _hpio.write_row(fname_output_txt, _cols, [_vals.get(c, 0.0) for c in _cols])
        else:
            # Current response only applicable to quasicircular MBHB signals. Save sky location even if not varying.
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, lisa_sky_lamda, lisa_sky_beta, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res, sampler.ntotal, neff ]]))

    # Comprehensive output (not yet provided)
    # Convert declination, inclination  parameters in sampler if needed
    if opts.save_samples and opts.output_file:
      import copy
      samples = copy.deepcopy(sampler._rvs)  # deep copy: avoid modifying structures and  having side effect on integrator, which loops over keys Expensive!
      # A POOLED record (--mc-error-replicas) is weighted BETWEEN blocks by the replica
      # evidences, and nothing below preserves those weights: convert it to an equal-weight
      # draw BEFORE anything consumes it -- including resample_samples_LISA, which picks a time
      # per row and so assumes the rows already are the posterior.  A no-op otherwise.
      samples = _export_rvs_equal_weight(samples, sampler, use_lnL=rvs_integrand_is_lnL)
      # Insert reference distance if it was marginalized over
      if "distance" not in samples:
        # Not distance output is the same as internal calculations: in *Mpc*
        samples["distance"] = np.full_like(
          samples["psi"],
          factored_likelihood.distMpcRef,  #*1e6*lal.PC_SI,
          )
      # Insert reference phase if it was marginalized over
      if "phi_orb" not in samples:
        samples["phi_orb"] = np.full_like(samples["psi"], 0.)
      if opts.inclination_cosine_sampler:
        samples["inclination"] = numpy.arccos(samples["inclination"].astype(numpy.float64))
      if opts.declination_cosine_sampler:
        samples["declination"] = numpy.pi/2 - numpy.arccos(samples["declination"].astype(numpy.float64))
      xmldoc = ligolw.Document()
      xmldoc.appendChild(ligolw.LIGO_LW())
      process.register_to_xmldoc(xmldoc, sys.argv[0], opts.__dict__)
      if not(opts.resample_time_marginalization):
        if not opts.time_marginalization:
            samples["t_ref"] += float(fiducial_epoch)
        else:
            samples["t_ref"] = float(fiducial_epoch)*numpy.ones(len(samples["psi"]))
      # rotate sky to recover physical coordinates.  Note on CPU
      if opts.internal_sky_network_coordinates:
        tmp_th, tmp_ph = my_rotation_cpu(np.pi/2 - samples['declination'],samples['right_ascension'])
        samples['declination'] = np.pi/2 - tmp_th
        samples['right_ascension'] = np.mod(tmp_ph, 2*np.pi)

      # recover true phase, polarization samples
      if opts.internal_rotate_phase:
        phi_orb_true = np.mod((samples['phi_orb'] + samples['psi'])/2., 2*np.pi)
        psi_true = np.mod((samples['phi_orb'] - samples['psi'])/2., 2*np.pi)   # keep as 2 pi range, to be consistent with past work
        samples['psi']= psi_true
        samples['phi_orb'] = phi_orb_true
      samples["polarization"] = samples["psi"]
      samples["coa_phase"] = samples["phi_orb"]
      if ("declination", "right_ascension") in sampler.params:
            samples["latitude"], samples["longitude"] = samples[("declination", "right_ascension")]
      else:
            samples["latitude"] = samples["declination"]
            samples["longitude"] = samples["right_ascension"]
      if "log_integrand" in samples:
        samples["loglikelihood"] = samples["log_integrand"] +  manual_avoid_overflow_logarithm
      else:
        samples["loglikelihood"] = numpy.log(samples["integrand"]) + manual_avoid_overflow_logarithm  # export with consistent offset
      if not opts.rom_integrate_intrinsic:
            # ILE mode: insert fixed model parameters
            samples["mass1"] = numpy.ones(samples["psi"].shape)*m1 # opts.mass1
            samples["mass2"] = numpy.ones(samples["psi"].shape)*m2 # opts.mass2
            samples["spin1x"] =numpy.ones(samples["psi"].shape)*P.s1x
            samples["spin1y"] =numpy.ones(samples["psi"].shape)*P.s1y
            samples["spin1z"] =numpy.ones(samples["psi"].shape)*P.s1z
            samples["spin2x"] =numpy.ones(samples["psi"].shape)*P.s2x
            samples["spin2y"] =numpy.ones(samples["psi"].shape)*P.s2y
            samples["spin2z"] =numpy.ones(samples["psi"].shape)*P.s2z
            samples["alpha4"] =numpy.ones(samples["psi"].shape)*P.eccentricity
            samples["alpha5"] =numpy.ones(samples["psi"].shape)*P.lambda1
            samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2
            # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so
            samples["alpha2"] = numpy.zeros(samples["psi"].shape)
            samples["alpha3"] = numpy.zeros(samples["psi"].shape)

      if opts.resample_time_marginalization:
        samples = resample_samples_LISA(samples, rholms, cross_terms, lisa_sky_lamda, lisa_sky_beta, P, modes, reference_distance)
        samples["loglikelihood" ] = samples["lnL_raw"]  # export the non-time-marginalized likelihood, if we are in the final stages
#        print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref']))
      xmlutils.append_samples_to_xmldoc(xmldoc, samples)
        # Extra metadata
      dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal}
#      if 'distance' in pinned_params:
#            dict_out['distance'] = pinned_params["distance"]
      converged_result = False
      if "convergence_test_results" in dict_return:
        converged_result = dict_return["convergence_test_results"]["normal_integral"]
      xmlutils.append_likelihood_result_to_xmldoc(xmldoc, log_res+manual_avoid_overflow_logarithm, neff=neff, converged=converged_result, **dict_out)
      fname_output_xml = opts.output_file +"_"+str(indx_event)+"_" + ".xml.gz"
      utils.write_filename(xmldoc, fname_output_xml, compress="gz")




    if opts.maximize_only and opts.output_file:
      # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization)
      if "integrand" in sampler._rvs:
        indx_guess = numpy.argmax(sampler._rvs["integrand"])   # start search near maximum-likelihood point. (WARNING: can be very close by)
      else:
        indx_guess = numpy.argmax(sampler._rvs["log_integrand"])
      P.radec=True
      P.phi = sampler._rvs["right_ascension"][indx_guess]
      P.theta = sampler._rvs["declination"][indx_guess]
      P.phiref = sampler._rvs["phi_orb"][indx_guess]
      P.incl = sampler._rvs["inclination"][indx_guess]
      P.psi = sampler._rvs["psi"][indx_guess]
      P.dist = sampler._rvs["distance"][indx_guess]*1e6*lal.PC_SI
      print( " ---- Best extrinsic paramers in MC   ---- ")
      P.print_params()
      lalsimutils.ChooseWaveformParams_array_to_xml([P],"notime_raw_maxpt_"+opts.output_file) # best point, not including time

      import scipy.optimize
      def fn_scaled(x):
          P.phi = float(x[0]*2*numpy.pi) # right ascension
          P.theta = float((x[1])*numpy.pi) # declination, really polar angle. NOT zero on equator
          P.tref = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) # ref. time (rel to epoch for data taking)
          P.phiref = float(x[3]*numpy.pi*2) # ref. orbital phase
          P.incl = float(x[4]*numpy.pi) # inclination
          P.psi = float(x[5]*numpy.pi) # polarization angle
          P.dist = x[6]*dmax* 1.e6 * lalsimutils.lsu_PC # luminosity distance
    
          return -1.0* factored_likelihood.FactoredLogLikelihood(
                    P, rholms,rholms_intp, cross_terms, cross_terms_V, opts.l_max)

    # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization)
      x0 =numpy.array( [ \
       sampler._rvs["right_ascension"][indx_guess]/(2*numpy.pi) , \
       (sampler._rvs["declination"][indx_guess]/numpy.pi),  \
       0.5, \
       (sampler._rvs["phi_orb"][indx_guess]/(2*numpy.pi)), \
       sampler._rvs["inclination"][indx_guess]/(numpy.pi),\
       sampler._rvs["psi"][indx_guess]/numpy.pi,\
       sampler._rvs["distance"][indx_guess]/dmax\
            ],dtype=numpy.float128)
      x0 = numpy.fmod(x0,numpy.ones(len(x0)))  # had BETTER be defined on this range!
      # Pick the best starting time. BRUTE FORCE METHOD: use grid
      def fn_scaled_t(t,x0):
        return fn_scaled( [x0[0], x0[1], t, x0[3], x0[4], x0[5], x0[6]])
      npts_guess = int(t_ref_wind*2/(0.5*1e-4))   # Need to have enough points to fully explore the peak, timing to sub-ms accuracy
      print( " Using ", npts_guess, " time points to select the best time, fixing the remaining extrinsic parameters ")
      tvals_scaled_guess = numpy.linspace(0,1,npts_guess)
      lnLvals = numpy.array([-1*fn_scaled_t(t,x0) for t in tvals_scaled_guess])   # note the two -1's cancel
    #    from matplotlib import pyplot as plt;  plt.plot(tvals_scaled_guess,lnLvals,'o'); plt.show(); numpy.savetxt("dump-lnL.dat", numpy.array([tvals_scaled_guess*2*t_ref_wind,lnLvals]).T)
      tbest = tvals_scaled_guess[numpy.argmax(lnLvals)]
      print( " ---- Best extrinsic paramers in MC, after time offset   ---- ")
      P.tref = fiducial_epoch + tbest
      P.print_params()
      lalsimutils.ChooseWaveformParams_array_to_xml([P],"withtime_raw_maxpt_"+opts.output_file) # best point, not including time
      x0[2] = tbest
      x0p = x0
    # Refine best starting time with a search
    #    t_scaled_est = scipy.optimize.brent(fn_scaled_t, brack=(tbest-0.01,tbest,tbest+0.01),args=(x0),maxiter=500)
    #    print "Scaled starting time estimate after, before ",  t_scaled_est, tbest
    #    x0p[2] = t_scaled_est
      t_best_now = lal.LIGOTimeGPS()   # create correct data type
      t_best_now = fiducial_epoch + float((x0p[2]-0.5)*2*t_ref_wind)
      t_best_now_s = int(numpy.floor(t_best_now))
      t_best_now_ns = int((t_best_now - t_best_now_s)*1e9)
      print( " Fitting best (geocentric) time : ", str(t_best_now_s)+ '.'+ str(t_best_now_ns) , " relative time ", (x0p[2]-0.5)*2*t_ref_wind)
      x0  = x0p
    # Full multi-d search for best point
      print( " Starting point [dimensionless] ", x0)
      print( " Starting point [physical] ", [x0[0]*2*numpy.pi, numpy.pi*(x0[1]),  t_ref_wind*2*(x0[2]-0.5), 2*numpy.pi*x0[3], numpy.pi*x0[4], numpy.pi*x0[5], x0[6]*dmax])
      lnLstart = -1*fn_scaled(x0)   # note the two -1's cancel. Compare to 'rho^2/2' value reported by code
      print( " lnL at start :", lnLstart, " [note can be lower than peak because of time offset]: best reported by ILE (including weights) is ", numpy.log(numpy.max(sampler._rvs["integrand"])))
      x = scipy.optimize.fmin(fn_scaled,x0, xtol=1e-5,ftol=1e-3,maxiter=opts.n_max,maxfun=opts.n_max)
      lnLmax = -1.0*fn_scaled(x)
      if lnLmax < lnLstart:
        print( " Maximization failed to improve our initial point ! ")
        x = x0p
        lnLmax = lnLstart
      t_best_now = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind)
      print( "Best lnL = ", lnLmax)
      print( "Best time =", t_best_now)
      print( "Best point", [x[0]*2*numpy.pi, numpy.pi*(x[1]),  t_ref_wind*2*(x[2]-0.5), 2*numpy.pi*x[3], numpy.pi*x[4], numpy.pi*x[5], x[6]*dmax])
      # Output best fit.  Note STANDARD output will occur as usual
      P_max = P.manual_copy()
      P_max.tref = t_best_now
      P_max.phi = float(x[0]*2*numpy.pi)
      P_max.theta = float(numpy.pi*(x[1]))
      P_max.phiref =float(x[3]*numpy.pi*2)
      P_max.incl = float(x[4]*numpy.pi)
      P_max.psi = float(x[5]*numpy.pi)
      P_max.dist = x[6]*dmax*1e6*lal.PC_SI
      P_max.m1 = lal.MSUN_SI*m1
      P_max.m2 = lal.MSUN_SI*m2
      print( " Sanity check: log likelihood for this set is ", factored_likelihood.FactoredLogLikelihood(
                    P_max, rholms, rholms_intp, cross_terms, cross_terms_V, opts.l_max))
      print( " ---- Best extrinsic paramers after polishing   ---- ")
      P_max.print_params()
      P_list = [P_max]
      lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"maxpt_"+opts.output_file)

    # Clear sampler _rvs, to avoid side effects when called again
    for key in sampler._rvs.keys():
      sampler._rvs[key] = []

    return res

def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec):
    nEvals=0
    # PROVENANCE RESET, ON ENTRY, BEFORE ANYTHING CAN FAIL.
    #
    # _rvs_is_pooled describes the record this call is about to build, and it is set by THIS
    # function (the replica block) rather than by a sampler, so no sampler-side per-pass reset
    # can clear it.  Clearing it only on the normal return is not enough: _reject_if_collapsed
    # RAISES after pooling, the caller's `except Exception` swallows that and moves to the next
    # event, and the marker survives.  The next ordinary fair draw is then read as "pooled",
    # _rvs_is_equal_weight goes False, and any consumer that weights rows applies importance
    # weights to rows that already carry them -- the w^2 defect, resurrected on the event after
    # any failure.
    #
    # On ENTRY rather than in a `finally`: entry is reached on every call by construction, and
    # it leaves the state correct even for a caller that never returns normally at all.
    sampler._rvs_is_pooled = False
    P = P_list[indx_event]
    # if pin-distance-to-sim, change the distance prior accordingly
    if opts.pin_distance_to_sim:
      pinned_params['distance']=P.dist/(1.e6 * lal.PC_SI)

    # Call external likleihood preparation if any
    if supplemental_ln_likelhood_prep:
      supplemental_ln_likelhood_prep(P=P,config=supplemental_ln_likelhood_parsed_ini)

    extra_waveform_kwargs = {}
    extra_waveform_kwargs['fd_alignment_postevent_time'] = 2 # 2 seconds after merger for ChooseFDModes
    if opts.internal_waveform_fd_L_frame:
      extra_waveform_kwargs['fd_L_frame'] = True
    if opts.internal_waveform_fd_no_condition:
      extra_waveform_kwargs['no_condition'] = True
    if opts.rom_group:
      extra_waveform_kwargs['rom_taper_start'] = True  # really for NRHyb3dq8 given discussion with Aasim for high-SNR sources, but valuable generally
    if opts.use_gwsignal_lmax_nyquist:
      extra_waveform_kwargs['lmax_nyquist'] = int(opts.use_gwsignal_lmax_nyquist)
    if opts.internal_waveform_extra_lalsuite_args:
      extra_args_dict = eval(opts.internal_waveform_extra_lalsuite_args)  # should only do this once and for all, not in loop!
      print(" Waveform interface: extra args passed ", extra_args_dict)
      extra_waveform_kwargs['extra_waveform_args'] = extra_args_dict
    # Precompute
    t_window = 0.15
    rholms_intp, cross_terms, cross_terms_V,  rholms,  guess_snr, rest=factored_likelihood.PrecomputeLikelihoodTerms(
            fiducial_epoch, t_window, P, data_dict, psd_dict, opts.l_max, fmax,
            False, inv_spec_trunc_Q, T_spec,
            NR_group=NR_template_group,NR_param=NR_template_param,
            use_gwsignal=opts.use_gwsignal,
            use_gwsignal_approx=opts.approximant,
            use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs)

    # skip nan ! Something horrible has happened
    if np.isnan(guess_snr):
      print("  --- NAN SNR GUESS, ABORTING THIS EVENT {} --".format(indx_event))
      raise Exception(" NAN SNR ")

    if opts.auto_logarithm_offset and guess_snr:
      # important: this only impacts *this* analysis
      # important: if we have a very loud signal, it is important to change this dynamically to avoid *underflow*, which sometimes happens otherwise if we fix the scale early on.
      print("    : naive overflow protection: updating lnL overflow based on SNR guess of {} ".format(guess_snr))
      print("    : reminder, diagnostics for lnLmax below are offset by this amount! ")
      # Note by changing the offset, we change how adapt-weight-exponent is being applied (because we've offset the range before squashing it!).  Use with care: pipeline auto-sets the weight exponent
      manual_avoid_overflow_logarithm = guess_snr**2/2 - 100  # more conservative than helper, if I were to auto-set it from --hint-snr.  So the integral peak is more likely positive
      # Don't scale so much that peak guess likelihood is negative:  we do fine for low-ampltiude sources
      if manual_avoid_overflow_logarithm < 0:
        manual_avoid_overflow_logarithm = 0
    elif opts.auto_logarithm_offset:
      print(" PROBLEM: guess_snr not being returned, but auto-tuning requested ")
      # reset to default.  Should not be needed, but weird python scoping error
      manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default 
    else:
      # reset to default.  Should not be needed, but weird python scoping error
      manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default 

    if opts.vectorized:
        lookupNKDict = {}
        lookupKNDict={}
        lookupKNconjDict={}
        ctUArrayDict = {}
        ctVArrayDict={}
        rholmArrayDict={}
        rholms_intpArrayDict={}
        epochDict={}
        for det in rholms_intp.keys():
            print( " Packing ", det)
            lookupNKDict[det],lookupKNDict[det], lookupKNconjDict[det], ctUArrayDict[det], ctVArrayDict[det], rholmArrayDict[det], rholms_intpArrayDict[det], epochDict[det] = factored_likelihood.PackLikelihoodDataStructuresAsArrays( rholms[det].keys(), rholms_intp[det], rholms[det], cross_terms[det],cross_terms_V[det])
            if opts.gpu and (not xpy_default is np):
                lookupNKDict[det] = cupy.asarray(lookupNKDict[det])
                rholmArrayDict[det] = cupy.asarray(rholmArrayDict[det])
                ctUArrayDict[det] = cupy.asarray(ctUArrayDict[det])
                ctVArrayDict[det] = cupy.asarray(ctVArrayDict[det])
                epochDict[det] = cupy.asarray(epochDict[det])



    # Likelihood
    if not opts.time_marginalization:

        def likelihood_function(right_ascension, declination, t_ref, phi_orb,
                inclination, psi, distance):

            dec = numpy.copy(declination).astype(numpy.float64)
            if opts.declination_cosine_sampler:
                dec = numpy.pi/2 - numpy.arccos(dec)
            incl = numpy.copy(inclination).astype(numpy.float64)
            if opts.inclination_cosine_sampler:
                incl = numpy.arccos(incl)

            # use EXTREMELY many bits
            lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128)
            i = 0
            for ph, th, tr, phr, ic, ps, di in zip(right_ascension, dec,
                    t_ref, phi_orb, incl, psi, distance):
                P.phi = ph # right ascension
                P.theta = th # declination
                P.tref = fiducial_epoch + tr # ref. time (rel to epoch for data taking)
                P.phiref = phr # ref. orbital phase
                P.incl = ic # inclination
                P.psi = ps # polarization angle
                P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance

                lnL[i] = factored_likelihood.FactoredLogLikelihood(
                        P, rholms_intp, cross_terms, cross_terms_V, opts.l_max)
                i+=1
            if return_lnL:
              return lnL - manual_avoid_overflow_logarithm 
            return numpy.exp(lnL - manual_avoid_overflow_logarithm)

    else: # Sum over time for every point in other extrinsic params
     if not (opts.rom_integrate_intrinsic or opts.vectorized):
        def likelihood_function(right_ascension, declination, phi_orb, inclination,
                psi, distance):
            dec = numpy.copy(declination).astype(numpy.float64)  # get rid of 'object', and allocate space
            if opts.declination_cosine_sampler:
                dec = numpy.pi/2 - numpy.arccos(dec)
            incl = numpy.copy(inclination).astype(numpy.float64)
            if opts.inclination_cosine_sampler:
                incl = numpy.arccos(incl)

            # use EXTREMELY many bits
            lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128)
            i = 0
            tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally

            for ph, th, phr, ic, ps, di in zip(right_ascension, dec,
                    phi_orb, inclination, psi, distance):
                P.phi = ph # right ascension
                P.theta = th # declination
                P.tref = fiducial_epoch  # see 'tvals', above
                P.phiref = phr # ref. orbital phase
                P.incl = ic # inclination
                P.psi = ps # polarization angle
                P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance


                lnL[i] = factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals,
                        P, rholms_intp, rholms, cross_terms, cross_terms_V,                   
                        opts.l_max,interpolate=opts.interpolate_time)
                i+=1
            if supplemental_ln_likelihood:
              lnL += supplemental_ln_likelihood(right_ascension, declination, phi_orb,inclination, psi, distance)
            if return_lnL:
              return lnL -manual_avoid_overflow_logarithm 
            return numpy.exp(lnL -manual_avoid_overflow_logarithm)

     elif opts.vectorized: # use array-based multiplications, fewer for loops
        if (not opts.gpu):
          if opts.distance_marginalization:
            print(" **Warning**: distance marginalization not being used, this is the old vectorized code path not the xpy path.  Check your code path; you may want --force-xpy ")
          def likelihood_function(right_ascension, declination, phi_orb, inclination,
                psi, distance):
#            global nEvals
            tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally
            dec = numpy.copy(declination).astype(numpy.float64)
            if opts.declination_cosine_sampler:
              dec = numpy.pi/2 - numpy.arccos(dec)
            incl = numpy.copy(inclination).astype(numpy.float64)
            if opts.inclination_cosine_sampler:
              incl = numpy.arccos(incl)

            # use EXTREMELY many bits
            lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128)
            P.phi = right_ascension.astype(float)  # cast to float
            P.theta = dec #declination.astype(float)
            P.tref = float(fiducial_epoch)
            P.phiref = phi_orb.astype(float)
            P.incl = incl #inclination.astype(float)
            P.psi = psi.astype(float)
            P.dist = (distance* 1.e6 * lalsimutils.lsu_PC).astype(float) # luminosity distance

            # rotate sky if needed
            if opts.internal_sky_network_coordinates:
                  P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi)
                  P.theta = np.pi/2 - P.theta
                  P.phi = xpy_default.mod(P.phi, 2*np.pi)

            # rotate phase if needed
            # make copies of arrays
            #   Sampling assumes P.phiref == phi+ \psi ,   P.psi == phi - psi
            if opts.internal_rotate_phase:
              phi_orb_true = (P.phiref + P.psi)/2.
              psi_true = (P.phiref - P.psi)/2.
              P.psi= psi_true
              P.phiref = phi_orb_true

            lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals,
                        P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max)
#            nEvals +=len(right_ascension)
            if supplemental_ln_likelihood:
              lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, P.dist) # use these variables so they are already float-type
            if return_lnL:
              return lnL -manual_avoid_overflow_logarithm
            return numpy.exp(lnL-manual_avoid_overflow_logarithm)
        else: # vectorized and gpu either available or being forced to use xpy code path
            print( " Using CUDA GPU likelihood, if cupy available ")
            if not opts.distance_marginalization:
              def likelihood_function(right_ascension, declination, phi_orb, inclination,
                psi, distance):
#                global nEvals
                tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally
                P.phi = xpy_asarray_already(right_ascension)  # cast to float
                if opts.declination_cosine_sampler:
                  P.theta = numpy.pi/2 - xpy_default.arccos(xpy_asarray_already(declination))
                else:
                  P.theta = xpy_asarray_already(declination)
                P.tref = float(fiducial_epoch)
                P.phiref = xpy_asarray_already(phi_orb)
                if opts.inclination_cosine_sampler:
                  P.incl = xpy_default.arccos(xpy_asarray_already(inclination))
                else:
                  P.incl = xpy_asarray_already(inclination)
                P.psi = xpy_asarray_already(psi)
                P.dist = xpy_asarray_already(distance* 1.e6 * lalsimutils.lsu_PC) # luminosity distance

                # rotate sky if needed
                if opts.internal_sky_network_coordinates:
                  P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi)
                  P.theta = np.pi/2 - P.theta
                  P.phi = xpy_default.mod(P.phi, 2*np.pi)

                # rotate phase if needed
                # make copies of arrays
                #   Sampling assumes P.phiref == phi+ \psi ,   P.psi == phi - psi
                if opts.internal_rotate_phase:
                  phi_orb_true = (P.phiref + P.psi)/2.
                  psi_true = (P.phiref - P.psi)/2.
                  P.psi= psi_true
                  P.phiref = phi_orb_true


                lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,
                        P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default)
#                nEvals +=len(right_ascension)
                if supplemental_ln_likelihood:
                  lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, P.dist,xpy=xpy_default) # use these variables so they are already float-type
                if return_lnL:
                  return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm)
                return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm))
            else:
              print( " Using direct distance marginalization  ")
              xmin = factored_likelihood.distMpcRef / dmax
              xmax = factored_likelihood.distMpcRef / dmin
              bmax = xpy_default.asarray(lookup_table["bmax"])
              sqrt_bmax = xpy_default.sqrt(bmax)
              bref = xpy_default.asarray(lookup_table["bref"])
              s_array = xpy_default.asarray(lookup_table["s_array"])
              smin = s_array[0]
              smax = s_array[-1]
              t_array = xpy_default.asarray(lookup_table["t_array"])
              tmax = t_array[-1]
              lnI_array = xpy_default.asarray(lookup_table["lnI_array"])

              intp = EvenBivariateLinearInterpolator(s_array[0], s_array[1] - s_array[0], t_array[0], t_array[1] - t_array[0], lnI_array)

              def exponent_max(x0, b):
                x0_expmax = xpy_default.clip(x0, a_min=xmin, a_max=xmax)
                return b * x0_expmax * (x0 - 0.5*x0_expmax)

              def b_to_t(b):
                # TODO: this function is duplicate with what is in util_InitMargTable
#                return np.arcsinh(b / bref)
                b_by_bref = b / bref
                return xpy_default.arcsinh(b_by_bref, out=b_by_bref)


              def x0_to_s(x0):
                # TODO: this function is duplicate with what is in util_InitMargTable
#                return np.arcsinh(np.sqrt(bmax) * (x0 - xmin)) - np.arcsinh(np.sqrt(bmax) * (xmax - x0))
                A = x0 - xmin
                A *= sqrt_bmax
                xpy_default.arcsinh(A, out=A)

                B = xmax - x0
                B *= sqrt_bmax
                xpy_default.arcsinh(B, out=B)

                A -= B
                return A

              def distmarg_loglikelihood(kappa_sq, rho_sq):
                x0 = kappa_sq / rho_sq
#                lnI = np.ones(shape=x0.shape) * -np.inf
                lnI = xpy_default.full_like(x0, xpy_default.NINF)
                s = x0_to_s(x0)
                t = b_to_t(rho_sq)
#                in_bounds = (s > smin) * (s < smax) * (t < tmax)
                in_bounds = (s > smin) & (s < smax) & (t < tmax)
                lnI[in_bounds] = intp(s[in_bounds], t[in_bounds])
                return exponent_max(x0, rho_sq) + lnI

              if lookup_table["phase_marginalization"]:

                print( " Using direct phase marginalization  ")
                for det in lookupNKDict:
                  if set((lm[0], lm[1]) for lm in lookupNKDict[det]) != {(2, 2), (2, -2)}:
                    raise Exception(
                        " Phase marginalization is implemented only for 2-2 modes, "
                        f"while the modes consired here are {lookupNKDict[det]}."
                    )

                def likelihood_function(right_ascension, declination, inclination, psi):
#                  global nEvals
                  tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally
                  P.phi = xpy_default.asarray(right_ascension, dtype=np.float64)  # cast to float
                  if opts.declination_cosine_sampler:
                    P.theta = numpy.pi/2 - xpy_default.arccos(xpy_default.asarray(declination,dtype=np.float64))
                  else:
                    P.theta = xpy_default.asarray(declination,dtype=np.float64)
                  P.tref = float(fiducial_epoch)
                  P.phiref = xpy_default.full_like(inclination, 0., dtype=np.float64)
                  if opts.inclination_cosine_sampler:
                    P.incl = xpy_default.arccos(xpy_default.asarray(inclination, dtype=np.float64))
                  else:
                    P.incl = xpy_default.asarray(inclination, dtype=np.float64)
                  P.psi = xpy_default.asarray(psi, dtype=np.float64)
                  P.dist = xpy_default.asarray(factored_likelihood.distMpcRef * 1.e6 * lalsimutils.lsu_PC,dtype=np.float64) # luminosity distance

                  # rotate sky if needed
                  if opts.internal_sky_network_coordinates:
                    P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi)
                    P.theta = np.pi/2 - P.theta
                    P.phi = xpy_default.mod(P.phi, 2*np.pi)

                  # rotate phase if needed
                  # make copies of arrays
                  #   Sampling assumes P.phiref == phi+ \psi ,   P.psi == phi - psi
                  if opts.internal_rotate_phase:
                    phi_orb_true = (P.phiref + P.psi)/2.
                    psi_true = (P.phiref - P.psi)/2.
                    P.psi= psi_true
                    P.phiref = phi_orb_true

                  lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,
                    P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood, phase_marginalization=True)
#                  nEvals +=len(right_ascension)
                  if supplemental_ln_likelihood:
                    lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, 0,xpy=xpy_default) # Same API
                  if return_lnL:
                    return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm)
                  return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm))

              else:

                def likelihood_function(right_ascension, declination, phi_orb, inclination, psi):
#                  global nEvals
                  tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally
                  P.phi = xpy_default.asarray(right_ascension, dtype=np.float64)  # cast to float
                  if opts.declination_cosine_sampler:
                    P.theta = numpy.pi/2 - xpy_default.arccos(xpy_default.asarray(declination,dtype=np.float64))
                  else:
                    P.theta = xpy_default.asarray(declination,dtype=np.float64)
                  P.tref = float(fiducial_epoch)
                  P.phiref = xpy_default.asarray(phi_orb, dtype=np.float64)
                  if opts.inclination_cosine_sampler:
                    P.incl = xpy_default.arccos(xpy_default.asarray(inclination, dtype=np.float64))
                  else:
                    P.incl = xpy_default.asarray(inclination, dtype=np.float64)
                  P.psi = xpy_default.asarray(psi, dtype=np.float64)
                  P.dist = xpy_default.asarray(factored_likelihood.distMpcRef * 1.e6 * lalsimutils.lsu_PC,dtype=np.float64) # luminosity distance

                  # rotate sky if needed
                  if opts.internal_sky_network_coordinates:
                    P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi)
                    P.theta = np.pi/2 - P.theta
                    P.phi = xpy_default.mod(P.phi, 2*np.pi)

                  # rotate phase if needed
                  # make copies of arrays
                  #   Sampling assumes P.phiref == phi+ \psi ,   P.psi == phi - psi
                  if opts.internal_rotate_phase:
                    phi_orb_true = (P.phiref + P.psi)/2.
                    psi_true = (P.phiref - P.psi)/2.
                    P.psi= psi_true
                    P.phiref = phi_orb_true

                  lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,
                    P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood)
#                  nEvals +=len(right_ascension)
                  if supplemental_ln_likelihood:
                    lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, 0,xpy=xpy_default) # Same API
                  if return_lnL:
                    return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm)
                  return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm))

     else: # integrate over intrinsic variables. Right now those variables ahave HARDCODED NAMES, alas
        def likelihood_function(right_ascension, declination, phi_orb, inclination,
                psi, distance,q):
            dec = numpy.copy(declination).astype(numpy.float64)
            if opts.declination_cosine_sampler:
                dec = numpy.pi/2 - numpy.arccos(dec)
            incl = numpy.copy(inclination).astype(numpy.float64)
            if opts.inclination_cosine_sampler:
                incl = numpy.arccos(incl)

            lnL = numpy.zeros(len(right_ascension),dtype=numpy.float128)
#            i = 0
            tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT))  # choose an array at the target sampling rate. P is inherited globally


#            t_start =lal.GPSTimeNow() 
            for ph, th, phr, ic, ps, di,qi in zip(right_ascension, dec,
                    phi_orb, incl, psi, distance,q):
                 # Reconstruct U,V using ROM fits.  PROBABLY should do this once for every q, rather than deep on the loop
                P.assign_param('q',qi)  # mass ratio
                rholms_intp_A, cross_terms_A, cross_terms_V_A, rholms_A, rest_A = factored_likelihood.ReconstructPrecomputedLikelihoodTermsROM(P, rest, rholms_intp, cross_terms, cross_terms_V, rholms,verbose=False)
                # proceed for rest
                P.phi = ph # right ascension
                P.theta = th # declination
                P.tref = fiducial_epoch  # see 'tvals', above
                P.phiref = phr # ref. orbital phase
                P.incl = ic # inclination
                P.psi = ps # polarization angle
                P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance


                lnL[i] = factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals,
                        P, rholms_intp_A, rholms_A, cross_terms_A, cross_terms_V_A,                   
                        opts.l_max,interpolate=opts.interpolate_time)
                if numpy.isnan(lnL[i]) or lnL[i]<-200:
                    lnL[i] = -200   # regularize  : a hack, for now, to deal with rare ROM problems. Only on the ROM logic fork
                i+=1
#            t_end =lal.GPSTimeNow() 
#            print " Cost per evaluation ", (t_end - t_start)/len(q)
#            print " Max lnL for this iteration ", numpy.max(lnL)
            if return_lnL:
              return lnL -manual_avoid_overflow_logarithm
            return numpy.exp(lnL - manual_avoid_overflow_logarithm)

    if opts.sampler_method == "adaptive_cartesian_gpu":
      # reset sampling parameter as needed (distance, inclination but not sky location)
      # distance (and inclination) can conceivably correlate with mass, and we don't want to truncate in distance/inclination prematurely from history effects
      # method ONLY implemented for adaptcart!
      if 'distance' in sampler.params:
        sampler.reset_sampling('distance')
      sampler.reset_sampling('inclination')
    elif use_gmm_args:   # standalone GMM or a portfolio with a GMM member (gmm_dict exists in both)
      if 'distance' in sampler.params:
         pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination'])
         if  pair_d_incl in gmm_dict:
            gmm_dict[pair_d_incl] = None # reset inclination/distance sampling
      single_incl = sampler_param_tuple(sampler, ['inclination'])
      if (opts.distance_marginalization) and single_incl in gmm_dict:
        gmm_dict[single_incl] = None # reset inclination sampling


    # Integrate
    args = likelihood_function.__code__.co_varnames[:likelihood_function.__code__.co_argcount]
    print( " --------> Arguments ", args)
    # if exactly zero likelihood function
    like_to_integrate = likelihood_function
    if opts.zero_likelihood:
              if opts.internal_use_lnL:
                   like_to_integrate = zero_like
              else:
                   like_to_integrate = unit_like
    if oracleRS:
         print("  ORACLE  - seeding sampling  with" , oracleRS.params_ordered)
         if hasattr(sampler, 'update_sampling_prior'):
             rvs_train = {}
             _, _, rv_oracle = oracleRS.draw_simplified(opts.n_chunk)
             for indx,p  in enumerate(sampler.params_ordered):
               rvs_train[p] = rv_oracle[:,indx]
             # train with equal weight - no likelihood information
             lnL_oracles  = np.zeros(opts.n_chunk)
             sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level)

    _maybe_load_av_state(sampler)
    _maybe_enable_anisotropic_bins(sampler)
    res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params)

    # L0 auto-rescue: on a very sharply-peaked (high-amplitude) point a cold AV can stall at
    # n_eff ~ 1 because it never draws near the tiny peak.  If so, seed a SECOND pass from this
    # same point's own highest-likelihood samples and re-run.  Opt-in via
    # --sampler-warmstart-retry-neff.  MUST run BEFORE the not(res) guard below: a degenerate
    # early termination returns (None,None,None,None) and is the strongest rescue trigger, so
    # raising on it first would skip exactly the case the rescue exists for.
    res, var, neff, dict_return = _maybe_l0_rescue(
        sampler, res, var, neff, dict_return,
        like_to_integrate, unpinned_params, pinned_params,
        lnL_offset=manual_avoid_overflow_logarithm)

    if not(res): # no resut
      raise ValueError(" No integral result returned")

    if not(opts.internal_use_lnL):
      log_res = numpy.log(res)
      sqrt_var_over_res =  numpy.sqrt(var)/res
    else:
      log_res = res
      sqrt_var_over_res =  numpy.exp(var/2 - log_res)

    # MC-error diagnostics, the collapse report/gate, and replica replication+pooling.  This
    # OWNS both collapse-gate calls (first run and pooled verdict) -- do not add a separate
    # _report_and_gate_collapse call here, or the pooled one gets bypassed.  It needs
    # log_res/sqrt_var_over_res, hence its position after they are computed; the main driver
    # has the same ordering inline.
    res, var, neff, log_res, sqrt_var_over_res, dict_return = _maybe_replicate_for_mc_error(
        sampler, res, var, neff, dict_return, log_res, sqrt_var_over_res,
        like_to_integrate, unpinned_params, pinned_params,
        lnL_offset=manual_avoid_overflow_logarithm)

    # Report results
    if opts.output_file:
        fname_output_txt = opts.output_file +"_"+str(indx_event)+"_" + ".dat"
        m1 =P.m1/lal.MSUN_SI
        m2 =P.m2/lal.MSUN_SI
        if opts.sim_xml: 
            event_id = opts.event
        else:
            event_id = -1
        if opts.event == None:
            event_id = -1
        if opts.save_eccentricity:
            # output format when eccentricity is being used
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.eccentricity, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
        elif not (P.lambda1>0 or P.lambda2>0):
          # output format when lambda is NOT used
          if not opts.pin_distance_to_sim:
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
          else:
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, pinned_params["distance"], log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
        else:
          if not(opts.export_eos_index):
            # Alternative output format if lambda is active
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.lambda1, P.lambda2, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"
          else:
            numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z,  P.lambda1, P.lambda2, P.eos_table_index, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]]))  #dict_return["convergence_test_results"]["normal_integral]"


    # Comprehensive output (not yet provided)
    # Convert declination, inclination  parameters in sampler if needed
    if opts.save_samples and opts.output_file:
      import copy
      samples = copy.deepcopy(sampler._rvs)  # deep copy: avoid modifying structures and  having side effect on integrator, which loops over keys Expensive!
      # A POOLED record (--mc-error-replicas) is weighted BETWEEN blocks by the replica
      # evidences, and nothing below preserves those weights: convert it to an equal-weight
      # draw BEFORE anything consumes it -- including the time resampler, which picks a time per
      # row and so assumes the rows already are the posterior.  A no-op otherwise.
      samples = _export_rvs_equal_weight(samples, sampler, use_lnL=rvs_integrand_is_lnL)
      # Insert reference distance if it was marginalized over
      if "distance" not in samples:
        # Not distance output is the same as internal calculations: in *Mpc*
        samples["distance"] = np.full_like(
          samples["psi"],
          factored_likelihood.distMpcRef,  #*1e6*lal.PC_SI,
          )
      # Insert reference phase if it was marginalized over
      if "phi_orb" not in samples:
        samples["phi_orb"] = np.full_like(samples["psi"], 0.)
      if opts.inclination_cosine_sampler:
        samples["inclination"] = numpy.arccos(samples["inclination"].astype(numpy.float64))
      if opts.declination_cosine_sampler:
        samples["declination"] = numpy.pi/2 - numpy.arccos(samples["declination"].astype(numpy.float64))
      xmldoc = ligolw.Document()
      xmldoc.appendChild(ligolw.LIGO_LW())
      process.register_to_xmldoc(xmldoc, sys.argv[0], opts.__dict__)
      if not(opts.resample_time_marginalization): 
        if not opts.time_marginalization:
            samples["t_ref"] += float(fiducial_epoch)
        else:
            samples["t_ref"] = float(fiducial_epoch)*numpy.ones(len(samples["psi"]))
      # rotate sky to recover physical coordinates.  Note on CPU
      if opts.internal_sky_network_coordinates:
        tmp_th, tmp_ph = my_rotation_cpu(np.pi/2 - samples['declination'],samples['right_ascension'])
        samples['declination'] = np.pi/2 - tmp_th
        samples['right_ascension'] = np.mod(tmp_ph, 2*np.pi)

      # recover true phase, polarization samples
      if opts.internal_rotate_phase:
        phi_orb_true = np.mod((samples['phi_orb'] + samples['psi'])/2., 2*np.pi)
        psi_true = np.mod((samples['phi_orb'] - samples['psi'])/2., 2*np.pi)   # keep as 2 pi range, to be consistent with past work
        samples['psi']= psi_true
        samples['phi_orb'] = phi_orb_true
      samples["polarization"] = samples["psi"]
      samples["coa_phase"] = samples["phi_orb"]        
      if ("declination", "right_ascension") in sampler.params:
            samples["latitude"], samples["longitude"] = samples[("declination", "right_ascension")]
      else:
            samples["latitude"] = samples["declination"]
            samples["longitude"] = samples["right_ascension"]
      if "log_integrand" in samples:
        samples["loglikelihood"] = samples["log_integrand"] +  manual_avoid_overflow_logarithm
      else:
        samples["loglikelihood"] = numpy.log(samples["integrand"]) + manual_avoid_overflow_logarithm  # export with consistent offset
      if not opts.rom_integrate_intrinsic:
            # ILE mode: insert fixed model parameters
            samples["mass1"] = numpy.ones(samples["psi"].shape)*m1 # opts.mass1
            samples["mass2"] = numpy.ones(samples["psi"].shape)*m2 # opts.mass2
            samples["spin1x"] =numpy.ones(samples["psi"].shape)*P.s1x
            samples["spin1y"] =numpy.ones(samples["psi"].shape)*P.s1y
            samples["spin1z"] =numpy.ones(samples["psi"].shape)*P.s1z
            samples["spin2x"] =numpy.ones(samples["psi"].shape)*P.s2x
            samples["spin2y"] =numpy.ones(samples["psi"].shape)*P.s2y
            samples["spin2z"] =numpy.ones(samples["psi"].shape)*P.s2z
            samples["alpha4"] =numpy.ones(samples["psi"].shape)*P.eccentricity
            samples["alpha5"] =numpy.ones(samples["psi"].shape)*P.lambda1
            samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2
            # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so
            samples["alpha2"] = numpy.zeros(samples["psi"].shape)
            samples["alpha3"] = numpy.zeros(samples["psi"].shape)

      if opts.resample_time_marginalization:
        samples = resample_samples(samples,lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict)
        samples["loglikelihood" ] = samples["lnL_raw"]  # export the non-time-marginalized likelihood, if we are in the final stages
#        print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref']))
      xmlutils.append_samples_to_xmldoc(xmldoc, samples)
        # Extra metadata
      dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal}
#      if 'distance' in pinned_params:
#            dict_out['distance'] = pinned_params["distance"]
      converged_result = False
      if "convergence_test_results" in dict_return:
        converged_result = dict_return["convergence_test_results"]["normal_integral"]
      xmlutils.append_likelihood_result_to_xmldoc(xmldoc, log_res+manual_avoid_overflow_logarithm, neff=neff, converged=converged_result, **dict_out)
      fname_output_xml = opts.output_file +"_"+str(indx_event)+"_" + ".xml.gz"
      utils.write_filename(xmldoc, fname_output_xml, compress="gz")




    if opts.maximize_only and opts.output_file:
      # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization)
      if "integrand" in sampler._rvs:
        indx_guess = numpy.argmax(sampler._rvs["integrand"])   # start search near maximum-likelihood point. (WARNING: can be very close by)
      else:
        indx_guess = numpy.argmax(sampler._rvs["log_integrand"])
      P.radec=True
      P.phi = sampler._rvs["right_ascension"][indx_guess]
      P.theta = sampler._rvs["declination"][indx_guess]
      P.phiref = sampler._rvs["phi_orb"][indx_guess]
      P.incl = sampler._rvs["inclination"][indx_guess]
      P.psi = sampler._rvs["psi"][indx_guess]
      P.dist = sampler._rvs["distance"][indx_guess]*1e6*lal.PC_SI
      print( " ---- Best extrinsic paramers in MC   ---- ")
      P.print_params()
      lalsimutils.ChooseWaveformParams_array_to_xml([P],"notime_raw_maxpt_"+opts.output_file) # best point, not including time

      import scipy.optimize
      def fn_scaled(x):
          P.phi = float(x[0]*2*numpy.pi) # right ascension
          P.theta = float((x[1])*numpy.pi) # declination, really polar angle. NOT zero on equator
          P.tref = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) # ref. time (rel to epoch for data taking)
          P.phiref = float(x[3]*numpy.pi*2) # ref. orbital phase
          P.incl = float(x[4]*numpy.pi) # inclination
          P.psi = float(x[5]*numpy.pi) # polarization angle
          P.dist = x[6]*dmax* 1.e6 * lalsimutils.lsu_PC # luminosity distance
    
          return -1.0* factored_likelihood.FactoredLogLikelihood(
                    P, rholms,rholms_intp, cross_terms, cross_terms_V, opts.l_max)

    # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization)
      x0 =numpy.array( [ \
       sampler._rvs["right_ascension"][indx_guess]/(2*numpy.pi) , \
       (sampler._rvs["declination"][indx_guess]/numpy.pi),  \
       0.5, \
       (sampler._rvs["phi_orb"][indx_guess]/(2*numpy.pi)), \
       sampler._rvs["inclination"][indx_guess]/(numpy.pi),\
       sampler._rvs["psi"][indx_guess]/numpy.pi,\
       sampler._rvs["distance"][indx_guess]/dmax\
            ],dtype=numpy.float128)
      x0 = numpy.fmod(x0,numpy.ones(len(x0)))  # had BETTER be defined on this range!
      # Pick the best starting time. BRUTE FORCE METHOD: use grid
      def fn_scaled_t(t,x0):
        return fn_scaled( [x0[0], x0[1], t, x0[3], x0[4], x0[5], x0[6]])
      npts_guess = int(t_ref_wind*2/(0.5*1e-4))   # Need to have enough points to fully explore the peak, timing to sub-ms accuracy
      print( " Using ", npts_guess, " time points to select the best time, fixing the remaining extrinsic parameters ")
      tvals_scaled_guess = numpy.linspace(0,1,npts_guess)
      lnLvals = numpy.array([-1*fn_scaled_t(t,x0) for t in tvals_scaled_guess])   # note the two -1's cancel
    #    from matplotlib import pyplot as plt;  plt.plot(tvals_scaled_guess,lnLvals,'o'); plt.show(); numpy.savetxt("dump-lnL.dat", numpy.array([tvals_scaled_guess*2*t_ref_wind,lnLvals]).T)
      tbest = tvals_scaled_guess[numpy.argmax(lnLvals)]
      print( " ---- Best extrinsic paramers in MC, after time offset   ---- ")
      P.tref = fiducial_epoch + tbest
      P.print_params()
      lalsimutils.ChooseWaveformParams_array_to_xml([P],"withtime_raw_maxpt_"+opts.output_file) # best point, not including time
      x0[2] = tbest
      x0p = x0
    # Refine best starting time with a search
    #    t_scaled_est = scipy.optimize.brent(fn_scaled_t, brack=(tbest-0.01,tbest,tbest+0.01),args=(x0),maxiter=500)
    #    print "Scaled starting time estimate after, before ",  t_scaled_est, tbest
    #    x0p[2] = t_scaled_est
      t_best_now = lal.LIGOTimeGPS()   # create correct data type
      t_best_now = fiducial_epoch + float((x0p[2]-0.5)*2*t_ref_wind)
      t_best_now_s = int(numpy.floor(t_best_now))
      t_best_now_ns = int((t_best_now - t_best_now_s)*1e9)
      print( " Fitting best (geocentric) time : ", str(t_best_now_s)+ '.'+ str(t_best_now_ns) , " relative time ", (x0p[2]-0.5)*2*t_ref_wind)
      x0  = x0p
    # Full multi-d search for best point
      print( " Starting point [dimensionless] ", x0)
      print( " Starting point [physical] ", [x0[0]*2*numpy.pi, numpy.pi*(x0[1]),  t_ref_wind*2*(x0[2]-0.5), 2*numpy.pi*x0[3], numpy.pi*x0[4], numpy.pi*x0[5], x0[6]*dmax])
      lnLstart = -1*fn_scaled(x0)   # note the two -1's cancel. Compare to 'rho^2/2' value reported by code
      print( " lnL at start :", lnLstart, " [note can be lower than peak because of time offset]: best reported by ILE (including weights) is ", numpy.log(numpy.max(sampler._rvs["integrand"])))
      x = scipy.optimize.fmin(fn_scaled,x0, xtol=1e-5,ftol=1e-3,maxiter=opts.n_max,maxfun=opts.n_max)
      lnLmax = -1.0*fn_scaled(x)
      if lnLmax < lnLstart:
        print( " Maximization failed to improve our initial point ! ")
        x = x0p
        lnLmax = lnLstart
      t_best_now = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind)
      print( "Best lnL = ", lnLmax)
      print( "Best time =", t_best_now)
      print( "Best point", [x[0]*2*numpy.pi, numpy.pi*(x[1]),  t_ref_wind*2*(x[2]-0.5), 2*numpy.pi*x[3], numpy.pi*x[4], numpy.pi*x[5], x[6]*dmax])
      # Output best fit.  Note STANDARD output will occur as usual
      P_max = P.manual_copy()
      P_max.tref = t_best_now
      P_max.phi = float(x[0]*2*numpy.pi)
      P_max.theta = float(numpy.pi*(x[1]))
      P_max.phiref =float(x[3]*numpy.pi*2)
      P_max.incl = float(x[4]*numpy.pi)
      P_max.psi = float(x[5]*numpy.pi)
      P_max.dist = x[6]*dmax*1e6*lal.PC_SI
      P_max.m1 = lal.MSUN_SI*m1
      P_max.m2 = lal.MSUN_SI*m2
      print( " Sanity check: log likelihood for this set is ", factored_likelihood.FactoredLogLikelihood(
                    P_max, rholms, rholms_intp, cross_terms, cross_terms_V, opts.l_max))
      print( " ---- Best extrinsic paramers after polishing   ---- ")
      P_max.print_params()
      P_list = [P_max]
      lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"maxpt_"+opts.output_file)

    # Clear sampler _rvs, to avoid side effects when called again
    for key in sampler._rvs.keys():
      sampler._rvs[key] = []

    return res


lnL_sofar = -np.inf
no_adapt_sky = False
for indx in numpy.arange(len(P_list)):
 try:
  if opts.LISA:
    res = analyze_event_LISA(P_list, indx, data_dict, psd_dict, fmax, opts)
  else:
    res = analyze_event(P_list, indx, data_dict, psd_dict, fmax, opts)
  # abort if horrible (nan event) - done with 'raise'
  lnL_sofar = np.max([lnL_sofar,res])
  if opts.force_reset_all:  # depends on integrator!  May not always be availble
    if opts.sampler_method == "adaptive_cartesian_gpu": 
      for name in sampler.params:
        sampler.reset_sampling(name)
    elif use_gmm_args:   # standalone GMM or a portfolio with a GMM member
      # reset the GMM dictionary
      for component in gmm_dict:
          gmm_dict[component] = None
    elif opts.sampler_method == 'AV':
      print(" AV always resets every iteration ! ")
    else:
      print(" force-reset-all not defined for this integrator ")
      sys.exit(1)
  if opts.no_adapt_after_first and (not no_adapt_sky):
    if lnL_sofar > 20:  # Use absolute threshold.  Expect this will give modest sky localization.
      # remove right_ascension, declination from adaptive parameters
      params_adapt = sampler.adaptive
      params_adapt  = list(set(params_adapt) - set(['right_ascension','declination']))
      sampler.adaptive = params_adapt
      # Disable saving of the integrand -- saves on memory and will reduce speed. 
      # Note this needs to be done in a few places
      pinned_params.update({"force_no_adapt":True,"save_intg":False, "igrand_threshold_deltalnL":20})  # massively reduce memory usage in logic branch, don't save all. Highly redundant sequence to self-document all related logic branches
 except Exception as exception_failure:
  print("  ===> FAILED ANALYSIS <==== ")
  print(exception_failure)
  if opts.internal_make_empty_file_on_error:
    fname_output_txt = opts.output_file +"_"+str(indx)+"_" + ".dat"
    open(fname_output_txt,'a').close()  # create empty file
  if opts.internal_hard_fail_on_error:
    sys.exit(1)
#  if len(exception_failure) >0:
#    if "CUBLAS" in exception_failure[0]:  # Hard fail if a cuda error!
#      sys.exit(1) 
  if ("CUDA" in str(exception_failure)) or ('CUBLAS' in str(exception_failure)) or ('cuda' in str(exception_failure)) or ('compilation' in (str(exception_failure)) ): # Hard fail if a cuda error with a catchable error code
    if (opts.internal_soft_fail_on_cuda_error):
      sys.exit(0)
    sys.exit(62)
  if 'Out of memory' in str(exception_failure):   # should never happen, most likely a crappy node or failure to set correct memory limit for ILE.
    sys.exit(63)
  str_err = " {} ".format(exception_failure)
  if ('Zero prior failure' in str_err) or ('effective samples' in str_err):
    # common failures that require reset of sampler
    #    - zero prior : very rare case where the prior is exactly zero for some reason. User error most likely (floors, etc). Should never happen
    #    - effective samples = nan : error with AC where the integrator gets confused. Reset
    # reset all variables sampling, so we don't contaminate subsequent versions
    #    - again, this should only be a problem if we are using multiply adaptive sampling, so it should NEVER happen on the first time through
    for name in sampler.params:
      sampler.reset_sampling(param)

  #a sketch of how I would do custom failure modes, but for sake of speed I'm just setting the above right now
  #for i,failure_mode in enumerate(opts.custom_fails):
  #  if failure_mode in str(exception_failure):
  #    sys.exit(opts.custom_fail_codes[i])

  print( " Probable reasons: SEOB nyquist or starting frequency limit or signal duration ")
  print( " Skipping the following binary! ")
  # Zero out extrinsic parameters -- these are CUDA-populated / meaningless, but could cause errors if populated
  P_list[indx].incl = P_list[indx].tref = P_list[indx].dist = P_list[indx].phiref = P_list[indx].psi =P_list[indx].theta = P_list[indx].phi =0
  P_list[indx].print_params()
  
