text:
#=
    InterpolatedRejectionSampling
    Copyright © 2019 Mark Wells <mwellsa@gmail.com>

    Distributed under terms of the AGPL-3.0 license.
=#

module InterpolatedRejectionSampling

using Interpolations

export irsample!, irsample

include("envelope.jl")

#const maxruns = 1_000_000_000
#const maxruns = 10

#@inline function get_sample!(cache ::Vector{Float64},
#                              pnt ::Union{Missing,Float64},
#                              support ::NTuple{2,Float64})
#    if ismissing(pnt)
#        return support[1] + popfirst!(cache)*support[2]
#    else
#        return pnt
#    end
#end

@inline function get_sample(support ::NTuple{2,Float64}) where D
    return support[1] + rand()*support[2]
end

@inline function get_sample(support ::NTuple{D,NTuple{2,Float64}}) where D
    return Float64[ntuple(i -> get_sample(support[i]), Val(D))...]
end

@inline function get_sample(slice ::AbstractVector{Union{Missing,Float64}},
                            support ::NTuple{D,NTuple{2,Float64}}) where D

    return Float64[ntuple(i -> ismissing(slice[i]) ? get_sample(support[i]) : slice[i], Val(D))...]
end

@inline function rsample(interp ::Interpolations.Extrapolation,
                         support ::NTuple{2,Float64},
                         maxprob ::Float64) where D
    while true
        samp = get_sample(support)
        if rand()*maxprob ≤ interp(samp)
            return samp
        end
    end
    error("unable to draw a sample after $maxruns runs")
end

@inline function rsample(interp ::Interpolations.Extrapolation,
                         support ::NTuple{D,NTuple{2,Float64}},
                         maxprob ::Float64) where D
    while true
        samp = get_sample(support)
        if rand()*maxprob ≤ interp(samp...)
            return samp
        end
    end
    error("unable to draw a sample after $maxruns runs")
end

@inline function rsample(interp ::Interpolations.Extrapolation,
                         slice ::AbstractVector{Union{Missing,Float64}},
                         support ::NTuple{D,NTuple{2,Float64}},
                         maxprob ::Float64) where D

    while true
        samp = get_sample(slice, support)
        if rand()*maxprob ≤ interp(samp...)
            return samp
        end
    end
    error("unable to draw a sample after $maxruns runs")
end

""" irsample(knots, probs, n)
inputs:
    knots ⟶  vector of knots
    probs ⟶  vector of probs
    n ⟶  number of samples to draw
"""
@inline function irsample(knots ::AbstractVector{Float64},
                          probs ::AbstractVector{Float64},
                          n ::Int)
    retval = Vector{Float64}(undef, n)
    interp = LinearInterpolation(knots, probs, extrapolation_bc=Throw())

    support = get_support(knots)
    maxprob = maximum(probs)

    for i in 1:n
        retval[i] = rsample(interp, support, maxprob)
    end
    return retval
end

""" irsample(knots, probs, n)
inputs:
    knots ⟶  tuple of knots (xknots, yknots, ...)
    probs ⟶  D-dimension grid of probability
    n ⟶  number of samples to draw
"""
@inline function irsample(knots ::NTuple{D,AbstractVector{Float64}},
                          probs ::AbstractArray{Float64,D},
                          n ::Int) where D
    retval = Matrix{Float64}(undef, D, n)
    interp = LinearInterpolation(knots, probs, extrapolation_bc=Throw())

    support = get_support(knots)
    maxprob = maximum(probs)
    for i in 1:n
        retval[:,i] = rsample(interp, support, maxprob)
    end
    return retval
end

""" irsample!(slices, knots, probs)
inputs:
    slices ⟶  D by N matrix where D is the dimensionality and N is the number of points, each column is a point (xi,yi,...), missing values will be drawn
    knots ⟶  tuple of knots (xknots, yknots, ...)
    probs ⟶  D-dimension grid of probability
"""
@inline function irsample!(slices ::Matrix{Union{Missing,Float64}},
                           knots ::NTuple{D,AbstractVector{Float64}},
                           probs ::AbstractArray{Float64,D}) where D
    interp = LinearInterpolation(knots, probs, extrapolation_bc=Throw())
    nslices = last(size(slices))
    mask = ismissing.(slices)
    slice_inds = findall([any(mask[:,i]) for i = 1:nslices])

    support = get_support(knots)
    maxprob = maximum(probs)
    for (i,ind) in enumerate(slice_inds)
        slice = selectdim(slices, 2, ind)
        slice .= rsample(interp, vec(slice), support, maxprob)
    end
end

end

#""" irsample!(slices, knots, probs)
#inputs:
#    slices ⟶  D by N matrix where D is the dimensionality and N is the number of points, each column is a point (xi,yi,...), missing values will be drawn
#    knots ⟶  tuple of knots (xknots, yknots, ...)
#    probs ⟶  D-dimension grid of probability
#"""
#@inline function irsample2!(slices ::Matrix{Union{Missing,Float64}},
#                           knots ::NTuple{D,AbstractVector{Float64}},
#                           probs ::AbstractArray{Float64,D};
#                           verbose ::Bool = false) where D
#
#    interp = LinearInterpolation(knots, probs, extrapolation_bc=Throw())
#
#    nslices = last(size(slices))
#    mask = ismissing.(slices)
#    slice_inds = findall([any(mask[:,i]) for i = 1:nslices])
#    nmissing = [count(mask[:,i]) for i = 1:nslices]
#
#    support, pmaxes = get_envelope(slices, knots, probs)
#
#    for run = 0:maxruns
#        #@info "slice_inds" slice_inds
#        _slices = selectdim(slices, 2, slice_inds)
#        _support = view(support, slice_inds)
#        _pmaxes = view(pmaxes, slice_inds)
#        if length(slice_inds) == 0
#            return
#        end
#        #if verbose & (mod(run,10) == 0)
#        if verbose
#            @info """
#            running "irsample!" on run $run with $(length(slice_inds)) to draw
#            """
#        end
#
#        #@info "generating cache " length(slice_inds) nmissing[slice_inds]
#        cache = rand(length(slice_inds) + sum(nmissing[slice_inds]))
#
#        samples = Matrix{Float64}(undef,D,length(slice_inds))
#        for ind in eachindex(_slices)
#            samples[ind] = get_sample!(cache, _slices, ind, _support)
#        end
#
#        variate = cache.*_pmaxes
#
#        inds = findall([variate[ind] ≤ interp(samples[:,ind]...) for ind = 1:length(slice_inds)])
#        # update slices
#        _slices[:,inds] .= samples[:,inds]
#        deleteat!(slice_inds, inds)
#
#        #slice_inds = findall([any(ismissing.(_slices[:,i])) for i = 1:size(_slices)[2]])
#        #_slices = selectdim(slices, 2, slice_inds)
#        #_nslices = length(slice_inds)
#    end
#    nmissing = sum(count.([ismissing.(s) for s in slices]))
#    if nmissing != 0
#        error("failed to draw all samples after maximum number of runs")
#    end
#end

#""" irsample!(slices, knots, probs)
#inputs:
#    slices ⟶  D by N matrix where D is the dimensionality and N is the number of points, each column is a point (xi,yi,...), missing values will be drawn
#    knots ⟶  tuple of knots (xknots, yknots, ...)
#    probs ⟶  D-dimension grid of probability
#"""
#function irsample!(slices ::AbstractMatrix{Union{Missing,Float64}},
#                   knots ::NTuple{D,AbstractVector{Float64}},
#                   probs ::AbstractArray{Float64,D};
#                   verbose ::Bool = false) where D
#
#    interp = LinearInterpolation(knots, probs, extrapolation_bc=Throw())
#    support, pmaxes = get_envelope(slices, knots, probs)
#    pmaxes .+= eps.(pmaxes)
#
#    _, nc = size(slices)
#
#    for run = 1:maxruns
#        mask = BitArray(ismissing.(slices))
#        nmissing = count(mask)
#        if verbose & mod(run,10) == 0
#            @show run
#            @show nmissing
#        end
#
#        if nmissing == 0
#            break
#        end
#
#        slice_inds = findall([any(mask[:,i]) for i = 1:nc])
#        nslices = length(slice_inds)
#
#        cache = rand(nslices + nmissing)
#
#        _slices = selectdim(slices,2,slice_inds)
#
#        _support = view(support,slice_inds)
#        _pmaxes = view(pmaxes,slice_inds)
#        samples = collect(_slices)
#
#        cinds = findall(ismissing.(samples))
#        for (i,cind) in enumerate(cinds)
#            dim = cind[1]
#            ind = cind[2]
#            s_min = _support[ind][dim][1]
#            s_span = _support[ind][dim][2]
#            samples[cind] = s_min + pop!(cache)*s_span
#        end
#
#        variate = cache.*_pmaxes
#
#        inds = findall([variate[i] ≤ interp(samples[:,i]...) for i in eachindex(variate)])
#        # update slices
#        _slices[:,inds] = samples[:,inds]
#    end
#
#    mask = BitArray(ismissing.(slices))
#    nmissing = count(mask)
#    if nmissing != 0
#        error("failed to draw all samples after maximum number of runs")
#    end
#end


text:
"""
    RigidBody <: AbstractModel

Six-degree of freedom rigid-body model with state variables ``x, y, z, \\phi,
\\theta, \\psi, u, v, w, p, q, r``, inputs ``m, I_{xx}, I_{yy}, I_{zz}, I_{xz},
I_{xy}, I_{yz}, F_x, F_y, F_z, M_x, M_y, M_z``, and zero parameters.
"""
struct RigidBody{N,TI,TF} <: AbstractModel
    state_indices::NTuple{N,TI}
    rate_indices::NTuple{N,TI}
    prescribed_values::NTuple{N,TF}
end

# --- Constructors --- #

"""
    RigidBody()

Initialize an object of type [`RigidBody`](@ref)
"""
RigidBody() = RigidBody{0,Int,Float64}((), (), ())

"""
    RigidBody(state_indices, rate_indices [, prescribed_values])

Initialize an object of type [`RigidBody`](@ref).  Use the state variables
corresponding to `state_indices` to set the state rates corresponding to
`rate_indices` to the values specified in `prescribed_values` (which default to zeros
if not specified). Otherwise, allow the state variables and their respective
rates to be defined by their rate equations.
"""
function RigidBody(state_indices::NTuple{N,TI}, rate_indices::NTuple{N,TI}) where {N,TI}

    prescribed_values = ntuple(i->0.0, N)

    return RigidBody(state_indices, rate_indices, prescribed_values)
end

# --- Traits --- #
number_of_states(::Type{<:RigidBody}) = 12
number_of_inputs(::Type{<:RigidBody}) = 13
number_of_parameters(::Type{<:RigidBody}) = 0

inplaceness(::Type{<:RigidBody}) = OutOfPlace()

rate_jacobian_type(::Type{RigidBody{0,TI,TF}}) where {TI,TF} = Identity()
rate_jacobian_type(::Type{RigidBody{N,TI,TF}}) where {N,TI,TF} = Linear()
state_jacobian_type(::Type{<:RigidBody}) = Nonlinear()
input_jacobian_type(::Type{<:RigidBody}) = Nonlinear()
parameter_jacobian_type(::Type{<:RigidBody}) = Nonlinear()
time_gradient_type(::Type{<:RigidBody}) = Zeros()

# --- Methods --- #

function get_residual(model::RigidBody, dstates, states, inputs, parameters, t)
    # prescribed state rates
    dx, dy, dz, dϕ, dθ, dψ, du, dv, dw, dp, dq, dr = dstates
    # extract states
    x, y, z, ϕ, θ, ψ, u, v, w, p, q, r = states
    # extract inputs
    m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz, Fx, Fy, Fz, Mx, My, Mz = inputs
    # rigid body kinematics
    xdot, ydot, zdot, ϕdot, θdot, ψdot = rigid_body_kinematics(x, y, z, ϕ, θ, ψ, u, v, w, p, q, r)
    # rigid body dynamics
    udot, vdot, wdot, pdot, qdot, rdot = rigid_body_dynamics(u, v, w, p, q, r, m,
        Ixx, Iyy, Izz, Ixz, Ixy, Iyz, Fx, Fy, Fz, Mx, My, Mz)
    # calculated state rates
    rates = SVector(xdot, ydot, zdot, ϕdot, θdot, ψdot, udot, vdot, wdot, pdot, qdot, rdot)
    # default residual
    residual = dstates - rates
    # replace elements of the residual vector with new constraints
    for (istate, irate, value) in zip(model.state_indices, model.rate_indices, model.prescribed_values)
        # replace state rate equation for `istate` with the prescribed constraint
        residual = setindex(residual, value - rates[irate], istate)
    end
    # return result
    return residual
end

# --- Performance Overloads --- #

#TODO: Add rate, state, input, and parameter jacobians

# --- Convenience Methods --- #

function set_states!(states, model::RigidBody; x, y, z, ϕ, θ, ψ, u, v, w, p, q, r)

    states[1] = x
    states[2] = y
    states[3] = z
    states[4] = ϕ
    states[5] = θ
    states[6] = ψ
    states[7] = u
    states[8] = v
    states[9] = w
    states[10] = p
    states[11] = q
    states[12] = r

    return states
end

function set_inputs!(inputs, model::RigidBody; m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz,
    Fx, Fy, Fz, Mx, My, Mz)

    inputs[1] = m
    inputs[2] = Ixx
    inputs[3] = Iyy
    inputs[4] = Izz
    inputs[5] = Ixz
    inputs[6] = Ixy
    inputs[7] = Iyz
    inputs[8] = Fx
    inputs[9] = Fy
    inputs[10] = Fz
    inputs[11] = Mx
    inputs[12] = My
    inputs[13] = Mz

    return inputs
end

set_parameters!(p, model::RigidBody) = p

function separate_states(model::RigidBody, states)

    return (x = states[1], y = states[2], z = states[3],
        ϕ = states[4], θ = states[5], ψ = states[6],
        u = states[7], v = states[8], w = states[9],
        p = states[10], q = states[11], r = states[12])
end

function separate_inputs(model::RigidBody, inputs)

    return (m = inputs[1],
        Ixx = inputs[2], Iyy = inputs[3], Izz = inputs[4],
        Ixz = inputs[5], Ixy = inputs[6], Iyz = inputs[7],
        Fx = inputs[8], Fy = inputs[9], Fz = inputs[10],
        Mx = inputs[11], My = inputs[12], Mz = inputs[13])
end

separate_parameters(model::RigidBody, parameters) = ()

# --- Internal Methods --- #

function rigid_body_kinematics(x, y, z, ϕ, θ, ψ, u, v, w, p, q, r)

    Vb = SVector(u, v, w)
    Ωb = SVector(p, q, r)

    sϕ, cϕ = sincos(ϕ)
    sθ, cθ = sincos(θ)
    sψ, cψ = sincos(ψ)

    # linear kinematics
    Rib = @SMatrix [cθ*cψ    cθ*sψ           -sθ;
         sϕ*sθ*cψ - cϕ*sψ sϕ*sθ*sψ + cϕ*cψ sϕ*cθ;
         cϕ*sθ*cψ + sϕ*sψ cϕ*sθ*sψ - sϕ*cψ cϕ*cθ]
    rdot = Rib' * Vb

    # angular kinematics
    ϕdot = p + (q*sϕ + r*cϕ)*sθ/cθ
    θdot = q*cϕ - r*sϕ
    ψdot = (q*sϕ + r*cϕ)/cθ

    return SVector(rdot..., ϕdot, θdot, ψdot)
end

function rigid_body_dynamics(u, v, w, p, q, r, m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz,
    Fx, Fy, Fz, Mx, My, Mz)

    F = SVector(Fx, Fy, Fz)
    M = SVector(Mx, My, Mz)

    Vb = SVector(u, v, w)
    Ωb = SVector(p, q, r)

    # linear dynamics
    vdot = F/m - cross(Ωb, Vb)

    # angular dynamics
    I = @SMatrix [Ixx -Ixy -Ixz; -Iyz Iyy -Iyz; -Ixz -Iyz Izz]
    ωdot = I \ (M - cross(Ωb, I*Ωb))

    return SVector(vdot..., ωdot...)
end

ext: jl


text:
# This file includes an example file to simulate chaotic maps

using ChaoticCommunications
using Plots

# Generate data
cmaps = [Logistic(), Tent(), Logistic(), Henon(), Bernoulli()]
plt = plot(layout=(length(cmaps), 1))
for (i, cmap) in enumerate(cmaps)
    tx = trajectory!(cmap, 100.)
    plot!(tx, marker=(:circle, 1), subplot=i)
end
display(plt)


text:
@testset "397.integer-replacement.jl" begin
    @test integer_replacement(8) == 3
    @test integer_replacement(7) == 4
    for (i, j) in zip([25, 19, 235, 2623, 53], [6, 6, 11, 14, 8])
        @test integer_replacement(i) == j
    end
end


text:
# a modifying sampler modifies one or more distribution samplers or other modifying samplers
# TODO perhaps better name is "MetaSampler"?

abstract ModifyingSampler <: Sampler

paramranges(s::ModifyingSampler) = paramranges(s.subsampler)

getparams(s::ModifyingSampler) = getparams(s.subsampler)

setparams!(s::ModifyingSampler, params) = setparams!(s.subsampler, params)

function sample(s::ModifyingSampler, support, cc::ChoiceContext)
	x, trace = sample(s.subsampler, support, cc)
	x, Dict{Symbol, Any}(:sub=>trace)
end

estimateparams!(s::ModifyingSampler, traces) = estimateparams!(s.subsampler, map(trace->trace[:sub], traces))

amendtrace(s::ModifyingSampler, trace, x) = amendtrace(s.subsampler, trace[:sub], x)

# modifying samplers without a single subsampler, or with more info to display, will need to override
function show(io::IO, s::ModifyingSampler, indentdepth::Int=1)
	print(io, getsamplertypename(s) * " ")
	# note, no newline - keeping the entire output on one line if possible - the root subsampler(s) will handle the newline
	show(io, s.subsampler, indentdepth)
end

include("mixture_sampler.jl")
include("adjust_parameters_to_support_sampler.jl")
include("align_minimum_support_sampler.jl")
include("truncate_to_support_sampler.jl")
include("transform_sampler.jl")
include("constrain_parameters_sampler.jl")
include("conditional_sampler.jl")
include("recursion_depth_sampler.jl")


text:
module AirtableTests

using Airtable
using Airtable.HTTP
using Airtable.JSON3
using ReTest
using Random

const CI_STRING = randstring()

@testset "Airtable.jl" begin
    @testset "Constructors" begin
        bs = AirBase("appphImnhJO8AXmmo")
        @test Airtable.id(bs) == "appphImnhJO8AXmmo"
        @test Airtable.path(bs) == "/v0/appphImnhJO8AXmmo"

        tab = AirTable("Table 1", bs)
        @test Airtable.id(tab) == "Table 1"
        @test Airtable.path(tab) == "/v0/appphImnhJO8AXmmo/Table 1"

        rec = AirRecord("some id", tab, Dict())
        @test Airtable.id(rec) == "some id"
        @test Airtable.path(rec) == "/v0/appphImnhJO8AXmmo/Table 1/some id"

    end

    @testset "Interface" begin
        tab = AirTable("Table 1", AirBase("appphImnhJO8AXmmo"))

        @test length(Airtable.query(tab; filterByFormula="{Keep}")) == 3

        toadd = [
            (; Name   = "TEST1",
               Notes  = "Some note",
               Status = "Todo",
               CI     = CI_STRING),
            (; Name   = "TEST2",
               Notes  = "Other note",
               Status = "Done",
               CI     = CI_STRING)
        ]

        resp = Airtable.post!(tab, toadd)
        @test resp isa Vector{AirRecord}
        @test length(resp) == 2

        for rec in resp
            sleep(1) # avoid over-taxing api

            @test rec[:Status] != "In progress"
            @test Airtable.patch!(rec, (; Status="In progress", CI=CI_STRING)).id == Airtable.id(rec)
            @test Airtable.get(rec)[:Status] == "In progress"

            @test_throws HTTP.ExceptionRequest.StatusError Airtable.patch!(rec, (; Status="Not valid"))
            Airtable.patch!(rec, (; Status = rec[:Status]))

            @test keys(rec) == (:Name, :Notes, :Status, :CI)
            @test Airtable.delete!(rec).deleted
            @test_throws HTTP.ExceptionRequest.StatusError Airtable.get(rec)
        end

        JSON3.write("add_records.json", (; records = [(; fields = rec) for rec in toadd]))
        resp = Airtable.post!(tab, open("add_records.json"))

        Airtable.patch!(tab, resp, [(; Status="In progress", CI=CI_STRING) for _ in 1:length(resp)])
        @test all([Airtable.get(rec)[:Status] == "In progress" for rec in resp])

    end
    # Cleanup
    dontkeep = Airtable.query(AirTable("Table 1", AirBase("appphImnhJO8AXmmo")); filterByFormula="AND({CI}='$CI_STRING', NOT({Keep}))")
    if !isempty(dontkeep)
        sleep(1)
        for rec in dontkeep
            Airtable.delete!(rec)
        end
    end
end

end


text:
module NCU

using QXZoo.Circuit
using QXZoo.GateOps
using QXZoo.GateMap
using QXZoo.DefaultGates

# =========================================================================== #
#
# =========================================================================== #

"""
    get_intermed_gate(U::GateOps.GateSymbol)

Generate sqrt U and adjoint sqrt U, given a previously cached gate U.
"""
function get_intermed_gate(U::GateOps.GateSymbol)
    su = GateOps.GateSymbol(U.label, U.rt_depth+1)
    asu = GateOps.GateSymbol(U.label, U.rt_depth+1, true)

    if ~haskey(GateMap.gates, su)

        SU = ()->sqrt(GateMap.gates[U]())
        ASU = ()->adjoint(SU())

        # Cache the gates
        GateMap.gates[su] = SU
        GateMap.gates[asu] = ASU
    end

    return su, asu
end

# =========================================================================== #
#
# =========================================================================== #

"""
    apply_ncu!(cct::Circuit.Circ, q_ctrl::Vector, q_aux::Vector, q_tgt, U::GateOps.GateSymbol)

Apply an n-qubit controlled gate operation on the given target qubit.
Ensure the 2-qubit gate corresponding with symbol c_U is registered with GateMap.gates before use.
Appends the GateCall operations to circuit

# Arguments
- `cct::Circuit.Circ`
- `ctrls::Vector`:
- `aux::Vector`:
- `tgt::Int`:
- `U::GateOps.GateSymbol`: symbol for 2-qubit gate to extend (c_x, c_z, etc)
"""
function apply_ncu!(cct::Circuit.Circ, q_ctrl::Vector, q_aux::Vector, q_tgt, c_U::GateOps.GateSymbol)

    # Check global Circuit cache for gate
    if ~haskey(GateMap.gates, c_U )
        error("Gate $(U) does not exist")
    end

    su,asu = get_intermed_gate(c_U)

    if length(q_ctrl) == 2
        cct << DefaultGates.c_u(su, q_tgt, q_ctrl[2])
        cct << DefaultGates.c_x(q_ctrl[2], q_ctrl[1])
        cct << DefaultGates.c_u(asu, q_tgt, q_ctrl[2])
        cct << DefaultGates.c_x(q_ctrl[2], q_ctrl[1])
        cct << DefaultGates.c_u(su, q_tgt, q_ctrl[1])

    elseif length(q_ctrl) == 3 # Optimisation for 3 control lines

        ssu,assu = get_intermed_gate(su)

        cct << DefaultGates.c_u(ssu, q_tgt, q_ctrl[1])
        cct << DefaultGates.c_x(q_ctrl[2], q_ctrl[1])
        cct << DefaultGates.c_u(assu, q_tgt, q_ctrl[2])
        cct << DefaultGates.c_x(q_ctrl[2], q_ctrl[1])

        cct << DefaultGates.c_u(ssu, q_tgt, q_ctrl[2])
        cct << DefaultGates.c_x(q_ctrl[3], q_ctrl[2])
        cct << DefaultGates.c_u(assu, q_tgt, q_ctrl[3])
        cct << DefaultGates.c_x(q_ctrl[3], q_ctrl[1])

        cct << DefaultGates.c_u(ssu, q_tgt, q_ctrl[3])
        cct << DefaultGates.c_x(q_ctrl[3], q_ctrl[2])
        cct << DefaultGates.c_u(assu, q_tgt, q_ctrl[3])
        cct << DefaultGates.c_x(q_ctrl[3], q_ctrl[1])
        cct << DefaultGates.c_u(ssu, q_tgt, q_ctrl[3])

    # Optimisation for Multi-controlled x with additional qubit availability
    elseif (length(q_ctrl)>=5) && (length(q_aux)>=length(q_ctrl)-2) && (c_U == DefaultGates.GateSymbols.c_x)

        apply_ncu!(cct, [q_ctrl[end], q_aux[ 1 + length(q_ctrl)-3 ]], Int[], q_tgt, c_U)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(cct, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], c_U)
        end

        apply_ncu!(cct, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], c_U)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(cct, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], c_U)
        end

        apply_ncu!(cct, [q_ctrl[end], q_aux[1 + length(q_ctrl) - 3]], Int[], q_tgt, c_U)
        for i in reverse(2:length(q_ctrl)-2)
            apply_ncu!(cct, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], c_U)
        end

        apply_ncu!(cct, [q_ctrl[1], q_ctrl[2]], Int[], q_aux[1], c_U)
        for i in 2:length(q_ctrl)-2
            apply_ncu!(cct, [q_ctrl[1 + i], q_aux[1 + (i-2)]], Int[], q_aux[1 + (i-1)], c_U)
        end

    elseif (length(q_ctrl)>=5) && (length(q_aux)>=length(q_ctrl)-2) && (c_U == DefaultGates.GateSymbols.c_z)

        cct << DefaultGates.h(q_tgt)
        apply_ncu!(cct, q_ctrl, q_aux, q_tgt, DefaultGates.GateSymbols.c_x)
        cct << DefaultGates.h(q_tgt)

    else
        # By cat'ing the unused qubits as part of this recursive decomposition steps, large
        # controlled circuits eventually converge to the auxiliary assisted routines,
        # assuming controlled-x operations.
        cct << DefaultGates.c_u(su, q_tgt, q_ctrl[end])
        apply_ncu!(cct, q_ctrl[1:end-1], vcat(q_aux, q_tgt), q_ctrl[end], DefaultGates.GateSymbols.c_x)
        cct << DefaultGates.c_u(asu, q_tgt, q_ctrl[end])
        apply_ncu!(cct, q_ctrl[1:end-1], vcat(q_aux, q_tgt), q_ctrl[end], DefaultGates.GateSymbols.c_x)
        apply_ncu!(cct, q_ctrl[1:end-1], vcat(q_aux, q_ctrl[end]), q_tgt, su)
    end
    ;
end

end


text:
header = "CC,model_flag,SP_Allowance_1,SP_Allowance_2,SP_Allowance_3,SP_Fee_1,SP_Fee_2,Sp_Fee_3,SP_Fee_4,Num_Carers_1,Num_Carers_2,Num_Carers_3,Num_Carers_4,Number of children/young people in placements (at 31 March),Number of disabled children/young people (at 31 March)^,White (at 31 March),Mixed (at 31 March),Asian/Asian British (at 31 March),Black/Black British (at 31 March),Other ethnic groups (at 31 March),Not known (at 31 March),Number of CYP considered to be at risk of CSE,Number of children/young people who went ‘missing’ while living with foster carers,Number of instances children/young people went ‘missing’ while living with foster carers,Less than 24 hours,One to six days,One week to 28 days,Longer than 28 days,Bullying,Involvement in offending behaviour,Substance misuse,Contact with family or friends,Alleged abuse within the foster home,Because they are an asylum seeker,Suspected sexual exploitation,Any other known reason,Reason unknown,Number of return home interviews conducted,Number of missing children/young people with no return home interview,Number of children experiencing unplanned endings,Number of children experiencing placement moves within 24 hours of unplanned ending notification,Transferred to LA fostering placement,Transferred to IFA fostering placement,Transferred to other type of looked after placement (e.g. children's home, adoption),No longer looked after,Staying put,Not staying put,Reason unknown,Number of young people in foster care who became 18,Number of young people in foster care who became 18 and remained living with their former foster carer,Groups assessed,Siblings assessed,Groups placed to plan,Siblings placed to plan,Groups NOT placed to plan,Siblings NOT placed to plan,Total approved carers ,White,Mixed,Asian/Asian British,Black/Black British,Other ethnic groups,Not known,Number of individual foster carers who had completed the TSD Standards,Number of individual foster carers who were awaiting sign-off on the TSD Standards,Number of individual foster carers who were in progress completing the TSD Standards,Number of individual foster carers who had not yet commenced the TSD Standards,Number of individual foster carers with TSD Standards training status unknown,Total approved households,Permanent foster care,Non-permanent foster care,Emergency,Family and Friends,Parent and child,Multi-dimensional treatment,Remand,Other,Fostering to adopt,Short breaks,Shared care,Not known,Number of households with exemptions to the usual fostering limit in place at any time during the year,Number of households approved for one CYP,Number of households approved for two CYP,Number of households approved for three CYP,Number of households approved for four or more CYP,Number of approved foster places,Number of filled places,Number of vacant places,Number of places not available,Number of initial enquiries from new prospective fostering households,Total number of applications received from prospective fostering households between 1 April and 31 March,Number of applications approved,Number of applications withdrawn by applicants,Number of applications withdrawn by services,Number of applications rejected after being presented to panel,Number of applications which remained in progress on 31 March,Total number of households approved in year,"Of which, number of family and friends households approved in year ~",Deregistered to facilitate a transfer to a local authority,Deregistered to facilitate a transfer to an independent fostering agency,Deregistration initiated by the fostering service,Deregistration initiated by the foster carer,Deregistration for unknown reason,"Total households deregistered between 1 April and 31 March""",Number of applications for a review to the IRM,Number of which resulted in a recommendation being made to review the decision,Number of allegations against foster carers by fostered children,Number of allegations against foster carers by other sources,Total number of allegations made,Number of children made allegations against foster carers,Number of foster carers subject to allegations,Physical abuse,Sexual abuse,Neglect,Emotional abuse,"Concerns resolved, no further action",Continued monitoring for agreed period,"Concern remained, refer to fostering panel for review",Less than 21 working days,22-30 working days (4-6 weeks),31-50 working days (6-10 weeks),More than 50 working days (more than 10 weeks),Number of referrals to the Disclosure and Barring Service (DBS),Number of referrals to the local authority child protection team,Number of Section 47 enquiries,Number of recorded incidents of physical restraint by foster carers,Number of CYP subject to restraint,Number of households that restrained CYP,CYP (complaints made),CYP (complaints not upheld),Foster carers (complaints made),Foster carers (complaints not upheld),CYP's social workers (complaints made),CYP's social workers (complaints not upheld),Other professionals (complaints made),Other professionals (complaints not upheld),CYP's parents (complaints made),CYP's parents (complaints not upheld),Others (complaints made),Others (complaints not upheld),Total complaints (complaints made),Total complaints (complaints not upheld)"

d = split( header, ",")