#! /usr/bin/env python
#
#  GOAL
#    Simple pipeline: just ILE+fitting iteration jobs
#    User must provide arglists for both ILE and fitting/iteration jobs
#    Assumes user has done all necessary setup (e.g., PSDs, data selection, picking channel names, etc)
#    Will add in fit assessment jobs later
#
#  WHY DO THIS?
#    We're often reanalyzing a single stretch of data again and again. (Or otherwise have some simple setup).
#    This saves time for large-scale injection studies or systematics reruns
#
#    Intend to replace naive iteration code with reproducible, extendable/scalable workflow
#
#  SEE ALSO
#     - original shell-script-based generation scripts
#     - create_event_parameter_pipeline  # long term plan
#     - create_postprocessing_event_dag.py
#     - create_event_dag_via_grid   # basic iteration code
#
# WORKFLOW STRUCTURE
#    - Required elements:
#          cache file  (specified in  --cache <cache> in cip_args)
#          PSD files (.... --psd  in cip_args)
#          channel names
#          event time
#    - Inputs
#         args_cip.txt   # no --fname, --fname-output-...  (overrridden)
#         args_ile.txt    # no --sim-xml,  --event (overrridden)      ... though you may want to add them for command-single.sh to run properly
#    - Workspaces
#          Creates one directory for each task, for each iteration
#          Main directory:
#              - output-ILE-<iteration>.xml.gz, starting with iteration 0 (=input)
#              - *.composite files, from each consolidation step
#              - all.net : union of all *.composite files 
#   
# EXAMPLES
#    python create_event_parameter_pipeline_BasicIteration.py --ile-args args_ile.txt --cip-args args.txt 
#
# TEST SEQUENCE
#        unixhome/Projects/LIGO-ILE-Applications/communications/20180806-Me-PipelineDevelopment
#    * Trivial iteration, fitting a constant likelihood at each step: 
#           rm -rf test_workflow_trivial; make test_workflow_trivial
#         


###
###  PLANS FOR IMPROVEMENT
###
#
# * Stress-testing final fit
#       - single-iteration runs which change lnL and/or coordinates
#       - note this can be done with a DRIVER, does not need to be hardcoded?
# * BCR/BCI (single-ifo vs multi-ifo run)
#       - runs which do PE on single IFOs and 
#       - again, this can be done with a DRIVER

import argparse
import sys
import os
import shutil
import numpy as np
import RIFT.lalsimutils as lalsimutils
import lalsimulation as lalsim
import lal
import functools
import itertools
import re

# Backend-neutral pipeline namespace (htcondor/glue/slurm) provided by dag_utils_generic
from RIFT.misc.dag_utils_generic import pipeline

import RIFT.misc.dag_utils_generic as dag_utils
from RIFT.misc.dag_utils_generic import mkdir
from RIFT.misc.dag_utils_generic import which


def parse_ile_args_for_bw(my_arg_str):
    # Copied verbatim from ILE
    from optparse import OptionParser, OptionGroup
    optp = OptionParser()
    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("-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-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-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=50*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=int,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("--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-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("--verbose",action='store_true')
    integration_params = OptionGroup(optp, "Integration Parameters", "Control the integration with these options.")
    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("--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.")
    integration_params.add_option("--no-adapt", action="store_true", help="Turn off adaptive sampling. Adaptive sampling is on by default.")
    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) 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("--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")
    optp.add_option_group(integration_params)

    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("--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")
    optp.add_option_group(intrinsic_params)

    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)


    opts, args = optp.parse_args(my_arg_str.split())  # parse full set, assume initial argument already stripped

    return opts

parser = argparse.ArgumentParser()
parser.add_argument("--working-directory",default="./")
parser.add_argument("--approx",default=None,action='append',help="add multiple --approx options to analyze many approximants.  Every approximant is evaluated on the SAME intrinsic grid, and their likelihoods are marginalized over at each grid point; see RIFT/misc/DESIGN_multiapprox_marginalization.md")
parser.add_argument("--approx-prior",default=None,action='append',help="APPROX=WEIGHT prior weight for one waveform model (repeatable).  Default: uniform.  These are the p(m) in L_marg(lambda) = sum_m p(m) L_m(lambda); they are NOT sampling weights.")
parser.add_argument("--approx-gwsignal",default=None,action='append',help="Waveform model that must be generated through gwsignal rather than lalsim (repeatable).  The generator route is PER MODEL: the phenom family has no time-domain mode generator in gwsignal ('generator does not provide a method to generate time-domain modes'), while SEOBNRv5* is only available there.  A single global --use-gwsignal therefore cannot serve an EOB-vs-phenom comparison, which is the comparison this builder exists for.")
parser.add_argument("--require-all-approx",action='store_true',help="Drop intrinsic points not successfully evaluated under EVERY approximant.  Without it, such points are marginalized over whichever subset survived, which changes the estimator point by point.")
parser.add_argument("--input-grid",default="overlap-grid.xml.gz")
parser.add_argument("--cip-args",default=None,help="filename of args_cip.txt file  which holds CIP arguments.  Should NOT conflict with arguments auto-set by this DAG ... in particular, i/o arguments will be modified. ")
parser.add_argument("--cip-args-list",default=None,help="filename of args_cip_list.file, which holds CIP arguments.  Overrides cip-args if present. One CIP_n.sub file is created for each line in the file, which is used for an integer m iterations, where m is the first item of each line (normally 'X' in CIP)")
parser.add_argument("--cip-explode-jobs",default=None,type=int,help="Number of CIP jobs to use to use in parallel to produce posterior samples for each iteration. Code will generate a fit first, save it, use this number of workers in parallel to generate samples, and then consolidates the samples together. Note that --n-output-samples and --n-eff are NOT adjusted ... if the user wants to adaptively fix the resolution, that needs to be controlled at a higher level")
parser.add_argument("--cip-explode-jobs-flat",action='store_true',help="Pass to use the same arguments for all worker jobs. The main job will be /bin/true.")
parser.add_argument("--puff-exe",default=None,help="util_ParameterPuffball.py")
parser.add_argument("--puff-args",default=None,help="util_ParameterPuffball arguments.  If not specified, puffball will not be performed  ")
parser.add_argument("--puff-cadence",default=None,type=int,help="Every n iterations (not including 0), the puffball code will be applied.  Puffball points will be done *in addition* to the usual results from the DAG.  (The puffball is based on perturbing points from that iteration, and this will roughly double that iteration in ILE job size).  Proposed value 2 (i.e., puff overlap-grid-2, ...-4, ...-6.  If not specified, puffball will not be performed ")
parser.add_argument("--puff-max-it",default=-1,type=int,help="Maximum iteration number that puffball is applied.  If negative, puffball is not applied ")
parser.add_argument("--last-iteration-extrinsic",action='store_true',help="Configure last iteration to extract *one* set of extrinsic parameters from each intrinsic point. [This is highly inefficient, but people like having one extrinsic point per intrinsic point.]  Requires --convert-args")
parser.add_argument("--last-iteration-extrinsic-nsamples",default=3000,type=int,help="Construct this number of extrinsic samples")
parser.add_argument("--last-iteration-extrinsic-samples-per-ile",default=5,type=int,help="Fair-draw this many extrinsic posterior samples from each terminal ILE result")
parser.add_argument("--last-iteration-extrinsic-samples-per-ile-internal",default=10,type=int,help="Minimum terminal ILE n_eff used to support the per-result fair draw")
parser.add_argument("--last-iteration-extrinsic-time-resampling",action='store_true',help="Recover posterior times when the terminal ILE uses time marginalization")
parser.add_argument("--last-iteration-export-marginal-distance-grid", action='store_true', help="Add argument to ILE_extr")
parser.add_argument("--last-iteration-export-distance-slices", default=0, type=int, help="If >0, the ILE_extr (extrinsic) stage exports K-row .dslice files: Plan-B fixed-distance extrinsic-marginalized likelihoods. Adds --export-distance-slices K (+ --internal-use-lnL) to ILE_extr and strips --distance-marginalization.")
parser.add_argument("--last-iteration-export-distance-slices-n-core", default=0, type=int, help="Passthrough to ILE --n-distance-slice-core for the extrinsic-stage .dslice export (0 = ILE default).")
parser.add_argument("--last-iteration-export-distance-slices-n-wing", default=0, type=int, help="Passthrough to ILE --n-distance-slice-wing for the extrinsic-stage .dslice export (0 = ILE default).")
parser.add_argument("--last-iteration-export-distance-slices-wing-delta-lnL", default=None, type=float, help="Passthrough to ILE --distance-slice-wing-delta-lnL for the extrinsic-stage .dslice export.")
parser.add_argument("--last-iteration-export-distance-slices-skip-threshold", default=None, type=float, help="Passthrough to ILE --distance-slice-skip-threshold for the extrinsic-stage .dslice export.")
parser.add_argument("--ile-args",default=None,help="filename of args_ile.txt file  which holds ILE arguments.  Should NOT conflict with arguments auto-set by this DAG ... in particular, i/o arguments will be modified")
parser.add_argument("--ile-exe",default=None,help="filename of ILE or equivalent executable. Will default to `which integrate_likelihood_extrinsic` in low-level code")
parser.add_argument("--ile-retries",default=0,type=int,help="Number of retry attempts for ILE jobs. (These can fail)")
parser.add_argument("--general-retries",default=0,type=int,help="Number of retry attempts for internal jobs (convert, CIP, ...). (These can fail, albeit more rarely, usually due to filesystem problems)")
parser.add_argument("--ile-n-events-to-analyze",default=1,type=int,help="If >1, you are using ILE_batchmode.  Structures the DAG correctly to account for batch cadence")
parser.add_argument("--ile-runtime-max-minutes",default=None,type=int,help="If not none, kills ILE jobs that take longer than the specified integer number of minutes. Do not use unless an expert")
parser.add_argument("--cip-exe",default=None,help="filename of CIP or equivalent executable. Will default to `which util_ConstructIntrinsicPosterior_GenericCoordinates` in low-level code")
parser.add_argument("--test-exe",default=None,help="filename of test code or equivalent executable.  Must have a --test-output <fname> argument.  Used for convergence testing or other termination. NOT ACTIVE; see 'convergence_test_samples.py' for example")
parser.add_argument("--plot-exe",default=None,help="filename of plot code or equivalent executable.  Will default to `which plot_posterior_corner.py`.  Default is to plot last set of samples")
parser.add_argument("--gridinit-exe",default=None,help="filename of initial grid creation or equivalent executable.  Will default to `which util_ManualOverlapGrid.py`.  Must create overlap-grid-0.xml.gz")
parser.add_argument("--bw-exe",default=None,help="BayesWave location (or wrapper)")
parser.add_argument("--bw-post-exe",default=None,help="BayesWave location (or wrapper)")
parser.add_argument("--test-args",default=None,help="filename of args_test.txt, which holds test arguments.  Note i/o arguments will be modified, so should NOT specify the samples files or the output file, just the test to be performed and any related arguments")
parser.add_argument("--convert-args",default=None,help="filename of args_convert.txt, which holds arguments to the convert function. Needed if you plan to export tide informaton")
parser.add_argument("--plot-args",default=None,help="filename of args_plot.txt, which holds plot arguments.  Note i/o arguments will be modified, so should NOT specify the samples files or the output file, just the test to be performed and any related arguments.  To create the posterior samples file, you probably need to run the convert job (e.g., via the tests)")
parser.add_argument("--gridinit-args",default=None,help="arguments for grid creation or equivalent executable.  If empty, will attempt to use the --input-grid argument")
parser.add_argument("--ile-batch",action='store_true',help="Different workflow: the ILE job is assumed to create a .composite file directly in the current working directory. Designed for fake ILE sets and long-term batch workflow. NOT IMPLEMENTED ")
parser.add_argument("--transfer-file-list",default=None,help="File containing list of *input* filenames to transfer, one name per file. Copied into transfer_files for condor directly.  If provided, also enables attempts to deduce files that need to be transferred for the pipeline to operate, as needed for OSG, etc")
parser.add_argument("--request-memory-ILE",default=4096,type=int,help="Memory request for condor (in Mb) for ILE jobs.")
parser.add_argument("--request-gpu-ILE",action='store_true',help="ILE will request a GPU. [Note: if you do this, your code will only run on GPU-enabled slots]")
parser.add_argument("--request-memory-CIP",default=16384,type=int,help="Memory request for condor (in Mb) for fitting jobs.")
parser.add_argument("--use-singularity",action='store_true',help="Attempts to use a singularity image in SINGULARITY_RIFT_IMAGE")
parser.add_argument("--use-osg",action='store_true',help="Attempts to set up an OSG workflow.  Must submit from an osg allowed submit machine")
parser.add_argument("--use-osg-simple-requirements",action='store_true',help="Uses an aggressive simplified requirements string for worker jobs")
parser.add_argument("--condor-local-nonworker",action='store_true',help="Uses local universe for non-worker condor jobs. Important to run in non-NFS location, as other jobs don't have file transfer set up.")
parser.add_argument("--frames-dir",default=None,help="If you are using synthetic dat")
parser.add_argument("--cache-file",default=None,help="If you are using real data and have CVMS frames, you want to transfer the cache file over.")
parser.add_argument("--use-cvmfs-frames",action='store_true',help="If true, require LIGO frames are present (usually via CVMFS). User is responsible for generating cache file compatible with it.")
parser.add_argument('--n-copies',default=2,type=int,help="Number of duplicates of each ILE job, for redundant Monte Carlo")
parser.add_argument('--n-iterations',default=3,type=int,help="Number of iterations to perform")
parser.add_argument("--start-iteration",default=0,type=int,help="starting iteration. If >0, does not copy over --input-grid. DOES rewrite sub files.  This allows you to change the arguments provided (e.g., use more iterations or settings at late times). Note this overwrites the .sub files ")
parser.add_argument('--n-samples-per-job',default=2000,type=int,help="Number of samples generated each iteration; also, number of ILE jobs run each iteration. Should increase with dimension of problem")
parser.add_argument('--neff-threshold',default=800,type=int,help="Number of samples generated each iteration")
parser.add_argument("--workflow",default='single',help="[single|fit+posterior|full] describes workflow layout used.  'Single' is a single node, running the fit and posterior for each iteration; 'full' produces many followup jobs to produce a reliable posterior")
parser.add_argument("--n-post-jobs",default=1,type=int,help="Number of posterior jobs. Used in posterior and fit+posterior workflows")
parser.add_argument("--use-bw-psd",action='store_true',help="Use BW PSD, attempting to use fiducial arguments as in LI for placement (i.e., signal at seglen -2). Assumes LI style data convention.  Necessary BW will be parsed out of ile-args.txt (required)")
# ---------------------------------------------------------------------------
# API compatibility with util_RIFT_pseudo_pipe.py.
#
# pseudo_pipe emits ONE command line, shaped for BasicIteration's option surface,
# which is far wider than this builder's: this workflow implements cross-model
# marginalization, not every feature of the single-model pipeline.
#
# These options are ACCEPTED so a pseudo_pipe command line parses, and REFUSED
# below if one is actually used.  They are never silently ignored -- a flag that
# parses and does nothing changes a run without saying so, and is invisible in a
# config diff.  test_multiapprox_pseudo_pipe.py re-derives these lists from the
# two builders' argparse surfaces so they cannot drift.
_API_ONLY_FLAGS = ['--calibration-reweighting', '--calibration-reweighting-osg', '--calmarg-pilot', '--cip-explode-jobs-dag', '--cip-explode-jobs-subdag', '--comov-distance-reweighting', '--condor-containerize-nonworker', '--condor-local-nonworker-igwn-prefix', '--condor-nogrid-nonworker', '--extrinsic-handoff', '--first-iteration-jumpstart', '--frame-rotation', '--ile-group-subdag', '--ile-group-subdag-check-work', '--last-iteration-export-distance-slices-all-fresh', '--last-iteration-export-distance-slices-randomize', '--last-iteration-extrinsic-batched-convert', '--search-reflected-sky-mode', '--use-eccentricity', '--use-eccentricity-squared-sampling', '--use-full-submit-paths', '--use-hyperbolic', '--use-osg-cip', '--use-tabular-eos-file']
_API_ONLY_VALUED = ['--bilby-ini-file', '--bilby-pickle-exe', '--bilby-pickle-file', '--cal-request-disk', '--calibration-reweighting-batchsize', '--calibration-reweighting-count', '--calibration-reweighting-exe', '--calibration-reweighting-extra-args', '--calibration-reweighting-initial-extra-args', '--calmarg-pilot-cadence', '--calmarg-pilot-max-it', '--calmarg-pilot-max-points', '--calmarg-pilot-top-fraction', '--cip-explode-jobs-last', '--cip-post-exe', '--cip-request-disk', '--comov-distance-reweighting-exe', '--convert-ascii2h5-exe', '--extrinsic-handoff-select', '--fetch-ext-grid-args', '--fetch-ext-grid-exe', '--general-request-disk', '--ile-condor-commands', '--ile-gpu-fanout', '--ile-n-events-to-analyze-first', '--ile-post-exe', '--ile-request-disk', '--last-iteration-export-distance-slices-wing-neff', '--last-iteration-export-distance-slices-wing-nmax', '--lisa-reference-time', '--n-eff', '--n-iterations-subdag-max', '--n-samples-per-job-threshold', '--reflected-sky-mode-exe', '--search-reflected-sky-mode-iteration', '--use-oauth-files']
for _flag in _API_ONLY_FLAGS:
    parser.add_argument(_flag, action='store_true', help=argparse.SUPPRESS)
for _flag in _API_ONLY_VALUED:
    parser.add_argument(_flag, default=None, help=argparse.SUPPRESS)
# ---------------------------------------------------------------------------
opts=  parser.parse_args()

# Refuse the API-only options if any was actually requested (see above).
_requested = [f for f in _API_ONLY_FLAGS if getattr(opts, f[2:].replace('-', '_'), False)]
_requested += [f for f in _API_ONLY_VALUED if getattr(opts, f[2:].replace('-', '_'), None) is not None]
if _requested:
    print(" create_event_parameter_pipeline_BasicMultiApproxIteration accepts the")
    print(" following for API compatibility with util_RIFT_pseudo_pipe.py, but does")
    print(" NOT implement them: " + ", ".join(sorted(_requested)))
    print(" Refusing rather than ignoring: a flag that parses and does nothing")
    print(" changes a run without saying so.  Run those settings single-model, or")
    print(" implement the option here.")
    sys.exit(1)

if not opts.approx or len(opts.approx) < 2:
    parser.error("the multi-approximant builder requires at least two --approx values")
if len(set(opts.approx)) != len(opts.approx):
    parser.error("duplicate --approx values are not allowed: {}".format(opts.approx))
if opts.approx_prior:
    _prior_labels = set()
    _prior_total = 0.0
    for _item in opts.approx_prior:
        if "=" not in _item:
            parser.error("--approx-prior wants APPROX=WEIGHT, got {}".format(_item))
        _label, _value = [x.strip() for x in _item.split("=", 1)]
        if _label in _prior_labels:
            parser.error("duplicate --approx-prior for {}".format(_label))
        try:
            _weight = float(_value)
        except ValueError:
            parser.error("--approx-prior weight is not numeric: {}".format(_item))
        if not np.isfinite(_weight) or _weight < 0:
            parser.error("--approx-prior weights must be finite and non-negative")
        _prior_labels.add(_label)
        _prior_total += _weight
    if _prior_total <= 0:
        parser.error("--approx-prior weights must have positive total mass")



local_worker_universe="vanilla"
if opts.condor_local_nonworker:
    local_worker_universe="local"

working_dir_inside = opts.working_directory
if opts.use_singularity or opts.use_osg:
    working_dir_inside = "./" # all files on the remote machine are in the current directory

singularity_image = None
if opts.use_singularity:
    print(" === USING SINGULARITY === ")
    singularity_image = os.environ["SINGULARITY_RIFT_IMAGE"]  # must be present to use singularity
    # SINGULARITY IMAGES ARE ON CVMFS, SO WE CAN AVOID THE SINGULARITY EXEC CALL
    # hardcoding a fiducial copy of lalapps_path2cache; beware about the executable name change
    os.environ['LALAPPS_PATH2CACHE'] = "/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py39/bin/lalapps_path2cache" #"singularity exec {singularity_image} lalapps_path2cache".format(singularity_image=singularity_image)
    print(singularity_image)

if (opts.cip_args is None) and (opts.cip_args_list is None):
    print(" No arguments provided for low-level job")
    sys.exit(0)

# Load args.txt. Remove first item.  Store
if not (opts.cip_args is None):
    with open(opts.cip_args) as f:
        cip_args_list = f.readlines()
    cip_args = ' '.join( [x.replace('\\','') for x in cip_args_list] )
    cip_args = ' '.join(cip_args.split(' ')[1:])
    # Some argument protection for later
    cip_args = cip_args.replace('[', ' \'[')
    cip_args = cip_args.replace(']', ']\'')
    cip_args=cip_args.rstrip()
    cip_args += ' --no-plots '  
    print("CIP", cip_args)

cip_args_lines = None
cip_args_prefixes = []  # so it works correctly even in flat mode
if not (opts.cip_args_list is None):
    with open(opts.cip_args_list) as f:
        cip_args_lines = f.readlines()
    # Pull off the leading token. It is usually an integer (number of iterations), but util_RIFT_pseudo_pipe.py
    # also emits 'Z' (run-to-convergence marker) and 'G' (grid-stage marker) prefixes; tolerate those here.
    cip_args_prefixes = [(x.split(' ')[0]) for x in cip_args_lines]
    cip_args_n = (cip_args_prefixes.copy())
    for indx in np.arange(len(cip_args_n)):
        if cip_args_prefixes[indx][0]=='Z':
            cip_args_n[indx] = 1  # one nominal iteration; this workflow has no run-to-convergence CIP stage
        elif cip_args_prefixes[indx][0] =='G':
            cip_args_n[indx] = int(cip_args_prefixes[indx][1:])  # integer after G assumed
        else:
            cip_args_n[indx]  = int(cip_args_n[indx])
    cip_args_lines = [' '.join(x.split(' ')[1:]) for x in cip_args_lines] # pull off the integer
    cip_args_lines = [x.replace('[', ' \'[').replace(']', ']\'').rstrip() for x in cip_args_lines]
    cip_args_lines = [x.lstrip() for x in cip_args_lines] # remove leading whitespace
    for line in cip_args_lines:
        print(" CIP ", line)

    if np.sum(cip_args_n) < opts.n_iterations:
        print(" Not enough CIP versions to complete all requested iterations; extending ")
        cip_args_n[-1]  = opts.n_iterations - np.sum(cip_args_n)

# Load args.txt. Remove first item.  Store
with open(opts.ile_args) as f:
    ile_args_list = f.readlines()
ile_args = ' '.join( [x.replace('\\','') for x in ile_args_list] )
ile_args = ' '.join(ile_args.split(' ')[1:])
# Some argument protection for later
ile_args = ile_args.replace('[', ' \'[')
ile_args = ile_args.replace(']', ']\'')
ile_args=ile_args.rstrip()
if opts.request_gpu_ILE:
    ile_args+=" --vectorized --gpu "  # append this
if opts.ile_n_events_to_analyze > 1:
    ile_args+= " --n-events-to-analyze " + str(opts.ile_n_events_to_analyze)
#ile_args += ' --soft-fail-event-range '  #important in case event is out of range. Simply exit NO ...not backwards compatible
print("ILE", ile_args)

puff_args=None
puff_cadence = None
puff_max_it = opts.puff_max_it
if opts.puff_args and opts.puff_cadence:
    puff_cadence = opts.puff_cadence
    # Load args.txt. Remove first item.  Store
    with open(opts.puff_args) as f:
        puff_args_list = f.readlines()
    puff_args = ' '.join( [x.replace('\\','') for x in puff_args_list] )
    puff_args = ' '.join(puff_args.split(' ')[1:])
# Some argument protection for later
    puff_args = puff_args.replace('[', ' \'[')
    puff_args = puff_args.replace(']', ']\'')
    puff_args=puff_args.rstrip()
    print("PUFF", puff_args)
    print("PUFF CADENCE", puff_cadence)


test_args=None
if opts.test_args:
    # Load args.txt. Remove first item.  Store
    with open(opts.test_args) as f:
        test_args_list = f.readlines()
    test_args = ' '.join( [x.replace('\\','') for x in test_args_list] )
    test_args = ' '.join(test_args.split(' ')[1:])
# Some argument protection for later
    test_args = test_args.replace('[', ' \'[')
    test_args = test_args.replace(']', ']\'')
    test_args=test_args.rstrip()
    # A convergence success is reported as exit status 1 so DAGMan can stop
    # scheduling later intrinsic iterations.  That is incompatible with this
    # builder's fixed terminal fork: the per-model evidence and extrinsic
    # products are deliberately created only after the configured iteration
    # loop.  util_RIFT_pseudo_pipe.py already adds --always-succeed whenever
    # --add-extrinsic is requested; enforce the same contract here so direct
    # builder callers cannot make the terminal branch unreachable.
    if (opts.last_iteration_extrinsic and
            re.search(r'(?<!\S)--always-succeed(?=\s|$)', test_args) is None):
        print(" Terminal extraction: disabling early convergence termination")
        test_args += " --always-succeed "
    print("CONVERGE", test_args)

plot_args=None
if opts.plot_args:
    # Load args.txt. Remove first item.  Store
    with open(opts.plot_args) as f:
        plot_args_list = f.readlines()
    plot_args = ' '.join( [x.replace('\\','') for x in plot_args_list] )
    plot_args = ' '.join(plot_args.split(' ')[1:])
# Some argument protection for later
    plot_args = plot_args.replace('[', ' \'[')
    plot_args = plot_args.replace(']', ']\'')
    plot_args=plot_args.rstrip()
    print("PLOT", plot_args)

convert_args=None
if opts.convert_args:
    # Load args.txt. Remove first item.  Store
    with open(opts.convert_args) as f:
        convert_args_list = f.readlines()
    convert_args = ' '.join( [x.replace('\\','') for x in convert_args_list] )
    convert_args = ' '.join(convert_args.split(' ')[1:])
# Some argument protection for later
    convert_args = convert_args.replace('[', ' \'[')
    convert_args = convert_args.replace(']', ']\'')
    convert_args=convert_args.rstrip()
    print("CONVERT", convert_args)

gridinit_args=None
if opts.gridinit_args:
    # Load args.txt. Remove first item.  Store
    with open(opts.gridinit_args) as f:
        gridinit_args_list = f.readlines()
    gridinit_args = ' '.join( [x.replace('\\','') for x in gridinit_args_list] )
    gridinit_args = ' '.join(gridinit_args.split(' ')[1:])
# Some argument protection for later
    gridinit_args = gridinit_args.replace('[', ' \'[')
    gridinit_args = gridinit_args.replace(']', ']\'')
    gridinit_args=gridinit_args.rstrip()
    print("GRIDINIT", gridinit_args)


# Copy seed grid into place as overlap-grid-0.xml.gz
it_start = opts.start_iteration
n_initial = opts.n_samples_per_job
if (it_start == 0) and not gridinit_args:
    shutil.copyfile(opts.input_grid,"overlap-grid-0.xml.gz")  # put in working directory !
    n_initial = len(lalsimutils.xml_to_ChooseWaveformParams_array("overlap-grid-0.xml.gz"))

transfer_file_names = []
if not (opts.transfer_file_list is None):
    transfer_file_names=[]
    # Load args.txt. Remove first item.  Store
    with open(opts.transfer_file_list) as f:
        for  line in f.readlines():
            transfer_file_names.append(line.rstrip())
    print(" Input files to transfer to job working directory (note!)", transfer_file_names)

###
### Fiducial fit job (=sanity check that code will run)
###
if (opts.cip_args is None):
    cip_args = cip_args_lines[0]
cmdname="%s/command-single_fit.sh" % opts.working_directory
cmd = open(cmdname, 'w')
arg_list = cip_args
exe = dag_utils.which("util_ConstructIntrinsicPosterior_GenericCoordinates.py")
cmd.write('#!/usr/bin/env bash\n')
cmd.close()
st = os.stat(cmdname)
import stat
os.chmod(cmdname, st.st_mode | stat.S_IEXEC)


cmdname="%s/command-single.sh" % opts.working_directory
cmd = open(cmdname, 'w')
arg_list = ile_args
exe = dag_utils.which("integrate_likelihood_extrinsic")
if opts.ile_n_events_to_analyze:
    exe = dag_utils.which("integrate_likelihood_extrinsic_batchmode")
cmd.write('#!/usr/bin/env bash\n')
cmd.write('# (run me only in working directories)\n')
cmd.write(exe + ' ' + arg_list + " --sim-xml overlap-grid-0.xml.gz")  # make it concrete, so I can test it
cmd.close()
st = os.stat(cmdname)
import stat
os.chmod(cmdname, st.st_mode | stat.S_IEXEC)
# Add argument to ile_args.txt to 
#    identify overlap grid input
#    identify output file names (?)
ile_args_orig = ile_args  # provides ability to strip out the output and replace it with alternate
ile_args+= ' --sim-xml ' + working_dir_inside + '/overlap-grid-$(macroiteration).xml.gz '


###
### DAG generation
###

dag = pipeline.CondorDAG(log=os.getcwd())


# Make PSD, if requested
if opts.use_bw_psd:
    print(" ===> Adding BW PSD to pipeline <=== ")
    opts_ile = parse_ile_args_for_bw(ile_args_orig)

    channel_names = {}
    for inst, chan in [c.split("=") for c in opts_ile.channel_name]:
        channel_names[inst] = chan

    flow_ifo_dict = {}
    if opts_ile.fmin_ifo:
        for inst, freq_str in [c.split("=") for c in opts_ile.fmin_ifo]:
           flow_ifo_dict[inst]=float(freq_str) 

    
    channel_dict_bw = {}
    for ifo in list(channel_names.keys()):
        fmin_here= opts_ile.fmin_template
        if ifo in list(flow_ifo_dict.keys()):
            fmin_here  = flow_ifo_dict[ifo]
        channel_dict_bw[ifo]= [channel_names[ifo], fmin_here]
    
    # make BW directory
    bw_dir = opts.working_directory+"/make_bw_psds"
    mkdir(bw_dir)
    for ifo in list(channel_names.keys()):
        mkdir(bw_dir+"/"+ifo)

    # seglen computation
    seglen  = opts_ile.data_end_time - opts_ile.data_start_time

    # Write bw sub files
    bw_job, bw_job_name = dag_utils.write_psd_sub_BW_monoblock(cache_file=opts_ile.cache_file,psd_length=seglen,srate=opts_ile.srate,event_time=opts_ile.event_time,log_dir=bw_dir+"/",exe=opts.bw_exe)
    bw_job.add_condor_cmd("initialdir",bw_dir+"/$(ifo)")
    bw_job.write_sub_file()

    # bw_job, bw_job_name = dag_utils.write_psd_sub_BW_step0(channel_dict=channel_dict_bw,cache_file=opts_ile.cache_file,psd_length=seglen,srate=opts_ile.srate,event_time=opts_ile.event_time,log_dir=bw_dir+"/",exe=opts.bw_exe)
    # bw2_job, bw2_job_name=dag_utils.write_psd_sub_BW_step1(channel_dict=channel_dict_bw,cache_file=opts_ile.cache_file,psd_length=seglen,srate=opts_ile.srate,event_time=opts_ile.event_time,log_dir=bw_dir+"/",exe=opts.bw_post_exe)
    # bw_job.add_condor_cmd("initialdir",bw_dir)
    # bw2_job.add_condor_cmd("initialdir",bw_dir)
    # bw2_job.write_sub_file()

    # write convert command
#    event_time_trunc ='%.2f' % float(opts_ile.event_time)   # format so only 3 digits after decimal 
#    event_time_trunc += '0'  # add last character, always zero?
#    convert_psd_job, convert_psd_job_name = dag_utils.write_convertpsd_sub(ifo="$(ifo)",file_input=bw_dir+"/" + str(event_time_trunc) + "_$(ifo)-PSD.dat",log_dir=bw_dir+"/")
    convert_psd_job, convert_psd_job_name = dag_utils.write_convertpsd_sub(ifo="$(ifo)",file_input=bw_dir+"/$(ifo)/$(ifo)_fairdraw_psd.dat",log_dir=bw_dir+"/$(ifo)/")
    convert_psd_job.add_condor_cmd('initialdir',working_dir_inside)
    convert_psd_job.write_sub_file()



# Which generator route does each model need?  Empty string for lalsim models,
# " --use-gwsignal " for the ones that only exist there.  Bound as a macro on
# every node that carries macroapprox, so one DAG can mix families.
_gwsignal_set = set(opts.approx_gwsignal or [])
for _m in _gwsignal_set:
    if _m not in opts.approx:
        print(" --approx-gwsignal names {}, which is not among --approx {}".format(_m, opts.approx)); sys.exit(1)
def gwsignal_macro(approx):
    return " --use-gwsignal " if approx in _gwsignal_set else " "

# Make directories for all iterations.
#   The CIP stage is model-INDEPENDENT inside the loop -- one fit over the
#   cross-model marginalized net -- so its directory is not approximant-tagged.
#   Everything upstream of the combine (ILE, consolidate) is per model.
for indx in np.arange(it_start,opts.n_iterations+1):
  mkdir(opts.working_directory+"/iteration_"+str(indx)+"_cip")
  mkdir(opts.working_directory+"/iteration_"+str(indx)+"_cip/logs")
  for approx in opts.approx:
    ile_dir = opts.working_directory+"/approx_{}_iteration_".format(approx)+str(indx)+"_ile"
    cip_dir = opts.working_directory+"/approx_{}_iteration_".format(approx)+str(indx)+"_cip"
    consolidate_dir = opts.working_directory+"/approx_{}_iteration_".format(approx)+str(indx)+"_con"
#    convert_dir = opts.working_directory+"/iteration_"+str(indx)+"_change"
    mkdir(ile_dir); mkdir(ile_dir+"/logs")
    mkdir(cip_dir);  mkdir(cip_dir+"/logs")
    mkdir(consolidate_dir); mkdir(consolidate_dir+"/logs")
#    mkdir(change_dir); mkdir(change_dir+"/logs")

    if opts.test_args:
        # model-independent: the convergence test compares successive
        # marginalized posteriors, not per-model ones
        test_dir = opts.working_directory+"/iteration_"+str(indx)+"_test"
        mkdir(test_dir); mkdir(test_dir+'/logs')
        
    if opts.plot_args:  # Overkill: currently only making plots on last iteration
        # model-independent, like the plot node itself: it summarises the run
        # once, and its node binds only the iteration macros
        plot_dir = opts.working_directory+"/iteration_"+str(indx)+"_plot"
        mkdir(plot_dir); mkdir(plot_dir+'/logs')



# ++++
# Create workflow tasks
# ++++

##   ILE job

ile_exe =opts.ile_exe
if ile_exe is None:
    # This builder's terminal stage relies on batchmode's bounded fair-draw
    # option even when a worker analyzes one intrinsic point.  Falling back to
    # integrate_likelihood_extrinsic for n_group==1 produces a command line
    # containing --fairdraw-extrinsic-output-n-max, which that executable does
    # not accept.
    ile_exe = dag_utils.which("integrate_likelihood_extrinsic_batchmode")
output_file_names = None
if opts.use_singularity or opts.use_osg:
    # ILE runs from approx_<model>_iteration_<N>_ile, while the shared grid is
    # written in the workflow root as overlap-grid-N.xml.gz.  `input-grid-N`
    # is not produced anywhere in this workflow.  This typo is invisible in a
    # shared-filesystem pool but makes Condor reject the transfer on OSG.
    transfer_file_names.append("../overlap-grid-$(macroiteration).xml.gz")
    #output_file_names = ','.join(["CME_out-$(macroevent)-$(cluster)-$(process).xml_{0}_.dat".format(x) for x in np.arange(opts.ile_n_events_to_analyze)])
    #print "OUTPUT FILES ", output_file_names
# The generator route is per MODEL, not per run.  A global --use-gwsignal (as
# util_RIFT_pseudo_pipe.py emits) sends every model through gwsignal, and the
# phenom family cannot be generated there at all -- its ILE jobs die with
# "ValueError: Invalid Argument" from gen_modes and contribute ZERO rows, so the
# run silently degrades to single-model.  Strip it and bind it per node instead.
if opts.approx_gwsignal:
    ile_args = ile_args.replace("--use-gwsignal", " ")
ile_args += " --approx $(macroapprox) $(macrogwsignal) "
# Puff ILE is the same per-model evaluation routed over a different grid.  It
# must be derived after the global route is stripped and the model macros are
# appended; capturing ile_args_orig above silently baked the primary model and
# its generator route into every puff node.
ile_args_forpuff = ile_args.replace(
    working_dir_inside + '/overlap-grid-$(macroiteration).xml.gz',
    working_dir_inside + '/puffball-$(macroiteration).xml.gz')
ile_job, ile_job_name = dag_utils.write_ILE_sub_simple(tag='ILE',log_dir=None,arg_str=ile_args,output_file="CME_out.xml",ncopies=opts.n_copies,exe=ile_exe,transfer_files=transfer_file_names,transfer_output_files=output_file_names,request_memory=opts.request_memory_ILE,request_gpu=opts.request_gpu_ILE,use_singularity=opts.use_singularity,singularity_image=singularity_image,use_osg=opts.use_osg,use_simple_osg_requirements=opts.use_osg_simple_requirements,frames_dir=opts.frames_dir,cache_file=opts.cache_file,use_cvmfs_frames=opts.use_cvmfs_frames,max_runtime_minutes=opts.ile_runtime_max_minutes)
# Modify: create macro for iteration 
#   - added on a per-node basis
# Modify: add macro argument for overlap grid to be used (kept in top-level directory)
# Modify: set 'initialdir'
#     NOTE: logs will be specified relative to this directory
ile_job.add_condor_cmd("initialdir",opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile")
# Modify output argument: change logs and working directory to be subdirectory for the run
ile_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).log")
ile_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).err")
ile_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).out")
ile_job.write_sub_file()

if not (opts.puff_args is None):
    transfer_file_names_puff = list(transfer_file_names[:-1])
    transfer_file_names_puff.append(opts.working_directory+"/puffball-$(macroiteration).xml.gz") # could also use ../puffball, because relative is ok to initialdir
    # Write a version of the ILE job that uses puffball inputs ... IDENTICAL to code above for standard ILE case,just different argument for input
    # Yes, it would logically be simpler to make overlap-grid.xml.gz larger ... but I want to keep 'real samples' and 'puffed samples' seperate.
    ilePuff_job, ilePuff_job_name = dag_utils.write_ILE_sub_simple(tag='ILE_puff',log_dir=None,arg_str=ile_args_forpuff,output_file="CME_out.xml",simple_unique=True,ncopies=opts.n_copies,exe=ile_exe,transfer_files=transfer_file_names_puff,request_memory=opts.request_memory_ILE,request_gpu=opts.request_gpu_ILE,use_singularity=opts.use_singularity,singularity_image=singularity_image,use_osg=opts.use_osg,use_simple_osg_requirements=opts.use_osg_simple_requirements,frames_dir=opts.frames_dir,cache_file=opts.cache_file,use_cvmfs_frames=opts.use_cvmfs_frames,max_runtime_minutes=opts.ile_runtime_max_minutes)
    ilePuff_job.add_condor_cmd("initialdir",opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile")
    ilePuff_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).log")
    ilePuff_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).err")
    ilePuff_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILE-$(macroevent)-$(cluster)-$(process).out")
    ilePuff_job.write_sub_file()

if  (opts.last_iteration_extrinsic):
    n_points_per_ILE = opts.last_iteration_extrinsic_samples_per_ile
    if n_points_per_ILE < 1:
        parser.error("--last-iteration-extrinsic-samples-per-ile must be positive")
    if opts.last_iteration_extrinsic_nsamples < 1:
        parser.error("--last-iteration-extrinsic-nsamples must be positive")
    n_extrinsic_batches = int(np.ceil(
        opts.last_iteration_extrinsic_nsamples/
        (1.0*opts.ile_n_events_to_analyze)))
    # The terminal work grid must cover every intrinsic event scheduled below.
    # Round up to a complete batch when the requested count is not divisible by
    # --ile-n-events-to-analyze; the final combiner applies the public output
    # cap and uniqueness policy after these interior evaluations.
    terminal_intrinsic_count = (n_extrinsic_batches *
                                opts.ile_n_events_to_analyze)
    # Match the maintained RIFT terminal-stage contract: never lower the ILE
    # convergence target supplied by the science configuration, and retain
    # enough effective samples to support the requested fair draw.
    n_eff_last = max(2*n_points_per_ILE,
                     opts.last_iteration_extrinsic_samples_per_ile_internal)
    configured_neff = re.findall(
        r'(?:^|\s)--n-eff(?:=|\s+)([0-9]+)(?=\s|$)', ile_args_orig)
    if configured_neff:
        n_eff_last = max(n_eff_last, int(configured_neff[-1]))
    # ILE job with modified output format.  The saved record MUST be a bounded
    # fair draw from the extrinsic posterior, not the raw importance-sampling
    # cache.  A difficult point can reach n-max with n_eff~1; raw
    # --save-samples then writes millions of rows and the downstream resampler
    # must materialize all of them merely to select five.  More importantly,
    # that intermediate is not itself a posterior sample set.
    #
    # --save-samples is still the switch that serializes the record, while
    # --fairdraw-extrinsic-output changes that record to an equal-weight draw
    # before serialization and bounds it at n_points_per_ILE.
    ile_args_extr = (ile_args +
        " --save-P 0.01 --save-samples --fairdraw-extrinsic-output "
        "--fairdraw-extrinsic-output-n-max {} --n-eff {} ".format(
            n_points_per_ILE, n_eff_last))
    # A time-marginalized likelihood otherwise exports the fiducial time for
    # every row.  Recover posterior times on the same terminal fair-draw path.
    _uses_time_marginalization = re.search(
        r'(?<!\S)--time-marginalization(?=\s|$)', ile_args_extr) is not None
    if ((opts.last_iteration_extrinsic_time_resampling or
         _uses_time_marginalization) and
            "--resample-time-marginalization" not in ile_args_extr):
        ile_args_extr += " --resample-time-marginalization "
    # Fair-drawn terminal rows must contain sampled distances.  If distance is
    # still marginalized, ILE fills the output with its reference distance;
    # its time-resampling helper also explicitly rejects this combination.
    # BasicIteration makes the same transition for terminal extraction.
    if re.search(r'(?<!\S)--distance-marginalization(?=\s|$)', ile_args_extr):
        print(" Fair-draw extraction: disabling distance marginalization")
        ile_args_extr = re.sub(
            r'(?<!\S)--distance-marginalization(?=\s|$)', ' ', ile_args_extr)
    # The extrinsic stage runs AFTER the fork, on this model's own terminal CIP
    # posterior -- not on the shared marginalized grid the iteration loop used.
    ile_args_extr = ile_args_extr.replace(
        working_dir_inside + '/overlap-grid-$(macroiteration).xml.gz',
        working_dir_inside + '/approx_$(macroapprox)_overlap-grid-$(macroiteration).xml.gz')
    #  - note we *disable* --no-adapt-after-first (if present), so each point is independent (e.g., in sky location)
    ile_args_extr = ile_args_extr.replace('--no-adapt-after-first','')
    if opts.last_iteration_export_marginal_distance_grid:
        ile_args_extr += " --export-marginal-distance-grid "
        ile_args_extr = re.sub(
            r'(?<!\S)--distance-marginalization(?=\s|$)', ' ', ile_args_extr)
    if opts.last_iteration_export_distance_slices and opts.last_iteration_export_distance_slices > 0:
        ile_args_extr += " --export-distance-slices {} ".format(opts.last_iteration_export_distance_slices)
        if opts.last_iteration_export_distance_slices_n_core:
            ile_args_extr += " --n-distance-slice-core {} ".format(opts.last_iteration_export_distance_slices_n_core)
        if opts.last_iteration_export_distance_slices_n_wing:
            ile_args_extr += " --n-distance-slice-wing {} ".format(opts.last_iteration_export_distance_slices_n_wing)
        if opts.last_iteration_export_distance_slices_wing_delta_lnL is not None:
            ile_args_extr += " --distance-slice-wing-delta-lnL {} ".format(opts.last_iteration_export_distance_slices_wing_delta_lnL)
        if opts.last_iteration_export_distance_slices_skip_threshold is not None:
            ile_args_extr += " --distance-slice-skip-threshold {} ".format(opts.last_iteration_export_distance_slices_skip_threshold)
        if "--internal-use-lnL" not in ile_args_extr:
            ile_args_extr += " --internal-use-lnL "
        ile_args_extr = ile_args_extr.replace("--distance-marginalization ", ' ')
    transfer_file_names_extr = list(transfer_file_names)
    if opts.use_singularity or opts.use_osg:
        # The terminal fork reads a model-specific CIP grid, not the shared
        # in-loop grid.  Condor expands these macros in transfer_input_files
        # and stages the basename used by ile_args_extr in the remote sandbox.
        transfer_file_names_extr[-1] = (
            "../approx_$(macroapprox)_overlap-grid-$(macroiteration).xml.gz")
    ileExtr_job, ileExtr_job_name = dag_utils.write_ILE_sub_simple(
        tag='ILE_extr', log_dir=None, arg_str=ile_args_extr,
        output_file="EXTR_out.xml", simple_unique=True, ncopies=1,
        exe=ile_exe, transfer_files=transfer_file_names_extr,
        request_memory=opts.request_memory_ILE*2,
        request_gpu=opts.request_gpu_ILE,
        use_singularity=opts.use_singularity,
        singularity_image=singularity_image, use_osg=opts.use_osg,
        use_simple_osg_requirements=opts.use_osg_simple_requirements,
        frames_dir=opts.frames_dir, cache_file=opts.cache_file,
        use_cvmfs_frames=opts.use_cvmfs_frames,
        max_runtime_minutes=opts.ile_runtime_max_minutes)
    ileExtr_job.add_condor_cmd("initialdir",opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile")
    ileExtr_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILEextr-$(macroevent)-$(cluster)-$(process).log")
    ileExtr_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILEextr-$(macroevent)-$(cluster)-$(process).err")
    ileExtr_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/ILEextr-$(macroevent)-$(cluster)-$(process).out")
    ileExtr_job.write_sub_file()
    
    # Convert task
    convert_args_extr =  " --convention LI  --export-cosmology --use-interpolated-cosmology " 
    if not (convert_args is None):
        convert_args_extr += convert_args
    convertExtr_job, convertExtr_job_name = dag_utils.write_convert_sub(tag='convert_extr',log_dir=None,arg_str=convert_args_extr,file_input=opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/EXTR_out-$(macroevent).xml_$(macroindx)_.xml.gz",file_output=opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/EXTR_out-$(macroevent).xml_$(macroindx)_.dat", out_dir=opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/",universe=local_worker_universe)
    convertExtr_job.add_condor_cmd("initialdir",opts.working_directory)
    convertExtr_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/convert-$(macroevent)-$(macroindx).log")
    convertExtr_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/convert-$(macroevent)-$(macroindx).err")
    convertExtr_job.write_sub_file()

    # The converted rows are already equal-weight fair draws.  Do not pass them
    # through util_ResampleILEOutputWithExtrinsic: that would apply the
    # importance weights a second time as well as reintroducing an unbounded
    # full-table read.
    cat_job, cat_job_name = dag_utils.write_cat_sub(
        # Match only convert_extr's EXTR_out-..._.dat products.  The narrower
        # suffix also prevents stale *.downsampled_dat.dat files from an old
        # DAG revision being silently concatenated after a rescue/rebuild.
        file_prefix='EXTR', file_postfix='_.dat',
        search_root=opts.working_directory+
                    '/approx_$(macroapprox)_iteration_$(macroiteration)_ile',
        file_output=opts.working_directory+
                    '/extrinsic_posterior_samples_$(macroapprox).dat',
        universe=local_worker_universe,
        expected_batches=n_extrinsic_batches,
        events_per_batch=opts.ile_n_events_to_analyze)
    cat_job.add_condor_cmd("initialdir",opts.working_directory)
    cat_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/cat-$(cluster)-$(process).log")
    cat_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/cat-$(cluster)-$(process).out")
    cat_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile/logs/cat-$(cluster)-$(process).err")
    cat_job.write_sub_file()


##   Consolidate job(s)
#   - consolidate output of single previous job
con_job, con_job_name = dag_utils.write_consolidate_sub_simple(tag='join',log_dir=None,arg_str='',base=opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_ile", target=opts.working_directory+'/approx_$(macroapprox)_consolidated_$(macroiteration)',universe=local_worker_universe)
# Modify: set 'initialdir' : should run in top-level direcory
#con_job.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_con")
# Modify output argument: change logs and working directory to be subdirectory for the run
con_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_con/logs/con-$(cluster)-$(process).log")
con_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_con/logs/con-$(cluster)-$(process).err")
con_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_con/logs/con-$(cluster)-$(process).out")
con_job.write_sub_file()

##   Unify jobs
#
#   TWO nets, because the workflow needs two different things:
#
#     all.net                 -- every model's composites pooled and
#                                MARGINALIZED over models at each shared
#                                intrinsic point.  This is what the in-loop CIP
#                                fits, so grid refinement follows the
#                                model-marginalized posterior.
#     approx_<A>_all.net      -- one model's composites only.  This feeds the
#                                terminal per-model CIP, which yields that
#                                model's posterior and evidence.
#
#   The marginalization itself is util_CleanILE.py's job: it keys on the
#   intrinsic parameters, so the same lambda evaluated under several models
#   collapses to one entry, combined LINEARLY in L.  See
#   RIFT/misc/DESIGN_multiapprox_marginalization.md.
clean_model_args = " --model-group-regex 'approx_(.+?)_consolidated' "
# Coverage must be judged against the models this run CONFIGURED.  A model whose
# composites are all empty is skipped before its label is recorded, so without
# this the pooled net silently becomes a lower-dimensional mixture and even
# --require-all-models sees nothing wrong.  Observed live: one model's ILE jobs
# all failed and util_CleanILE reported "combination over 1 models" as ordinary
# status.
clean_model_args += " --expect-models '{}' ".format(",".join(opts.approx))
if opts.approx_prior:
    for item in opts.approx_prior:
        if "=" not in item:
            print(" --approx-prior wants APPROX=WEIGHT, got ", item); sys.exit(1)
        label = item.split("=")[0].strip()
        if label not in opts.approx:
            print(" --approx-prior names {}, which is not among --approx {}".format(label, opts.approx)); sys.exit(1)
        clean_model_args += " --model-prior '{}' ".format(item)
    missing = [a for a in opts.approx if a not in [i.split("=")[0].strip() for i in opts.approx_prior]]
    if missing:
        print(" --approx-prior given but missing weights for ", missing, "; specify every approximant or none"); sys.exit(1)
if opts.require_all_approx:
    clean_model_args += " --require-all-models "

unify_job, unify_job_name = dag_utils.write_unify_sub_simple(tag='unify',log_dir='',arg_str=clean_model_args, base=opts.working_directory, target=opts.working_directory+'/all.net',universe=local_worker_universe,script_name='unify.sh',fail_on_error=True)
unify_job.add_condor_cmd("initialdir",opts.working_directory)
unify_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/unify-$(cluster)-$(process).log")
unify_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/unify-$(cluster)-$(process).err")
unify_job.set_stdout_file(opts.working_directory+'/all.net')
unify_job.write_sub_file()

##   Per-model unify, for the terminal forked CIP.  No --model-group-regex:
##   within one model the evaluations ARE replicas of one quantity, so the flat
##   ntot-weighted pool is the correct combination.
unify_model_job, unify_model_job_name = dag_utils.write_unify_sub_simple(tag='unify_model',log_dir='',arg_str='', base=opts.working_directory, target=opts.working_directory+'/approx_$(macroapprox)_all.net',universe=local_worker_universe,script_name='unify_model.sh',glob_pattern='approx_$(macroapprox)_*.composite',fail_on_error=True)
unify_model_job.add_condor_cmd("initialdir",opts.working_directory)
unify_model_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_con/logs/unifymodel-$(cluster)-$(process).log")
unify_model_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_con/logs/unifymodel-$(cluster)-$(process).err")
unify_model_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_all.net")
unify_model_job.write_sub_file()




##   Fit job: default case
cip_args_base = cip_args
out_dir_base = opts.working_directory
cip_exe = opts.cip_exe
if not(opts.cip_explode_jobs is None):
    if not opts.cip_explode_jobs_flat:
        cip_args_base += " --fit-save-gp " + opts.working_directory + "/iteration_$(macroiteration)_cip/my_fit"
    else:
        cip_exe = "/bin/true"
    out_dir_base += "/iteration_$(macroiteration)_cip/"
cip_job, cip_job_name = dag_utils.write_CIP_sub(tag='CIP',log_dir=None,arg_str=cip_args_base,request_memory=opts.request_memory_CIP,input_net=opts.working_directory+'/all.net',output='overlap-grid-$(macroiterationnext)',out_dir=out_dir_base,exe=cip_exe,universe=local_worker_universe)
# Modify: set 'initialdir'
cip_job.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_cip")
# Modify output argument: change logs and working directory to be subdirectory for the run
cip_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).log")
cip_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).err")
cip_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).out")
cip_job.write_sub_file()
if not(opts.cip_explode_jobs is None):
    print(" Exploding stage 1, with  ",opts.cip_explode_jobs, " workers producing samples ")
    cip_args_base = cip_args_base.replace('fit-save-gp','fit-load-gp')
    cip_args_base = cip_args_base.replace('my_fit', 'my_fit.pkl')  # yes, asymmetric arguments
    cip_job_worker, cip_job_worker_name = dag_utils.write_CIP_sub(tag='CIP_worker',log_dir=None,arg_str=cip_args_base,request_memory=opts.request_memory_CIP,input_net=opts.working_directory+'/all.net',output='overlap-grid-$(macroiterationnext)-$(process)',out_dir=out_dir_base,exe=opts.cip_exe,ncopies=opts.cip_explode_jobs,universe=local_worker_universe)
    # Modify: set 'initialdir'
    cip_job_worker.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_cip")
    # Modify output argument: change logs and working directory to be subdirectory for the run
    cip_job_worker.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).log")
    cip_job_worker.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).err")
    cip_job_worker.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).out")
    cip_job_worker.write_sub_file()

    # Create worker join job
    join_cip_job,join_cip_job_name = dag_utils.write_joingrids_sub(input_pattern=opts.working_directory+"/iteration_$(macroiteration)_cip/output-grid-*.xml.gz",target_dir=opts.working_directory,output_base="overlap-grid-$(macroiterationnext)",n_explode=opts.cip_explode_jobs,log_dir=opts.working_directory+"/iteration_$(macroiteration)_cip/logs",universe=local_worker_universe)
    join_cip_job.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_cip")
    # Modify output argument: change logs and working directory to be subdirectory for the run
    join_cip_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/join-$(cluster)-$(process).log")
    join_cip_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/join-$(cluster)-$(process).err")
    join_cip_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/join-$(cluster)-$(process).out")
    join_cip_job.write_sub_file()



##   Terminal per-model CIP.
#
#   The iteration loop fits the MARGINALIZED net, because that is what should
#   steer grid placement.  The terminal stage instead FORKS: each model gets
#   its own fit over its own net, yielding that model's posterior and -- via
#   --fname-output-integral -- that model's evidence ln Z_m.  Those are the
#   deliverables, and they are what the final combination is weighted by.
#   Combining at the end rather than the middle keeps each model's posterior
#   inspectable, which is what a systematics study actually wants to look at.
# Which CIP arguments did the run actually FINISH on?
#
# `cip_args` is cip_args_lines[0], but the iteration loop advances through the
# list via cip_job_list, so with a multi-line --cip-args-list the last iteration
# uses a later entry.  Building the terminal fit from cip_args would give the
# per-model posterior and evidence different coordinates, fit method or
# resolution than the iteration that produced the grid they are fitting -- and
# the evidences are what the final mixture is weighted by.
cip_args_terminal = cip_args
if cip_args_lines is not None:
    _per_iteration = []
    for _i in np.arange(len(cip_args_lines)):
        _per_iteration += int(cip_args_n[_i])*[cip_args_lines[_i]]
    if _per_iteration:
        _last = min(int(opts.n_iterations), len(_per_iteration)) - 1
        cip_args_terminal = _per_iteration[_last]

cip_terminal_job = None
if opts.last_iteration_extrinsic:
    # The iteration fit's posterior size is an adaptation setting, not the
    # terminal work-grid contract.  In particular pseudo_pipe commonly leaves
    # `--n-output-samples 5000 --posterior-unique-draw` in the final CIP line,
    # while a science terminal stage schedules 20000 intrinsic events.  The
    # high event indices then soft-exit successfully as out of range and the
    # exact cat manifest can never be completed.  Interior duplicate points are
    # valid importance-work items; exact uniqueness belongs only to the final
    # exported posterior (util_CombineApproximantPosteriors.py).
    cip_args_terminal = re.sub(
        r'(?<!\S)--n-output-samples(?:=\S+|\s+\S+)', ' ',
        cip_args_terminal)
    cip_args_terminal = re.sub(
        r'(?<!\S)--posterior-unique-draw(?=\s|$)', ' ',
        cip_args_terminal)
    cip_args_terminal += " --n-output-samples {} ".format(
        terminal_intrinsic_count)
    cip_terminal_job, cip_terminal_job_name = dag_utils.write_CIP_sub(tag='CIP_terminal',log_dir=None,arg_str=cip_args_terminal,request_memory=opts.request_memory_CIP,input_net=opts.working_directory+'/approx_$(macroapprox)_all.net',output='approx_$(macroapprox)_overlap-grid-$(macroiterationnext)',out_dir=opts.working_directory,exe=opts.cip_exe,universe=local_worker_universe)
    cip_terminal_job.add_condor_cmd("initialdir",opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_cip")
    cip_terminal_job.set_log_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_cip/logs/cipterm-$(cluster)-$(process).log")
    cip_terminal_job.set_stderr_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_cip/logs/cipterm-$(cluster)-$(process).err")
    cip_terminal_job.set_stdout_file(opts.working_directory+"/approx_$(macroapprox)_iteration_$(macroiteration)_cip/logs/cipterm-$(cluster)-$(process).out")
    cip_terminal_job.write_sub_file()


##   Cross-model combination of the per-model extrinsic posteriors, weighted by
##   p(m) Z_m.  This is the one place the models are recombined AFTER the fork.
combine_job = None
if opts.last_iteration_extrinsic:
    _combined_sample_target = (n_extrinsic_batches *
                               opts.ile_n_events_to_analyze *
                               n_points_per_ILE)
    combine_args = " --output " + opts.working_directory + "/extrinsic_posterior_samples.dat "
    combine_args += " --n-output-samples {} ".format(_combined_sample_target)
    for approx in opts.approx:
        combine_args += " --model {}:{}/extrinsic_posterior_samples_{}.dat:{}/approx_{}_overlap-grid-$(macroiteration)+annotation.dat ".format(
            approx, opts.working_directory, approx, opts.working_directory, approx)
    if opts.approx_prior:
        for item in opts.approx_prior:
            combine_args += " --model-prior '{}' ".format(item)
    combine_job, combine_job_name = dag_utils.write_convert_sub(
        tag='combine_models', log_dir=None, arg_str=combine_args,
        file_input=None, file_output=None, out_dir=opts.working_directory,
        exe=which("util_CombineApproximantPosteriors.py"), universe=local_worker_universe)
    combine_job.add_condor_cmd("initialdir",opts.working_directory)
    combine_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/combine-$(cluster)-$(process).log")
    combine_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/combine-$(cluster)-$(process).err")
    combine_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/combine-$(cluster)-$(process).out")
    combine_job.write_sub_file()


##   puffball job: default case
if puff_args and puff_cadence:
    # The puffball perturbs the grid the CIP just wrote.  This read
    # 'input-grid-N.xml.gz', a name nothing in this workflow produces, so the
    # puffball could never run -- it dies with FileNotFoundError as soon as the
    # iteration count is high enough to reach it.  BasicIteration uses
    # overlap-grid-N, which is what CIP actually writes.
    puff_job, puff_job_name = dag_utils.write_puff_sub(tag='PUFF',log_dir=None,arg_str=puff_args,request_memory=opts.request_memory_ILE,input_net=opts.working_directory+'/overlap-grid-$(macroiterationnext).xml.gz',output=opts.working_directory+'/puffball-$(macroiterationnext)',out_dir=opts.working_directory,exe=opts.puff_exe,universe=local_worker_universe)
    # Modify: set 'initialdir' to CIP WORKING DIR 
    puff_job.add_condor_cmd("initialdir",opts.working_directory)
    # Modify output argument: change logs and working directory to be subdirectory for the run
    puff_job.set_log_file(opts.working_directory+"/puff-$(cluster)-$(process).log")
    puff_job.set_stderr_file(opts.working_directory+"/puff-$(cluster)-$(process).err")
    puff_job.set_stdout_file(opts.working_directory+"//puff-$(cluster)-$(process).out")
    puff_job.write_sub_file()



cip_job_list = None
if  (cip_args_lines is None):
    # Write the default cip_job into cip_job_list n times
    if opts.cip_explode_jobs is None:
        cip_job_list =opts.n_iterations* [cip_job ]
    else:
        cip_job_list =opts.n_iterations* [ [cip_job, cip_job_worker] ]
else:  
    # we have different cip jobs for different iteration numbers
    cip_job_list = []
    for indx in np.arange(len(cip_args_lines)):
        out_dir_base = opts.working_directory
        cip_args_extra = ""
        cip_args_truncate=''
        cip_exe = opts.cip_exe
        if not(opts.cip_explode_jobs is None):
            if not opts.cip_explode_jobs_flat:
                cip_args_extra += " --fit-save-gp " + opts.working_directory + "/iteration_$(macroiteration)_cip/my_fit"
            out_dir_base += "/iteration_$(macroiteration)_cip/"
            # set n_eff for primary job to be small. ONLY used for the primary non-worker job
            cip_args_truncate = " --n-eff 5 --n-max 10000 "  # cap n_eff and number of iterations for non-worker jobs. Nonzero to avoid accidental crashes.
        # Write the appropriate CIP jobs. [note only one CIP per iteration, so unique
        cip_job, cip_job_name = dag_utils.write_CIP_sub(tag='CIP_'+str(indx),log_dir=None,arg_str=cip_args_lines[indx]+cip_args_extra+cip_args_truncate,request_memory=opts.request_memory_CIP,input_net=opts.working_directory+'/all.net',output='overlap-grid-$(macroiterationnext)',out_dir=out_dir_base,exe=cip_exe,universe=local_worker_universe)
        # Modify: set 'initialdir'
        cip_job.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_cip")
        # Modify output argument: change logs and working directory to be subdirectory for the run
        cip_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).log")
        cip_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).err")
        cip_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cip-$(cluster)-$(process).out")
        cip_job.write_sub_file()


        if not(opts.cip_explode_jobs is None):
           print( " Exploding stage 2, with  ",opts.cip_explode_jobs, " workers producing samples ")
           cip_args_extra = cip_args_extra.replace('fit-save-gp','fit-load-gp')
           cip_args_extra = cip_args_extra.replace('my_fit','my_fit.pkl')
           cip_job_worker, cip_job_worker_name = dag_utils.write_CIP_sub(tag='CIP_worker'+str(indx),log_dir=None,arg_str=cip_args_lines[indx]+cip_args_extra,request_memory=opts.request_memory_CIP,input_net=opts.working_directory+'/all.net',output='overlap-grid-$(macroiterationnext)-$(process)',out_dir=out_dir_base,exe=opts.cip_exe,ncopies=opts.cip_explode_jobs,universe=local_worker_universe)
           # Modify: set 'initialdir'
           cip_job_worker.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_cip")
           # Modify output argument: change logs and working directory to be subdirectory for the run
           cip_job_worker.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cipworker-$(cluster)-$(process).log")
           cip_job_worker.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cipworker-$(cluster)-$(process).err")
           cip_job_worker.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/cipworker-$(cluster)-$(process).out")
           cip_job_worker.write_sub_file()


       # augment the CIP job list
        if opts.cip_explode_jobs is None:
            cip_job_list = cip_job_list + cip_args_n[indx]*[cip_job]
        else:
            print(" Exploding stage 3 ")
            cip_job_list = cip_job_list + cip_args_n[indx]*[[cip_job,cip_job_worker]]

##   Test job (terminate, convergence
if opts.test_args:
    ##  Convert job : make results accessible after every iteration.  (Only performed if the tests are active, to make my life easier)
    convert_job, convert_job_name =dag_utils.write_convert_sub(tag='convert',log_dir=None,arg_str=convert_args,file_input=opts.working_directory+'/overlap-grid-$(macroiteration).xml.gz', file_output=opts.working_directory+'/posterior_samples-$(macroiteration).dat' ,out_dir=opts.working_directory,exe=opts.test_exe,universe=local_worker_universe)
    convert_job.add_condor_cmd("initialdir",opts.working_directory)
    convert_job.set_log_file(opts.working_directory+"/iteration_$(macroiterationlast)_test/logs/convert-$(cluster)-$(process).log")
    convert_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiterationlast)_test/logs/convert-$(cluster)-$(process).err")
    convert_job.write_sub_file()


    test_job, test_job_name = dag_utils.write_test_sub(tag='test',log_dir=None,arg_str=test_args,samples_files=[ opts.working_directory+'/posterior_samples-$(macroiteration).dat', opts.working_directory+'/posterior_samples-$(macroiterationlast).dat'] ,out_dir=opts.working_directory,exe=opts.test_exe,universe=local_worker_universe)
    # Modify: set 'initialdir'
    test_job.add_condor_cmd("initialdir",opts.working_directory+"/iteration_$(macroiteration)_test")
    # Modify output argument: change logs and working directory to be subdirectory for the run
    test_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_test/logs/test-$(cluster)-$(process).log")
    test_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_test/logs/test-$(cluster)-$(process).err")
    test_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_test/logs/test-$(cluster)-$(process).out")
    test_job.write_sub_file()


if opts.plot_args and opts.test_args:
    # default: last two iterations
    #  ... unless there is extrinsic samples, then use those
    #  ... note this has to work with PEsummary as well, which only takes one set of samples
#    samples_files =['posterior_samples-$(macroiteration).dat','posterior_samples-$(macroiterationlast).dat'] 
#    if opts.last_iteration_extrinsic:
#        samples_files =[]
    # User will be responsible for passing argument strings.  We are not hardcoding in the argument format
    samples_files =[]
    plot_job, plot_job_name = dag_utils.write_plot_sub(tag='plot',log_dir=None,arg_str=plot_args,samples_files=samples_files,out_dir=opts.working_directory,exe=opts.plot_exe,universe=local_worker_universe)
    plot_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_plot/logs/test-$(cluster)-$(process).log")
    plot_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_plot/logs/test-$(cluster)-$(process).err")
    plot_job.set_stdout_file(opts.working_directory+"/iteration_$(macroiteration)_plot/logs/test-$(cluster)-$(process).out")
    plot_job.write_sub_file()


if opts.gridinit_args:
    gridinit_job, gridinit_job_name = dag_utils.write_init_sub(tag='grid',log_dir=None,arg_str=gridinit_args,out_dir=opts.working_directory,exe=opts.gridinit_exe,universe=local_worker_universe)
    gridinit_job.set_log_file(opts.working_directory+"/init-$(cluster)-$(process).log")
    gridinit_job.set_stderr_file(opts.working_directory+"init--$(cluster)-$(process).err")
    gridinit_job.set_stdout_file(opts.working_directory+"/init-$(cluster)-$(process).out")
    gridinit_job.write_sub_file()


##   Convert job(s)
#    - relocate input grid to correct location.  (MAKE TRIVIAL: do via output arguments in previous stage!)
##   Assessment job (no-op for now)




# ++++
# Create workflow 
# ++++

# Create workflow
#   - Create grid node as needed
#   - Loop over iterations
#      - if iteration0, use seed grid (should already be copied in place)
#      - if not iteration 0, grid should be in place (from previous stage)
#      - Loop over events, make ILE node per event
#      - create consolidate job, make it depend on all events in that iteration
#      - create fit job, make it depend on consolidate job
#    

parent_fit_node = None
test_node_list = []

if opts.gridinit_args:
   grid_node = pipeline.CondorDAGNode(gridinit_job) 
   dag.add_node(grid_node)
   parent_fit_node = grid_node  # this must happen before the first ILE jobs

convert_psd_node_list = []
if opts.use_bw_psd:
    print(" ===> Adding BW PSD nodes to dag <=== ")

#    bw1_node =pipeline.CondorDAGNode(bw2_job)
#    bw1_node.add_parent(bw0_node)
#    dag.add_node(bw1_node)

    # Add nodes to convert the PSD to the correct format
    # 
    for ifo in list(channel_names.keys()):
        bw0_node =pipeline.CondorDAGNode(bw_job)
        bw0_node.add_macro('ifo',ifo)
        # Create argument string
        bw_arg_str = ''
        channel_name, channel_flow = channel_dict_bw[ifo]
        bw_arg_str += " --ifo "+ifo
        bw_arg_str += " --"+ifo+"-channel "+ifo+":"+channel_name
        ifo_cache_name =  working_dir_inside +"/for_bw_"+ifo+".cache"  # bw requires single-IFO cache files, argh!!!
        ifo_char = ifo[0]
        os.system("grep ^"+ifo_char +  " " +opts_ile.cache_file+ " > " + ifo_cache_name)
        bw_arg_str += " --"+ifo+"-cache "+ifo_cache_name
        bw_arg_str += " --"+ifo+"-flow "+str(channel_flow)
        bw_arg_str += " --"+ifo+"-timeslide 0.0"
        bw0_node.add_macro('macroargument0',bw_arg_str)
        bw0_node.set_retry(opts.ile_retries)  # these fail all the time
        dag.add_node(bw0_node)

        # Need correct inheritance: this will be annoying b/c loop below does not allow for it
        # ...so do the PSD convert steps *serially* for now.  (The individual PSDs can be done in parallel)
        convert_psd_node =pipeline.CondorDAGNode(convert_psd_job)
        convert_psd_node.add_parent(bw0_node)
#        if not(parent_fit_node is None):
#            convert_psd_node.add_parent(parent_fit_node)
        parent_fit_node = convert_psd_node
        convert_psd_node.add_macro('ifo',ifo)
        convert_psd_node.set_retry(opts.ile_retries)  # these fail all the time
        dag.add_node(convert_psd_node)
        convert_psd_node_list.append(convert_psd_node)

n_group = opts.ile_n_events_to_analyze

# ONE iteration = every model's ILE on the SAME grid -> per-model consolidate
# -> one marginalized net -> ONE fit -> one new grid, which every model's ILE
# reads next iteration.  The models are coupled at the combine step, not run
# side by side: that coupling is the point of this builder.  See
# RIFT/misc/DESIGN_multiapprox_marginalization.md.
for it in np.arange(it_start,opts.n_iterations):
  con_nodes_this_iteration = []
  for approx in opts.approx:
    # Create consolidate job (per model: composites stay separated by model,
    # because the terminal per-model CIP needs them apart)
    con_node = pipeline.CondorDAGNode(con_job)
    con_node.add_macro("macroiteration",it)
    con_node.add_macro("macroapprox",approx)
    con_node.set_retry(opts.general_retries)

    # Create one node per job
    n_jobs_this_time = opts.n_samples_per_job
    if it ==it_start:
        n_jobs_this_time = n_initial
    indx_max = int(np.ceil((1.0*n_jobs_this_time)/n_group))
    for event in np.arange(indx_max): #np.arange(n_jobs_this_time):
        # Add task per ILE operation
        ile_node = pipeline.CondorDAGNode(ile_job)
#        ile_node.set_priority(JOB_PRIORITIES["ILE"])
        ile_node.set_retry(opts.ile_retries)
        ile_node.add_macro("macroevent", event*n_group)
        ile_node.add_macro("macroapprox",approx)
        ile_node.add_macro("macrogwsignal", gwsignal_macro(approx))
        ile_node.add_macro("macroiteration", it)
        if not(parent_fit_node is None):
            ile_node.add_parent(parent_fit_node)
        if it == it_start:
            for node in convert_psd_node_list:  # for every PSD conversion job, make sure PSD is present before we run the first iteration!
                ile_node.add_parent(node)
        con_node.add_parent(ile_node) # consolidate depends on all of the individual jobs
        dag.add_node(ile_node)

    if puff_args and puff_cadence:
     if it>0 and it <= puff_max_it  and (it-1)%puff_cadence ==0:  # we made a puffball last iteration, so run it through ILE now
        print(" ILE jobs for puffball on iteration ", it)
        for event in np.arange(indx_max):
            ile_node = pipeline.CondorDAGNode(ilePuff_job)  # only difference is here: uses puffball, which by construction is the same size/ perturbed points
            ile_node.set_retry(opts.ile_retries)
            ile_node.add_macro("macroevent", event*n_group)
            ile_node.add_macro("macroapprox",approx)
            ile_node.add_macro("macrogwsignal", gwsignal_macro(approx))
            ile_node.add_macro("macroiteration", it)
            if not(parent_fit_node is None):
                ile_node.add_parent(parent_fit_node)
            con_node.add_parent(ile_node) # consolidate depends on all of the individual jobs
            dag.add_node(ile_node)

    # add con job
    dag.add_node(con_node)
    con_nodes_this_iteration.append(con_node)

  # ---- everything below is per ITERATION, not per model --------------------
  # One unify: pools every model's composites and marginalizes over models at
  # each shared intrinsic point.  It must wait for ALL models to consolidate.
  if True:
    unify_node = pipeline.CondorDAGNode(unify_job)
    unify_node.add_macro("macroiteration",it)
    for node in con_nodes_this_iteration:
        unify_node.add_parent(node)
    unify_node.set_retry(opts.general_retries)
    dag.add_node(unify_node)

    # Create fit node, which depends on the marginalized net
    cip_job = cip_job_list[it]
    if isinstance(cip_job,list):
        cip_worker_job = cip_job[1]
        cip_job=cip_job[0]
    fit_node = pipeline.CondorDAGNode(cip_job)
    fit_node.add_macro("macroiteration", it)
    fit_node.add_macro("macroiterationnext", it+1)
    fit_node.set_category("CIP")
    fit_node.add_parent(unify_node)  # only fit if we have results from the previous iteration
    fit_node.set_retry(opts.general_retries)
    dag.add_node(fit_node) 
    parent_fit_node = fit_node

    if not(opts.cip_explode_jobs is None):
        print(" Exploding workers out ")
        # Create exploded worker job nodes
        worker_node =pipeline.CondorDAGNode(cip_worker_job)
        worker_node.add_macro("macroiteration", it)
        worker_node.add_macro("macroiterationnext", it+1)
        worker_node.set_category("CIP_worker")
        worker_node.add_parent(parent_fit_node)  # only fit if we have results from the previous iteration
        worker_node.set_retry(opts.general_retries)

        # Create job to consolidate worker outputs
        join_node =pipeline.CondorDAGNode(join_cip_job)
        join_node.add_macro("macroiteration", it)
        join_node.add_macro("macroiterationnext", it+1)
        join_node.set_category("join_cip")
        join_node.add_parent(worker_node)
        join_node.set_retry(opts.general_retries)
        
        dag.add_node(worker_node)
        dag.add_node(join_node)
        parent_fit_node=join_node

    # Check if puffball being created, and if so create it *immediately* after CIP job creates the overlap-grid file
    if puff_args and puff_cadence:
      if it > -1 and it <= puff_max_it and it%puff_cadence ==0:
        print(" Puffball for iteration ", it)
        puff_node = pipeline.CondorDAGNode(puff_job)
        puff_node.add_macro("macroiteration", it)
        puff_node.add_macro("macroiterationnext", it+1)
        puff_node.set_category("PUFF")
        puff_node.set_retry(opts.general_retries)
        if not (parent_fit_node is None):
            puff_node.add_parent(parent_fit_node)  # only fit if we have results from the previous iteration
        dag.add_node(puff_node) 
        
        parent_fit_node = puff_node

        

    # Create convert node, which depends on fit node, *if* tests are being performed
    if opts.test_args:
        convert_node=pipeline.CondorDAGNode(convert_job)
        convert_node.add_macro("macroiteration", it+1)  # convert the NEWLY-PRODUCED iteration
        convert_node.add_macro("macroiterationlast", it)  # use log files in the previous directory
        convert_node.add_parent(parent_fit_node)
        convert_node.set_category("CONVERT")
        convert_node.set_retry(opts.general_retries)
        dag.add_node(convert_node)

        parent_fit_node = convert_node



    if opts.test_args and it>0:
        # Cannot run test on first iteration
        test_node = pipeline.CondorDAGNode(test_job)
        test_node.add_macro("macroiteration", it+1)   # test the NEWLY-PRODUCED iteration against the old
        test_node.add_macro("macroiterationlast", it)
        test_node.add_parent(parent_fit_node)
        test_node.set_category("CONVERGE")
        dag.add_node(test_node)
        test_node_list.append(test_node)

        parent_fit_node=test_node


# Create export stages for extrinsic samples
#
# THE FORK.  Up to here every model shared one grid and one marginalized fit.
# Here each model gets: its own net -> its own terminal CIP (posterior + ln Z)
# -> its own extrinsic ILE chain -> its own extrinsic posterior.  A final
# combine node mixes them with weights p(m) Z_m.  Keeping the per-model
# posteriors on disk is deliberate: a systematics study wants to look at them.
if opts.last_iteration_extrinsic:
    # The terminal CIP reads the grid the loop finished on and writes the
    # per-model posterior at the same index the extrinsic stage reads.
    it = opts.n_iterations

    n_jobs_extrinsic = n_extrinsic_batches
    cat_nodes = []
    loop_final_node = parent_fit_node

    for approx in opts.approx:
        # Per-model net: this model's composites only, pooled flat (replicas of
        # ONE model), with no cross-model marginalization.
        unify_model_node = pipeline.CondorDAGNode(unify_model_job)
        unify_model_node.add_macro("macroiteration", it-1)
        unify_model_node.add_macro("macroapprox", approx)
        unify_model_node.set_retry(opts.general_retries)
        if not(loop_final_node is None):
            unify_model_node.add_parent(loop_final_node)
        dag.add_node(unify_model_node)

        # Per-model terminal fit: posterior AND evidence.
        cipterm_node = pipeline.CondorDAGNode(cip_terminal_job)
        cipterm_node.add_macro("macroiteration", it-1)
        cipterm_node.add_macro("macroiterationnext", it)
        cipterm_node.add_macro("macroapprox", approx)
        cipterm_node.set_category("CIP")
        cipterm_node.add_parent(unify_model_node)
        cipterm_node.set_retry(opts.general_retries)
        dag.add_node(cipterm_node)

        # Create nodes for followup tasks
        cat_node = pipeline.CondorDAGNode(cat_job)

        # Perform final ILE run on all points, saving bounded fair draws.
        # Need to perform number of events CONSISTENT WITH TARGET SAMPLE SIZE
        #    - *not* always same as number of ILE events being analyzed
        #    - *assumes* grid files have sufficiently large numbers of samples to allow this! (as in many other cases)
        for event in np.arange(n_jobs_extrinsic):
            # Add task per ILE operation
            ile_node = pipeline.CondorDAGNode(ileExtr_job)
            ile_node.set_retry(opts.ile_retries)
            ile_node.add_macro("macroevent", event*n_group)
            ile_node.add_macro("macroiteration", it)
            ile_node.add_macro("macroapprox", approx)
            ile_node.add_macro("macrogwsignal", gwsignal_macro(approx))
            ile_node.add_parent(cipterm_node)
            dag.add_node(ile_node)

            # Convert each already-fair-drawn output file.
            for indx in np.arange(n_group):
                convert_node = pipeline.CondorDAGNode(convertExtr_job)
                convert_node.add_macro("macroevent", event*n_group)
                convert_node.add_macro("macroiteration", it)
                convert_node.add_macro("macroindx",indx)
                convert_node.add_macro("macroapprox", approx)
                convert_node.set_retry(opts.ile_retries)  # this can fail too
                convert_node.add_parent(ile_node)

                # Make cat job 
                cat_node.add_parent(convert_node)
                cat_node.set_retry(opts.ile_retries)  # this can fail too
                cat_node.add_macro("macroiteration", it)  # needed to identify log file location
                cat_node.add_macro("macroapprox", approx)

                # Add nodes
                dag.add_node(convert_node)

        dag.add_node(cat_node)
        cat_nodes.append(cat_node)

    # Recombine: mixture over models weighted by p(m) Z_m.
    combine_node = pipeline.CondorDAGNode(combine_job)
    combine_node.add_macro("macroiteration", it)
    combine_node.set_retry(opts.general_retries)
    for node in cat_nodes:
        combine_node.add_parent(node)
    dag.add_node(combine_node)
    parent_fit_node = combine_node

# Create final node for overall plots.  (Note: default setup is designed to enable plots of the last two iterations *at each step* but this seems like overkill)
if plot_args:
        # Cannot run test on first iteration
        plot_node = pipeline.CondorDAGNode(plot_job)
        plot_node.add_macro("macroiteration", it)
        plot_node.add_macro("macroiterationlast", it-1)
        plot_node.add_parent(parent_fit_node)
        plot_node.set_category("PLOT")
        dag.add_node(plot_node)

dag_name="marginalize_intrinsic_parameters_BasicIterationWorkflow"
dag.set_dag_file(dag_name)
dag.write_concrete_dag()

# convergence_test_samples.py returns 1 when convergence is established.
# Match BasicIteration: without a terminal extraction this stops the remaining
# intrinsic loop and reports a successful DAG, rather than leaving a failed
# convergence node that merely blocks all of its children.
if opts.test_args:
    for node in test_node_list:
        dag.add_abort_on(node, exit_code=1, return_value=0)
