text:
@testset "utils" begin
    @test contains(basename(Playground.julia_url(Playground.NIGHTLY)), "julia-latest")
    @test contains(basename(Playground.julia_url(v"0.4+")), "julia-0.4-latest")
    @test contains(basename(Playground.julia_url(v"0.4")), "julia-0.4.0")
    @test_throws ArgumentError Playground.julia_url(Base.nextpatch(Playground.NIGHTLY))

    @test contains(basename(Playground.julia_url("nightly")), "julia-latest")
    @test ismatch(r"julia-.*-latest", basename(Playground.julia_url("release")))
    @test contains(basename(Playground.julia_url("0.4")), "julia-0.4.0")

    # Testing coverage
    versions = (v"0.5", v"0.6", v"0.7-")
    platforms = (:Darwin, :Linux, :Windows)
    archs = (32, 64)
    @testset "julia $v on $p-$a" for v in versions, p in platforms, a in archs
        if p === :Darwin && a == 32
            @test_throws ErrorException Playground.julia_url(v, p, a)
        else
            @test !isempty(Playground.julia_url(v, p, a))
        end
    end
end


text:
using Pkg
Pkg.activate("c:/git/JLBoost")
@time using DataFrames
@time using JDF
@time using JLBoost, LossFunctions

###############################################################################
# fitting fm
###############################################################################
if !isdir("c:/data/GiveMeSomeCredit/cs-training.jdf")
	using CSV
	a = CSV.read("c:/data/GiveMeSomeCredit/cs-training.csv", missingstrings=["NA"])
	rename!(a, Symbol("")=>:column1)
	create_missing!(a, :MonthlyIncome)
	create_missing!(a, :NumberOfDependents)
	type_compress!(a, compress_float=true)
	savejdf("c:/data/GiveMeSomeCredit/cs-training.jdf", a)
end

a = loadjdf("c:/data/GiveMeSomeCredit/cs-training.jdf")

using Calculus

softmax(w) = 1/(1 + exp(-w))
logloss(y) = w -> logloss(y, w)
logloss(y, w) = -(y*log( 1/(1 + exp(-w))) + (1-y)*log(1- 1/(1 + exp(-w))))

t = a.SeriousDlqin2yrs
ws = fill(0.0,nrow(a))

using CuArrays
ct = cu(t)
cws = cu(ws)

logloss.(ct, cws)

ct .+ cws

using BenchmarkTools
@benchmark Calculus.derivative.(logloss.($t), $ws)

using LossFunctions: deriv
loss = LogitLogLoss()
@benchmark deriv(loss, $t, $ws)

using ForwardDiff
@benchmark ForwardDiff.derivative.(logloss.($t), $ws)

using ForwardDiff

d = transpose(hcat(t, ws))
@benchmark ForwardDiff.gradient.(x->logloss(x...), eachcol(d))

using Zygote: @adjoint, gradient

@adjoint logloss(y, w) = logloss(y, w), Δ -> (Δ*-(log(1/(1+exp(-w))) - (log(1 - 1/(1+exp(-w))))), Δ*( -1/(1+exp(w)) - y + 1))

@time gradient(logloss, 1, 0.5)

f(t, ws) = gradient(logloss, t, ws)[2]

@benchmark f.(t, ws)

# Winner is LossFunctions

using CuArrays
ct = cu(t)
cws = cu(ws)

g(ct, cws) = deriv(LogitLogloss(), ct, cws)

@benchmark g.(loss, ct, cws)

#####################################################################

using Zygote: @adjoint, gradient
gg(ct, cws) = gradient(logloss, ct, cws)[2]
@time gg.(ct, cws)

logloss.(ct, cws)

|
d1(y,w) = @. -1 / (1 + exp(w)) - y + 1

using BenchmarkTools
CuArrays.allowscalar(false)
@benchmark d1(ct, cws)

ext: jl


text:
@inline Base.isempty(::EmptyCell) = true
@inline Base.isempty(::AbstractCell) = false
@inline iserror(c::Cell) = c.datatype == "e"
@inline iserror(::AbstractCell) = false
@inline row_number(c::EmptyCell) = row_number(c.ref)
@inline column_number(c::EmptyCell) = column_number(c.ref)
@inline row_number(c::Cell) = row_number(c.ref)
@inline column_number(c::Cell) = column_number(c.ref)
@inline relative_cell_position(c::Cell, rng::CellRange) = relative_cell_position(c.ref, rng)
@inline relative_cell_position(c::EmptyCell, rng::CellRange) = relative_cell_position(c.ref, rng)
@inline relative_column_position(c::Cell, rng::ColumnRange) = relative_column_position(c.ref, rng)
@inline relative_column_position(c::EmptyCell, rng::ColumnRange) = relative_column_position(c.ref, rng)

Base.:(==)(c1::Cell, c2::Cell) = c1.ref == c2.ref && c1.datatype == c2.datatype && c1.style == c2.style && c1.value == c2.value && c1.formula == c2.formula
Base.hash(c::Cell) = hash(c.ref) + hash(c.datatype) + hash(c.style) + hash(c.value) + hash(c.formula)

Base.:(==)(c1::EmptyCell, c2::EmptyCell) = c1.ref == c2.ref
Base.hash(c::EmptyCell) = hash(c.ref) + 10

function Cell(c::EzXML.Node)
    # c (Cell) element is defined at section 18.3.1.4
    # t (Cell Data Type) is an enumeration representing the cell's data type. The possible values for this attribute are defined by the ST_CellType simple type (§18.18.11).
    # s (Style Index) is the index of this cell's style. Style records are stored in the Styles Part.

    @assert EzXML.nodename(c) == "c" "`Cell` Expects a `c` (cell) XML node."

    ref = CellRef(c["r"])

    # type
    if haskey(c, "t")
        t = c["t"]
    else
        t = ""
    end

    # style
    if haskey(c, "s")
        s = c["s"]
    else
        s = ""
    end

    # iterate v and f elements
    local v::String = ""
    local f::String = ""
    local found_v::Bool = false
    local found_f::Bool = false

    for c_child_element in EzXML.eachelement(c)
        if EzXML.nodename(c_child_element) == "v"

            # we should have only one v element
            if found_v
                error("Unsupported: cell $(ref) has more than 1 `v` elements.")
            else
                found_v = true
            end

            v = EzXML.nodecontent(c_child_element)
        elseif EzXML.nodename(c_child_element) == "f"

            # we should have only one f element
            if found_f
                error("Unsupported: cell $(ref) has more than 1 `f` elements.")
            else
                found_f = true
            end

            f = EzXML.nodecontent(c_child_element)
        end
    end

    return Cell(ref, t, s, v, f)
end

@inline getdata(ws::Worksheet, empty::EmptyCell) = missing

const RGX_INTEGER = r"^\-?[0-9]+$"

"""
    getdata(ws::Worksheet, cell::Cell) :: CellValue

Returns a Julia representation of a given cell value.
The result data type is chosen based on the value of the cell as well as its style.

For example, date is stored as integers inside the spreadsheet, and the style is the
information that is taken into account to chose `Date` as the result type.

For numbers, if the style implies that the number is visualized with decimals,
the method will return a float, even if the underlying number is stored
as an integer inside the spreadsheet XML.

If `cell` has empty value or empty `String`, this function will return `missing`.
"""
function getdata(ws::Worksheet, cell::Cell) :: CellValueType

    if iserror(cell)
        return missing
    end

    if cell.datatype == "inlineStr"
        error("Excel data type inlineStr is not supported.")
    end

    if cell.datatype == "s"

        if isempty(cell.value)
            return missing
        end

        # use sst
        str = sst_unformatted_string(ws, cell.value)

        if isempty(str)
            return missing
        else
            return str
        end

    elseif (isempty(cell.datatype) || cell.datatype == "n")

        if isempty(cell.value)
            return missing
        end

        if !isempty(cell.style) && styles_is_datetime(ws, cell.style)
            # datetime
            return _celldata_datetime(cell.value, isdate1904(ws))

        elseif !isempty(cell.style) && styles_is_float(ws, cell.style)

            # float
            return parse(Float64, cell.value)

        else
            # fallback to unformatted number
            if occursin(RGX_INTEGER, cell.value)  # if contains only numbers
                v_num = parse(Int, cell.value)
            else
                v_num = parse(Float64, cell.value)
            end

            return v_num
        end
    elseif cell.datatype == "b"
        if cell.value == "0"
            return false
        elseif cell.value == "1"
            return true
        else
            error("Unknown boolean value: $(cell.value).")
        end
    elseif cell.datatype == "str"
        # plain string
        if isempty(cell.value)
            return missing
        else
            return cell.value
        end
    end

    error("Couldn't parse data for $cell.")
end

function _celldata_datetime(v::AbstractString, _is_date_1904::Bool) :: Union{Dates.DateTime, Dates.Date, Dates.Time}

    # does not allow empty string
    @assert !isempty(v) "Cannot convert an empty string into a datetime value."

    if occursin(".", v) || v == "0"
        time_value = parse(Float64, v)
        @assert time_value >= 0

        if time_value <= 1
            # Time
            return excel_value_to_time(time_value)
        else
            # DateTime
            return excel_value_to_datetime(time_value, _is_date_1904)
        end
    else
        # Date
        return excel_value_to_date(parse(Int, v), _is_date_1904)
    end
end

# Converts Excel number to Time.
# `x` must be between 0 and 1.
# To represent Time, Excel uses the decimal part
# of a floating point number. `1` equals one day.
function excel_value_to_time(x::Float64) :: Dates.Time
    @assert x >= 0 && x <= 1
    return Dates.Time(Dates.Nanosecond(round(Int, x * 86400) * 1E9 ))
end

time_to_excel_value(x::Dates.Time) :: Float64 = Dates.value(x) / ( 86400 * 1E9 )

# Converts Excel number to Date. See also XLSX.isdate1904.
function excel_value_to_date(x::Int, _is_date_1904::Bool) :: Dates.Date
    if _is_date_1904
        return Dates.Date(Dates.rata2datetime(x + 695056))
    else
        return Dates.Date(Dates.rata2datetime(x + 693594))
    end
end

function date_to_excel_value(date::Dates.Date, _is_date_1904::Bool) :: Int
    if _is_date_1904
        return Dates.datetime2rata(date) - 695056
    else
        return Dates.datetime2rata(date) - 693594
    end
end

# Converts Excel number to DateTime.
# The decimal part represents the Time (see `_time` function).
# The integer part represents the Date.
# See also XLSX.isdate1904.
function excel_value_to_datetime(x::Float64, _is_date_1904::Bool) :: Dates.DateTime
    @assert x >= 0

    local dt::Dates.Date
    local hr::Dates.Time

    dt_part = trunc(Int, x)
    hr_part = x - dt_part

    dt = excel_value_to_date(dt_part, _is_date_1904)
    hr = excel_value_to_time(hr_part)

    return dt + hr
end

function datetime_to_excel_value(dt::Dates.DateTime, _is_date_1904::Bool) :: Float64
    return date_to_excel_value(Dates.Date(dt), _is_date_1904) + time_to_excel_value(Dates.Time(dt))
end


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using Multivectors
using LinearAlgebra
using StaticArrays

module KPKG3
  using Multivectors
  @generate_basis("0+++")
end
using .KPKG3

@testset "KVectors" begin
  e₁, e₂, e₃, e₄ = alle( KPKG3, 4)[1:4]
  a = e₁(1.0); b = e₂(2.0); c = e₃(3.0); d = e₄(4.0)
  @test typeof(a+a) <: Blade
  @test (typeof(a+b) <: Blade) == false
  @test typeof(a+b) <: KVector
  B = a+b
  B2 = b+c
  @test a+b+c+b+a == e₁(2.0) + e₂(4.0) + e₃(3.0) == 2.0a+2.0b+c
  @test a+ (-b) == a - b
  @test B+B2 == B2+B == a+2.0b+c
  @test (a+b)+(c+d) == a+b+c+d
  @test 5.0*(a+b+c+d) == 5.0*b+5.0*c+5.0*d+5.0*a
  @test grade(a+b+c) == 1
  @test grade(KVector(a*b)) == 2
  @test dual(dual(B)) == B
  @test zero(KVector(a))+b == KVector(b) == b + zero(KVector(a))
  @test zero(B)*B == zero(B) == B*zero(B) == B*zero(KVector{Float64, grade(B), 1})
  @test 3.0*B == B+B+B == B*3.0
  @test iszero(B∧B)
  @test 2.0∧B∧B2 == B∧b∧2.0 - c∧B∧2.0
  @test grade(B∧B2) == 2
  @test grade(B,1) == B
  @test reverse(B) == B
  @test reverse(B∧B2) == -B∧B2
  @test B-B == -a - -B - b

  @test dual(prune(B-B)) == prune(B-B)
  @test iszero(KVector(0.0e₂))
  @test iszero(1.0e₁-1.0e₁)
  @test dual(B) == !B

#!me passes with Blades v0.1.1+  @test normalize(KVector(-2.2(e₂∧e₃))) == KVector(normalize(-2.2(e₂∧e₃))) == Multivectors.normalize_safe(KVector(-2.2(e₂∧e₃)))
end

module KG3
  using Multivectors
  @generate_basis("+++",false,true,true)
end
using .KG3

@testset "More KVectors" begin
  e₁, e₂, e₃ = alle(KG3, 3)[1:3]
  𝐼 = alle(KG3,3)[end]

  a = sort_basis(1.0e₁ + 3.0e₃)

  @test eltype(a) <: Blade
  @test size(a) == (2,)
  @test a == KVector([1.0,0.0,3.0], KG3) |> Multivectors.prune
  @test sort_basis(1.0e₂) == 1.0e₂
  @test convert(KVector, 1.0e₂) == KVector(1.0e₂)
  @test first(a) == a[1]
  @test firstindex(a) == 1
  @test a[end] == a[2]
  @test isnull(a) == false
  @test length(a) == 2
  @test isempty(a) == false
  @test [i for i in a] == map(i->i, a) |> collect == (i for i in a) |> collect
  B = -1.0(e₁∧e₂) + 2.0(e₁∧e₃)
  x = 0.0
  for i in a
    x = x+scalar(i)
  end
  @test x == mapreduce(scalar, +, a)

  @test promote_rule(typeof(KVector(1e₁)), typeof(1e₁)) == KVector

  @test conj(a) == a
  @test conj(B) == -B
  @test KVector(a) == a
  @test KVector([1,2,3], 𝐼) == 1e₁+2e₂+3e₃
  @test pseudoscalar(a) == pseudoscalar(a[1])
  @test grade(⟂(a)∧a) == grade(pseudoscalar(a))
  @test a/2.0 == a*0.5
  @test !(a == a∧e₂(1.0))

  @test coords(a) == scalar.(sort_basis(a+0.0e₂))
  @test coords(a[1]) == [scalar(a[1]), 0.0, 0.0]
  @test Multivectors.prune(KVector(coords(a) .* basis_1blades(a))) == a
  @test norm(basis_1vector(a)) == sqrt(3.0)
  @test Multivectors.norm_sqr(a) == mapreduce(aᵢ->aᵢ*aᵢ, +, a)
  @test norm(Multivectors.normalize_safe(a)) == norm(normalize(a))

  @test grade(⋆(a)) == grade(pseudoscalar(a))-grade(a)

  # warning these relys on LinearAlebra.:⋅ (dot product)
  # which only works on KVectors of grade p with matching sorted basis p-vectors present
  # for a robust inner product you should be using Multivectors.jl
  @test det(1.0e₁+1.0e₂, 1.0e₁+1.0e₂) == 1.0
  @test det(1.0e₁+1.0e₂, 1.0e₁+1.0e₂+1.0e₃) == 0.0
  @test cos(KVector([1.0,0.0,0.0], 𝐼), KVector([1.0,1.0,0.0], 𝐼)) ≈ cos(π/4)

  fᵢ = [rand()*e₁ + rand()*e₂ + rand()*e₃ for i in 1:3]
  fⁱ = reciprocal_frame(fᵢ)
  @test [a⋅b for a in fᵢ for b in fⁱ] |> sum ≈ 3.0

  # this is a neat property of reciprocal frames fᵢfⁱ == n via geometric product,  an integer!
  @test grade(mapreduce(*, +, fᵢ, fⁱ), 0) ≈ 3.0
  @test norm(grade(mapreduce(*, +, fᵢ, fⁱ), 0)) ≈ 3.0
  # not neccessarily parrallel
  @test !(norm.(map(∧, fᵢ, fⁱ)) |> sum ≈ 0.0)
  # only if source frame was orthogonal
  eᵢ = 1.0 .* basis_1blades(e₁)
  eⁱ = reciprocal_frame(eᵢ)
  @test norm.(map(∧,eᵢ, eⁱ)) |> sum ≈ 0.0

end

@testset "Null KVectors" begin
  e₁, e₂, e₃ = alle(KG3, 3)[1:3]
  n = KVector{Float64,1,0}()
  nn = KVector{Float64,1,0}()
  a = KVector{Float64,1,1}([1.0e₁])

  @test iszero(n)
  @test isnull(n+n)
  @test n+a == a == a+n
  @test isnull(n*n)
  @test iszero(dual(n))
  @test iszero(⟂(n))
end


text:
abstract type AbstractBorderedLinearSolver <: AbstractLinearSolver end

# the following stuctures, say `struct BDLS;...;end` rely on the hypotheses:
# - the constructor must provide BDLS() and BDLS(::AbstractLinearSolver)
# - the method (ls::BDLS)(J, dR, dzu, dzp, R, n, ξu, ξp; shift = nothing) must be provided

# call for using BorderedArray input, specific to Arclength Continuation
(lbs::AbstractBorderedLinearSolver)(J, dR,
							dz::BorderedArray, R,
							n::T, theta::T;
							shift::Ts = nothing) where {T, Ts} = (lbs)(J, dR,
																	dz.u, dz.p,
																	R, n,
																	theta / length(dz.u), one(T) - theta;
																	shift = shift)

####################################################################################################
"""
$(TYPEDEF)

This struct is used to provide the bordered linear solver based on the Bordering Method.

$(TYPEDFIELDS)
"""
@with_kw struct BorderingBLS{S <: Union{AbstractLinearSolver, Nothing}, Ttol} <: AbstractBorderedLinearSolver
	"Linear solver used for the Bordering method."
	solver::S = nothing

	"Tolerance for checking precision"
	tol::Ttol = 1e-12

	"Check precision of the linear solve?"
	checkPrecision::Bool = false
end

# dummy constructor to simplify user passing options to continuation
BorderingBLS(ls::AbstractLinearSolver) = BorderingBLS(solver = ls)

# solve in dX, dl
# ┌                           ┐┌  ┐   ┌   ┐
# │ (shift⋅I + J)     dR      ││dX│ = │ R │
# │  ξu * dz.u'    ξp * dz.p  ││dl│   │ n │
# └                           ┘└  ┘   └   ┘
function (lbs::BorderingBLS)(  J, dR,
								dzu, dzp::T, R, n::T,
								ξu::Tξ = 1, ξp::Tξ = 1; shift::Ts = nothing)  where {T, Tξ, Ts}
	# the following parameters are used for the pseudo arc length continuation
	# ξu = θ / length(dz.u)
	# ξp = 1 - θ

	# we make this branching to avoid applying a zero shift
	if isnothing(shift)
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR)
	else
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR; a₀ = shift)
	end

	dl = (n - dot(dzu, x1) * ξu) / (dzp * ξp - dot(dzu, x2) * ξu)

	# dX = x1 .- dl .* x2
	axpy!(-dl, x2, x1)

	# we check the precision of the solution from the bordering algorithm
	# mainly for debugging purposes
	if lbs.checkPrecision
		# at this point, x2 is not used anymore, we can use it for computing the residual
		# hence x2 = J*x1 + dl*dR - R
		x2 = apply(J, x1)
		axpy!(dl, dR, x2)
		axpy!(-1, R, x2)
		if norm(x2) > lbs.tol || abs(n - ξp * dzp * dl -ξu * dot(dzu, x1)) > lbs.tol
			@warn "BorderingBLS did not achieve tolerance"
		end
	end
	return x1, dl, true, (it1, it2)
end
####################################################################################################
# this interface should work for Sparse Matrices as well as for Matrices
"""
$(TYPEDEF)
This struct is used to  provide the bordered linear solver based on inverting the full matrix.
$(TYPEDFIELDS)
"""
struct MatrixBLS{S <: Union{AbstractLinearSolver, Nothing}} <: AbstractBorderedLinearSolver
	"Linear solver used to invert the full matrix."
	solver::S
end

# dummy constructor to simplify user passing options to continuation
MatrixBLS() = MatrixBLS(nothing)

# case of a scalar additional linear equation
function (lbs::MatrixBLS)(J, dR,
						dzu, dzp::T, R::AbstractVecOrMat, n::T,
						ξu::T = T(1), ξp::T = T(1);
						shift::Ts = nothing)  where {T <: Number, S, Ts}

	if isnothing(shift)
		A = J
	else
		A = J + shift * I
	end

	A = hcat(A, dR)
	A = vcat(A, hcat((dzu .* ξu)', dzp * ξp))

	# solve the equations and return the result
	rhs = vcat(R, n)
	res = A \ rhs
	return res[1:end-1], res[end], true, 1
end

####################################################################################################
# composite type to save the bordered linear system with expression
# ┌         ┐
# │  J    a │
# │  b'   c │
# └         ┘
# It then solved using Matrix Free algorithm applied to the full operator and not just J as for MatrixFreeBLS
#
struct MatrixFreeBLSmap{Tj, Ta, Tb, Tc, Ts}
	J::Tj
	a::Ta
	b::Tb
	c::Tc
	shift::Ts
end

function (lbmap::MatrixFreeBLSmap{Tj, Ta, Tb, Tc, Ts})(x::BorderedArray{Ta, Tc}) where {Tj, Ta, Tb, Tc <: Number, Ts}
	# This implements the case where Tc is a number, ie there is one scalar constraint in the
	# bordered linear system
	out = similar(x)
	copyto!(out.u, apply(lbmap.J, x.u))
	axpy!(x.p, lbmap.a, out.u)
	if isnothing(lbmap.shift) == false
		axpy!(lbmap.shift, x.u, out.u)
	end
	out.p = dot(lbmap.b, x.u)  + lbmap.c  * x.p
	return out
end

function (lbmap::MatrixFreeBLSmap)(x::AbstractArray)
	# This implements the case where Tc is a number, ie there is one scalar constraint in the
	# bordered linear system
	out = similar(x)
	xu = @view x[1:end-1]
	xp = x[end]
	# copyto!(out.u, apply(lbmap.J, x.u))
	if isnothing(lbmap.shift)
		out[1:end-1] .= apply(lbmap.J, xu) .+ xp .* lbmap.a
	else
		out[1:end-1] .= apply(lbmap.J, xu) .+ xp .* lbmap.a .+ xu .* lbmap.shift
	end
	out[end] = dot(lbmap.b, xu)  + lbmap.c  * xp
	return out
end

"""
$(TYPEDEF)

This struct is used to  provide the bordered linear solver based a matrix free operator for the full system in `(x,p)`.

$(TYPEDFIELDS)
"""
struct MatrixFreeBLS{S <: Union{AbstractLinearSolver, Nothing}} <: AbstractBorderedLinearSolver
	"Linear solver used to solve the extended linear system"
	solver::S
	"What is the structure used to hold `(x,p)`. If `true`, this is achieved using `BorderedArray`. If `false`, a `Vector` is used."
	useBorderedArray::Bool
end

# dummy constructor to simplify user passing options to continuation
MatrixFreeBLS(useBorderedArray::Bool = true) = MatrixFreeBLS(nothing, useBorderedArray)
MatrixFreeBLS(::Nothing) = MatrixFreeBLS()
MatrixFreeBLS(S::AbstractLinearSolver) = MatrixFreeBLS(S, ~(S isa GMRESIterativeSolvers))

extractVecBLS(x::AbstractVector) = @view x[1:end-1]
extractVecBLS(x::BorderedArray)  = x.u

extractParBLS(x::AbstractVector) = x[end]
extractParBLS(x::BorderedArray)  = x.p

# We restrict to bordered systems where the added component is scalar
function (lbs::MatrixFreeBLS{S})(J, 	dR,
								dzu, 	dzp::T, R, n::T,
								ξu::Tξ = 1, ξp::Tξ = 1; shift = nothing) where {T <: Number, Tξ, S}
	~isnothing(shift) && @warn "Shift is not implemented for the bordered linear solver MatrixFreeBLS"
	linearmap = MatrixFreeBLSmap(J, dR, rmul!(copy(dzu), ξu), dzp * ξp, shift)
	# what is the vector type used?
	rhs = lbs.useBorderedArray ? BorderedArray(copy(R), n) : vcat(R, n)
	sol, cv, it = lbs.solver(linearmap, rhs)
	return extractVecBLS(sol), extractParBLS(sol), cv, it
end
####################################################################################################
# Linear Solvers based on a bordered solver
# !!!! This one is used as a linear Solver, not as a Bordered one
####################################################################################################
"""
$(TYPEDEF)

This structure is used to provide the following linear solver. To solve (1) J⋅x = rhs, one decomposes J using Matrix by blocks and then use a bordering strategy to solve (1).

$(TYPEDFIELDS)

!!! warn "Warning"
    The solver only works for `AbstractMatrix`
"""
struct LSFromBLS{Ts} <: AbstractLinearSolver
	"Linear solver used to solve the smaller linear systems."
	solver::Ts
end

LSFromBLS() = LSFromBLS(BorderingBLS(DefaultLS(useFactorization = false)))

function (l::LSFromBLS)(J, rhs)
	F = factorize(J[1:end-1, 1:end-1])
	x1, x2, flag, it = l.solver(F, Array(J[1:end-1,end]), Array(J[end,1:end-1]), J[end, end], (@view rhs[1:end-1]), rhs[end])
	return vcat(x1,x2), flag, sum(it)
end

function  (l::LSFromBLS)(J, rhs1, rhs2)
	F = factorize(J[1:end-1,1:end-1])
	x1, x2, flag1, it1 = l.solver(F, Array(J[1:end-1,end]), Array(J[end,1:end-1]), J[end, end], (@view rhs1[1:end-1]), rhs1[end])

	y1, y2, flag2, it2 = l.solver(F, Array(J[1:end-1,end]), Array(J[end,1:end-1]), J[end, end], (@view rhs2[1:end-1]), rhs2[end])

	return vcat(x1,x2), vcat(y1,y2), flag1 & flag2, (1, 1)
end


text:
using NamedDims
using NamedDims: names
using Test
using Statistics

@testset "Base" begin
    a = [10 20; 31 40]
    nda = NamedDimsArray(a, (:x, :y))

    @testset "$f" for f in (sum, prod, maximum, minimum, extrema)
        @test f(nda) == f(a)
        @test f(nda; dims=:x) == f(nda; dims=1) == f(a; dims=1)

        @test names(f(nda; dims=:x)) == (:x, :y) == names(f(nda; dims=1))
    end

    @testset "$f" for f in (cumsum, cumprod, sort)
        @test f(nda; dims=:x) == f(nda; dims=1) == f(a; dims=1)

        @test names(f(nda; dims=:x)) == (:x, :y) == names(f(nda; dims=1))

        @test f([1, 4, 3]) == f(NamedDimsArray([1, 4, 3], :vec))
        @test_throws UndefKeywordError f(nda)
        @test_throws UndefKeywordError f(a)
    end

    @testset "sort!" begin
        a = [1 9; 7 3]
        nda = NamedDimsArray(a, (:x, :y))

        # Vector case
        veca = [1, 9, 7, 3]
        sort!(NamedDimsArray(veca, :vec); order=Base.Reverse)
        @test issorted(veca; order=Base.Reverse)

        # Higher-dim case: `dims` keyword in `sort!` requires Julia v1.1+
        if VERSION > v"1.1-"
            sort!(nda, dims=:y)
            @test issorted(a[2, :])
            @test_throws UndefKeywordError sort!(nda)

            sort!(nda; dims=:x, order=Base.Reverse)
            @test issorted(a[:, 1]; order=Base.Reverse)
        end
    end

    @testset "$f!" for (f,f!) in zip((sum, prod, maximum, minimum), (sum!, prod!, maximum!, minimum!))
        a = [10 20; 31 40]
        nda = NamedDimsArray(a, (:x, :y)) # size (2,2)

        nda1 = sum(nda, dims=1)           # size (1,2)
        nda2 = sum(nda, dims=2)           # size (2,1)
        @testset "ndims==2" begin
            @test f!(nda1, nda) == f!(nda1, a) == f(a, dims=1)
            @test f!(nda2, nda) == f!(nda2, a) == f(a, dims=2)

            @test NamedDims.names(f!(nda1, nda)) == (:x, :y) == NamedDims.names(f!(nda1, a))
            @test NamedDims.names(f!(nda2, nda)) == (:x, :y) == NamedDims.names(f!(nda2, a))

            @test_throws DimensionMismatch f!(nda1, transpose(nda)) # names wrong way around
        end
        @testset "ndims==1 too" begin
            ndx = NamedDimsArray([3,4], :x)
            ndy = NamedDimsArray([5,6], :y)
            nd_ = NamedDimsArray([7,8], :_)

            @test f!(ndx, nda) == f!([0,0], nda) == dropdims(f(a, dims=2), dims=2)
            @test f!(nd_, nda) == f!(ndx, a)

            @test f!(ndy', nda) == f!([0 0], nda) == f(a, dims=1)

            @test_throws DimensionMismatch f!(ndy, nda) # name y on wrong dimension
        end
    end

    @testset "eachslice" begin
        if VERSION > v"1.1-"
            slices = [[111 121; 211 221], [112 122; 212 222]]
            a = cat(slices...; dims=3)
            nda = NamedDimsArray(a, (:a, :b, :c))

            @test (
                sum(eachslice(nda; dims=:c)) ==
                sum(eachslice(nda; dims=3)) ==
                sum(eachslice(a; dims=3)) ==
                slices[1] + slices[2]
            )
            @test_throws ArgumentError eachslice(nda; dims=(1, 2))
            @test_throws ArgumentError eachslice(a; dims=(1, 2))

            @test_throws UndefKeywordError eachslice(nda)
            @test_throws UndefKeywordError eachslice(a)

            @test (
                names(first(eachslice(nda; dims=:b))) ==
                names(first(eachslice(nda; dims=2))) ==
                (:a, :c)
            )
        end
    end

    @testset "mapslices" begin
        a = [10 20; 31 40]
        nda = NamedDimsArray(a, (:x, :y))

        @test (
            mapslices(join, nda; dims=:x) ==
            mapslices(join, nda; dims=1) ==
            ["1031" "2040"]
        )
        @test (
            mapslices(join, nda; dims=:y) ==
            mapslices(join, nda; dims=2) ==
            reshape(["1020", "3140"], Val(2))
        )
        @test (
            mapslices(join, nda; dims=(:x, :y)) ==
            mapslices(join, nda; dims=(1, 2)) ==
            reshape(["10312040"], (1, 1))
        )
        @test_throws UndefKeywordError mapslices(join, nda)
        @test_throws UndefKeywordError mapslices(join, a)

        @test (
            names(mapslices(join, nda; dims=:y)) ==
            names(mapslices(join, nda; dims=2)) ==
            (:x, :y)
        )
    end

    @testset "mapreduce" begin
        a = [10 20; 31 40]
        nda = NamedDimsArray(a, (:x, :y))

        @test mapreduce(isodd, |, nda) == true == mapreduce(isodd, |, a)
        @test (
            mapreduce(isodd, |, nda; dims=:x) ==
            mapreduce(isodd, |, nda; dims=1) ==
            [true false]
        )
        @test (
            mapreduce(isodd, |, nda; dims=:y) ==
            mapreduce(isodd, |, nda; dims=2) ==
            [false true]'
        )
        @test (
            names(mapreduce(isodd, |, nda; dims=:y)) ==
            names(mapreduce(isodd, |, nda; dims=2)) ==
            (:x, :y)
        )
    end

    @testset "zero" begin
        a = [10 20; 31 40]
        nda = NamedDimsArray(a, (:x, :y))

        @test zero(nda) == [0 0; 0 0] == zero(a)
        @test names(zero(nda)) == (:x, :y)
    end

    @testset "count" begin
        a = [true false; true true]
        nda = NamedDimsArray(a, (:x, :y))

        @test count(nda) == count(a) == 3
        @test_throws Exception count(nda; dims=:x)
        @test_throws Exception count(a; dims=1)
    end

    @testset "push!, pop!, etc" begin
        ndv = NamedDimsArray([10, 20, 30], (:i,))

        @test length(push!(ndv, 40)) == 4
        @test names(pushfirst!(ndv, 0)) == (:i,)
        @test ndv == 0:10:40

        @test pop!(ndv) == 40
        @test popfirst!(ndv) == 0
        @test ndv == [10, 20, 30]
    end

    @testset "append!, empty!" begin
        ndv = NamedDimsArray([10, 20, 30], (:i,))
        ndv45 = NamedDimsArray([40, 50], (:i,))
        ndv0 = NamedDimsArray([0, 0], (:zero,))

        @test length(append!(ndv, ndv45)) == 5
        @test names(append!(ndv, [60,70])) == (:i,)

        @test_throws DimensionMismatch append!(ndv, ndv0)
        @test ndv == 10:10:70 # error was thrown before altering

        @test names(empty!(ndv)) == (:i,)
        @test length(ndv) == 0
    end

    @testset "map, map!" begin
        nda = NamedDimsArray([11 12; 21 22], (:x, :y))

        @test names(map(+, nda, nda, nda)) == (:x, :y)
        @test names(map(+, nda, parent(nda), nda)) == (:x, :y)
        @test names(map(+, parent(nda), nda)) == (:x, :y)

        # this method only called based on first two arguments:
        @test names(map(+, parent(nda), parent(nda), nda)) == (:_, :_)

        # one-arg forms work without adding anything... except on 1.0...
        @test names(map(sqrt, nda)) == (:x, :y)
        @test foreach(sqrt, nda) === nothing

        # map! may return a different wrapper of the same data, like sum!
        semi = NamedDimsArray(rand(2,2), (:x, :_))
        @test names(map!(sqrt, rand(2,2), nda)) == (:x, :y)
        @test names(map!(sqrt, semi, nda)) == (:x, :y)

        zed = similar(nda, Float64)
        @test map!(sqrt, zed, nda) == sqrt.(nda)
        @test zed[1,1] == sqrt(nda[1,1])

        # mismatching names
        @test_throws DimensionMismatch map(+, nda, transpose(nda))
        @test_throws DimensionMismatch map(+, nda, parent(nda), nda, transpose(nda))
        @test_throws DimensionMismatch map!(sqrt, semi, transpose(nda))

        @test foreach(+, semi, nda) === nothing
        @test_throws DimensionMismatch foreach(+, semi, transpose(nda))
    end

    @testset "filter" begin
        nda = NamedDimsArray([11 12; 21 22], (:x, :y))
        ndv = NamedDimsArray(1:7, (:z,))

        @test names(filter(isodd, ndv)) == (:z,)
        @test names(filter(isodd, nda)) == (:_,)
    end

    @testset "collect(generator)" begin
        nda = NamedDimsArray([11 12; 21 22], (:x, :y))
        ndv = NamedDimsArray([10, 20, 30], (:z,))

        @test names([sqrt(x) for x in nda]) == (:x, :y)

        @test names([x^i for (i,x) in enumerate(ndv)]) == (:z,)
        @test names([x^i for (i,x) in enumerate(nda)]) == (:x, :y)

        # Iterators.product -- has all names
        @test names([x+y for x in nda, y in ndv]) == (:x, :y, :z)
        @test names([x+y for x in nda, y in 1:5]) == (:x, :y, :_)
        @test names([x+y for x in 1:5, y in ndv]) == (:_, :z)
        four = [x*y/z^p for p in 1:2, x in ndv, y in 1:2, z in nda]
        @test names(four) == (:_, :z, :_, :x, :y)

        # Iterators.flatten -- no obvious name to use
        @test names([x+y for x in nda for y in ndv]) == (:_,)

        if VERSION >= v"1.1"
            # can't see inside eachslice generators, until:
            # https://github.com/JuliaLang/julia/pull/32310
            @test names([sum(c) for c in eachcol(nda)]) == (:_,)
        end
    end

end  # Base

@testset "Statistics" begin
    a = [10 20; 30 40]
    nda = NamedDimsArray(a, (:x, :y))
    @testset "$f" for f in (mean, std, var, median)
        @test f(nda) == f(a)
        @test f(nda; dims=:x) == f(nda; dims=1) == f(a; dims=1)

        @test names(f(nda; dims=:x)) == (:x, :y) == names(f(nda; dims=1))
    end
end


text:
module TestPlotting

using Mimi
using Test

using Mimi: plot_comp_graph

@defcomp LongComponent begin
    x = Parameter(index=[time])
    y = Parameter()
    z = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        v.z[ts] = p.x[ts] + p.y
    end
end

@defcomp ShortComponent begin
    a = Parameter()
    b = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        v.b[ts] = p.a * ts.t
    end
end

m = Model()
set_dimension!(m, :time, 2000:3000)
nsteps = Mimi.dim_count(m.md, :time)

add_comp!(m, ShortComponent) #; first=2100)
add_comp!(m, LongComponent) #; first=2000)

set_param!(m, :ShortComponent, :a, 2.)
set_param!(m, :LongComponent, :y, 1.)
connect_param!(m, :LongComponent, :x, :ShortComponent, :b, zeros(nsteps))

run(m)

graph = plot_comp_graph(m)
@test typeof(graph) == Mimi.Compose.Context

end #module


text:
using CuArrays, Test

include("util.jl")

using Random
Random.seed!(1)

using CUDAnative
using CUDAdrv

# GPUArrays has a testsuite that isn't part of the main package.
# Include it directly.
import GPUArrays
gpuarrays = pathof(GPUArrays)
gpuarrays_root = dirname(dirname(gpuarrays))
include(joinpath(gpuarrays_root, "test", "testsuite.jl"))

testf(f, xs...; kwargs...) = TestSuite.compare(f, CuArray, xs...; kwargs...)

# pick a suiteable device (by available memory,
# but also by capability if testing needs to be thorough)
candidates = [(dev=dev,
                cap=capability(dev),
                mem=CuContext(ctx->CUDAdrv.available_memory(), dev))
              for dev in devices()]
thorough = parse(Bool, get(ENV, "CI_THOROUGH", "false"))
if thorough
    sort!(candidates, by=x->(x.cap, x.mem))
else
    sort!(candidates, by=x->x.mem)
end
pick = last(candidates)
@info("Testing using device $(name(pick.dev)) (compute capability $(pick.cap), $(Base.format_bytes(pick.mem)) available memory) on CUDA driver $(CUDAdrv.version()) and toolkit $(CUDAnative.version())")
device!(pick.dev)

@testset "CuArrays" begin

CuArrays.allowscalar(false)

CuArrays.enable_timings()

@testset "GPUArrays test suite" begin
  TestSuite.test(CuArray)
end

include("base.jl")
include("memory.jl")
include("blas.jl")
include("rand.jl")
include("fft.jl")
include("sparse.jl")
include("solver.jl")
include("sparse_solver.jl")
include("dnn.jl")
include("tensor.jl")
include("forwarddiff.jl")

if haskey(ENV, "CI")
  CuArrays.memory_status()
  CuArrays.pool_timings()
  CuArrays.alloc_timings()
end

CuArrays.reset_timers!()

end


text:
np = nr * (ns + 4) + 1
p = randn(Float64, np) .* 1.e-2;
p[1:nr] .+= 0.8;  # w_b
p[nr*(ns+1)+1:nr*(ns+2)] .+= 0.8;  # w_out
p[nr*(ns+2)+1:nr*(ns+4)] .+= 0.1;  # w_b | w_Ea
p[end] = 0.1;  # slope

function p2vec(p)
    slope = p[end] .* 1.e1
    w_b = p[1:nr] .* (slope * 10.0)
    w_b = clamp.(w_b, 0, 50)

    w_out = reshape(p[nr+1:nr*(ns+1)], ns, nr)
    @. w_out[1, :] = clamp(w_out[1, :], -3.0, 0.0)
    @. w_out[end, :] = clamp(abs(w_out[end, :]), 0.0, 3.0)

    if p_cutoff > 0.0
        w_out[findall(abs.(w_out) .< p_cutoff)] .= 0.0
    end

    w_out[ns-1:ns-1, :] .=
        -sum(w_out[1:ns-2, :], dims = 1) .- sum(w_out[ns:ns, :], dims = 1)

    w_in_Ea = abs.(p[nr*(ns+1)+1:nr*(ns+2)] .* (slope * 100.0))
    w_in_Ea = clamp.(w_in_Ea, 0.0, 300.0)

    w_in_b = abs.(p[nr*(ns+2)+1:nr*(ns+3)])

    w_in_ocen = abs.(p[nr*(ns+3)+1:nr*(ns+4)])
    w_in_ocen = clamp.(w_in_ocen, 0.0, 1.5)

    if p_cutoff > 0.0
        w_in_ocen[findall(abs.(w_in_ocen) .< p_cutoff)] .= 0.0
    end

    w_in = vcat(clamp.(-w_out, 0.0, 4.0), w_in_Ea', w_in_b', w_in_ocen')
    return w_in, w_b, w_out
end

function display_p(p)
    w_in, w_b, w_out = p2vec(p)
    println("\n species (column) reaction (row)")
    println("w_in | Ea | b | n_ocen | lnA | w_out")
    show(stdout, "text/plain", round.(hcat(w_in', w_b, w_out'), digits = 2))
    # println("\n w_out")
    # show(stdout, "text/plain", round.(w_out', digits=3))
    println("\n")
end
# display_p(p);

function getsampletemp(t, T0, beta)
    if beta < 100
        T = T0 + beta / 60 * t  # K/min to K/s
    else
        tc = [999.0, 1059.0] .* 60.0
        Tc = [beta, 370.0, 500.0] .+ 273.0
        HR = 40.0 / 60.0
        if t <= tc[1]
            T = Tc[1]
        elseif t <= tc[2]
            T = minimum([Tc[1] + HR * (t - tc[1]), Tc[2]])
        else
            T = minimum([Tc[2] + HR * (t - tc[2]), Tc[3]])
        end
    end
    return T
end

const R = -1.0 / 8.314e-3  # universal gas constant, kJ/mol*K
@inbounds function crnn!(du, u, p, t)
    logX = @. log(clamp(u, lb, 10.0))
    T = getsampletemp(t, T0, beta)
    w_in_x = w_in' * vcat(logX, R / T, log(T), ocen)
    du .= w_out * (@. exp(w_in_x + w_b))
end

tspan = [0.0, 1.0];
u0 = zeros(ns);
u0[1] = 1.0;
prob = ODEProblem(crnn!, u0, tspan, p, abstol = lb)

condition(u, t, integrator) = u[1] < lb * 10.0
affect!(integrator) = terminate!(integrator)
_cb = DiscreteCallback(condition, affect!)

alg = AutoTsit5(TRBDF2(autodiff = true));
# sense = BacksolveAdjoint(checkpointing=true; autojacvec=ZygoteVJP());
sense = ForwardSensitivity(autojacvec = true)
function pred_n_ode(p, i_exp, exp_data)
    global T0, beta, ocen = l_exp_info[i_exp, :]
    global w_in, w_b, w_out = p2vec(p)
    ts = @view(exp_data[:, 1])
    tspan = [ts[1], ts[end]]
    sol = solve(
        prob,
        alg,
        tspan = tspan,
        p = p,
        saveat = ts,
        sensalg = sense,
        maxiters = maxiters,
        # callback = _cb,
    )

    if sol.retcode == :Success
        nothing
    else
        @sprintf("solver failed beta: %.0f ocen: %.2f", beta, exp(ocen))
    end
    if length(sol.t) > length(ts)
        # @show exp_data[:, 1]
        # @show sol.t
        return  sol[:, 1:length(ts)]
    else
        return sol
    end
end

function loss_neuralode(p, i_exp)
    exp_data = l_exp_data[i_exp]
    pred = Array(pred_n_ode(p, i_exp, exp_data))
    masslist = sum(clamp.(@view(pred[1:end-1, :]), 0, Inf), dims = 1)'
    gaslist = clamp.(@views(pred[end, :]), 0, Inf)

    loss = mae(masslist, @view(exp_data[1:length(masslist), 3]))
    if ocen < 1000.0
        loss += mae(gaslist, 1 .- @view(exp_data[1:length(masslist), 3]))
    end
    return loss
end
@time loss = loss_neuralode(p, 1)
# using BenchmarkTools
# @benchmark loss = loss_neuralode(p, 1)
# @benchmark grad = ForwardDiff.gradient(x -> loss_neuralode(x, 1), p)


text:
const UINT64_MAX = typemax(UInt64)
const SIZE_MAX = typemax(UInt32)
const crypto_pwhash_argon2i_MEMLIMIT_MAX = ( ( SIZE_MAX >= 4398046510080 ) ? 4398046510080 : ( SIZE_MAX >= Cuint(2147483648) ) ? Cuint(2147483648) : Cuint(32768) )
const crypto_pwhash_argon2id_MEMLIMIT_MAX = ( ( SIZE_MAX >= 4398046510080 ) ? 4398046510080 : ( SIZE_MAX >= Cuint(2147483648) ) ? Cuint(2147483648) : Cuint(32768) )


text:
"""
Calculate value, gradient, and hessian of the variational ELBO.
"""
module DeterministicVI

using Base.Threads: threadid, nthreads

import ..Config
using ..BivariateNormals: BivariateNormalDerivatives, BvnComponent,
                          GalaxySigmaDerivs, get_bvn_cov, eval_bvn_pdf!,
                          get_bvn_derivs!, transform_bvn_derivs!
using ..Model
using ..Model: SkyPatch, BvnBundle, populate_fsm!
import ..Celeste: Const, @aliasscope, @unroll_loop
import ..Infer
using ..SensitiveFloats
import ..Log
using ..Transform
import DataFrames
import Optim
import ForwardDiff.Dual
using StaticArrays
import Base: convert

export ElboArgs, generic_init_source, catalog_init_source, init_sources,
       VariationalParams, elbo, ElboIntermediateVariables

function init_thread_pool!(pool::Vector, create)
    if length(pool) != Base.Threads.nthreads()
        empty!(pool)
        for i = 1:Base.Threads.nthreads()
            push!(pool, create())
        end
    end
end

"""
Return a default-initialized VariationalParams instance.
"""
function generic_init_source(init_pos::Vector{Float64})
    ret = Vector{Float64}(length(CanonicalParams))
    ret[ids.is_star] = 0.5
    ret[ids.pos] = init_pos
    ret[ids.flux_loc] = log(2.0)
    ret[ids.flux_scale] = 1e-3
    ret[ids.gal_fracdev] = 0.5
    ret[ids.gal_ab] = 0.5
    ret[ids.gal_angle] = 0.0
    ret[ids.gal_scale] = 1.0
    ret[ids.k] = 1.0 / size(ids.k, 1)
    ret[ids.color_mean] = 0.0
    ret[ids.color_var] =  1e-2
    ret
end

"""
Return VariationalParams instance initialized form a catalog entry
"""
function catalog_init_source(ce::CatalogEntry; max_gal_scale=Inf)
    # TODO: sync this up with the transform bounds
    ret = generic_init_source(ce.pos)

    # TODO: don't do this thresholding for background sources,
    # just for sources that are being optimized
    ret[ids.is_star[1]] = ce.is_star ? 0.8: 0.2
    ret[ids.is_star[2]] = ce.is_star ? 0.2: 0.8

    ret[ids.flux_loc[1]] = log(max(0.1, ce.star_fluxes[3]))
    ret[ids.flux_loc[2]] = log(max(0.1, ce.gal_fluxes[3]))

    function get_color(color_var, color_mean)
        color_var > 0 && color_mean > 0 ? min(max(log(color_var / color_mean), -9.), 9.) :
            color_var > 0 && color_mean <= 0 ? 3.0 :
                color_var <= 0 && color_mean > 0 ? -3.0 : 0.0
    end

    function get_colors(raw_fluxes)
        [get_color(raw_fluxes[c+1], raw_fluxes[c]) for c in 1:4]
    end

    ret[ids.color_mean[:, 1]] = get_colors(ce.star_fluxes)
    ret[ids.color_mean[:, 2]] = get_colors(ce.gal_fluxes)

    ret[ids.gal_fracdev] = min(max(ce.gal_frac_dev, 0.015), 0.985)

    ret[ids.gal_ab] = ce.is_star ? .8 : min(max(ce.gal_ab, 0.015), 0.985)
    ret[ids.gal_angle] = ce.gal_angle
    ret[ids.gal_scale] = ce.is_star ? 0.2 : min(max_gal_scale, max(ce.gal_scale, 0.2))

    ret
end

function init_sources(target_sources::Vector{Int}, catalog::Vector{CatalogEntry})
    ret = Vector{Vector{Float64}}(length(catalog))
    for s in 1:length(catalog)
        ret[s] = catalog_init_source(catalog[s])
    end
    for s in target_sources
        ret[s][:] = generic_init_source(catalog[s].pos)
    end
    ret
end

include("deterministic_vi/elbo_args.jl")
include("deterministic_vi/elbo_kl.jl")
include("deterministic_vi/source_brightness.jl")
include("deterministic_vi/elbo_objective.jl")
include("deterministic_vi/ConstraintTransforms.jl")
include("deterministic_vi/ElboMaximize.jl")

"""
Infers one light source. This routine is intended to be called in parallel,
once per target light source.

Arguments:
    images: a collection of astronomical images
    neighbors: the other light sources near `entry`
    entry: the source to infer
"""
function infer_source(config::Config,
                      images::Vector{Image},
                      neighbors::Vector{CatalogEntry},
                      entry::CatalogEntry)
    if length(neighbors) > 100
        msg = string("object at RA, Dec = $(entry.pos) has an excessive",
                     "number ($(length(neighbors))) of neighbors")
        Log.warn(msg)
    end

    # It's a bit inefficient to call the next 5 lines every time we optimize_f.
    # But, as long as runtime is dominated by the call to maximize!, that
    # isn't a big deal.
    cat_local = vcat([entry], neighbors)
    vp = init_sources([1], cat_local)
    patches = Infer.get_sky_patches(images, cat_local)
    Infer.load_active_pixels!(config, images, patches)

    ea = ElboArgs(images, patches, [1])
    f_evals, max_f, max_x, nm_result = ElboMaximize.maximize!(ea, vp)
    return vp[1]
end

# legacy wrapper
function infer_source(images::Vector{Image},
                      neighbors::Vector{CatalogEntry},
                      entry::CatalogEntry)
    infer_source(Config(), images, neighbors, entry)
end

end


text:
export ReflectGate, reflect

"""
    ReflectGate{N, T, Tr} <: PrimitiveBlock{N}

Reflection operator to target state `psi`.

# Definition

```math
|ψ⟩ → 2 |s⟩⟨s| - 1
```
"""
struct ReflectGate{N, T, Tr <: ArrayReg{1, T}} <: PrimitiveBlock{N}
    psi::Tr
end

"""
    ReflectGate(r::ArrayReg{1})

Create a [`ReflectGate`](@ref) with a quantum register `r`.
"""
ReflectGate(r::ArrayReg{1, T}) where T = ReflectGate{nqubits(r), T, typeof(r)}(r)

"""
    ReflectGate(r::AbstractVector)

Create a [`ReflectGate`](@ref) with a quantum state vector `v`.
"""
ReflectGate(v::AbstractVector{<:Complex}) = ReflectGate(ArrayReg(v))

"""
    reflect(r::ArrayReg)

Create a [`ReflectGate`](@ref) with an [`ArrayReg`](@ref).
"""
reflect(r::ArrayReg) = reflect(statevec(r))

"""
    reflect(v::AbstractVector{<:Complex})

Create a [`ReflectGate`](@ref) with an quantum state vector `v`.

# Example

```jldoctest; setup=:(using YaoBlocks; using YaoArrayRegister)
julia> reflect(rand_state(3))
┌ Warning: `ReflectGate` will be moved to `YaoExtensions.jl` in the next release.
└ @ YaoBlocks ~/.julia/packages/YaoBlocks/XXC9Y/src/primitive/reflect_gate.jl:19
reflect(ArrayReg{1, Complex{Float64}, Array...})
```
"""
reflect(v::AbstractVector{<:Complex}) = ReflectGate(v)

function apply!(r::ArrayReg, g::ReflectGate{N, T, <:ArrayReg}) where {N, T}
    v = state(g.psi)
    r.state .= 2 .* (v' * r.state) .* v - r.state
    return r
end

# target type is the same with block's
function mat(::Type{T}, r::ReflectGate{N, T}) where {N, T}
    v = statevec(r.psi)
    return 2 * v * v' - IMatrix(length(v))
end

# different
function mat(::Type{T1}, r::ReflectGate{N, T2}) where {N, T1, T2}
    M = mat(T2, r)
    return copyto!(similar(M, T1), M)
end

Base.:(==)(A::ReflectGate, B::ReflectGate) = A.psi == B.psi
Base.copy(r::ReflectGate) = ReflectGate(r.psi)

Yao.isreflexive(::ReflectGate) = true
Yao.ishermitian(::ReflectGate) = true
Yao.isunitary(::ReflectGate) = true


text:
using MacroTools

include("transformed_distribution.jl")
include("relabeled_distribution.jl")

#TODO: Remove Distribution{T} from structs; use a type U instead

# `Arg` types: represent arguments to distributions
abstract type Arg end

# Represents a argument at a certain index
struct SimpleArg <: Arg
    i :: Int8
end

struct TransformedArg <: Arg
    f_args :: Vector{Arg} # Is it better to do Union{TransformedArg, SimpleArg}?
    orig_f :: Function
    arg_passer :: Function
end

# Each j is an Arg, and has a yes/no answer about whether
# it supports gradients. The Arg may correspond to multiple i's,
# which are user-facing parameters.
all_indices(arg::SimpleArg) = [arg.i]
all_indices(arg::TransformedArg) = vcat([all_indices(a) for a in arg.f_args]...)

eval_arg(x::Any, args) = x
eval_arg(x::SimpleArg, args) = args[x.i]
eval_arg(x::TransformedArg, args) = x.arg_passer(x.orig_f, [eval_arg(a, args) for a in x.f_args]...)

# DistWithArgs
struct DistWithArgs{T}
    base :: Distribution{T}
    arglist # Contains Args and other values
end

struct CompiledDistWithArgs{T} <: Distribution{T}
    base :: Distribution{T}
    n_args :: Int8
    arg_grad_bools
    arglist
end

function compile_dist_with_args(d::DistWithArgs{T}, n::Int8)::CompiledDistWithArgs{T} where T
    base_arg_grads = has_argument_grads(d.base)
    arg_grad_bools = fill(true, n)
    for (i, arg) in enumerate(d.arglist)
        if typeof(arg) <: Arg
            for idx in all_indices(arg)
                arg_grad_bools[idx] = arg_grad_bools[idx] && base_arg_grads[i]
            end
        end
    end
    CompiledDistWithArgs{T}(d.base, n, arg_grad_bools, d.arglist)
end

function logpdf(d::CompiledDistWithArgs{T}, x::T, args...) where T
    concrete_args = [eval_arg(arg, args) for arg in d.arglist]
    logpdf(d.base, x, concrete_args...)
end

function logpdf_grad(d::CompiledDistWithArgs{T}, x::T, args...) where T
    self_has_output_grad = has_output_grad(d)
    self_has_arg_grads = has_argument_grads(d)
    concrete_args = [eval_arg(arg, args) for arg in d.arglist]
    base_has_arg_grads = has_argument_grads(d.base)
    base_grads = logpdf_grad(d.base, x, concrete_args...)

    base_arg_grads = [g for (i, g) in enumerate(base_grads[2:end]) if base_has_arg_grads[i]]
    argvec = collect(args)
    eval_arg_grads = hcat([ReverseDiff.gradient(xs -> eval_arg(arg, xs), argvec) for (i, arg) in enumerate(d.arglist) if base_has_arg_grads[i]]...)

    retval = [base_grads[1]]
    for i in 1:d.n_args
        if self_has_arg_grads[i]
            push!(retval, eval_arg_grads[i,:]' * base_arg_grads)
        else
            push!(retval, nothing)
        end
    end
    retval
end

function random(d::CompiledDistWithArgs{T}, args...)::T where T
    concrete_args = [eval_arg(arg, args) for arg in d.arglist]
    random(d.base, concrete_args...)
end

is_discrete(d::CompiledDistWithArgs{T}) where T = is_discrete(d.base)

(d::CompiledDistWithArgs{T})(args...) where T = random(d, args...)

function has_output_grad(d::CompiledDistWithArgs{T}) where T
    has_output_grad(d.base)
end

function has_argument_grads(d::CompiledDistWithArgs{T}) where T
    d.arg_grad_bools
end

# ACTUAL DSL

# A function call within a @dist body
function dist_call(d::Distribution{T}, args) where T
    DistWithArgs{T}(d, args)
end
function dist_call(f, args)
    if any(x -> (typeof(x) <: DistWithArgs), args)
        f(args...)
    elseif any(x -> (typeof(x) <: Arg), args)
        # Create a transformed distribution
        f_args = []
        new_f_arg_names = []
        actual_arg_exprs = []

        for (i, x) in enumerate(args)
            if typeof(x) <: Arg
                push!(f_args, x)
                name = gensym()
                push!(new_f_arg_names, name)
                push!(actual_arg_exprs, name)
            else
                push!(actual_arg_exprs, Meta.quot(x))
            end
        end
        arg_passer = eval(:((f, $(new_f_arg_names...)) -> f($(actual_arg_exprs...))))
        TransformedArg(f_args, f, arg_passer)
    else
        f(args...)
    end
end

#TODO: Make this actually compile a new type
macro dist(fnexpr)
  # First, pull out arguments and body
  fndef = splitdef(longdef(fnexpr))
  #
  # if !@capture(fnexpr, argexpr_ -> body_)
  #     error("@dist expression must be of the form @dist (args...) -> body")
  # end
  # arguments = (typeof(argexpr) == Symbol) ? (argexpr,) : argexpr.args # [] #argexpr

  arguments = fndef[:args]
  body = fndef[:body]

  name_to_index = Dict{Symbol, Int8}()
  for (i, arg) in enumerate(arguments)
      (argname, _, _, _) = splitarg(arg)
      name_to_index[argname] = i
  end

  function process_node(node)
      if typeof(node) == Symbol && haskey(name_to_index, node)
          return :($(SimpleArg(name_to_index[node])))
      end
      if @capture(node, f_(xs__)) && f != :dist_call
          return :(Gen.dist_call($(f), ($(xs...),)))
      end
      node
  end

  dwa_expr = MacroTools.postwalk(process_node, body)
  :($(esc(fndef[:name])) = compile_dist_with_args($(esc(dwa_expr)), Int8($(length(arguments)))))
end

# Addition
Base.:+(b::DistWithArgs{T}, a::Real) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, x -> x + a, x -> x - a, x -> (1.0,)), b.arglist)
Base.:+(b::DistWithArgs{T}, a::Arg) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, +, -, (x, a) -> (1.0, -1.0)), (a, b.arglist...))
Base.:+(a::Real, b::DistWithArgs{T}) where T <: Real = b + a
Base.:+(a::Arg, b::DistWithArgs{T}) where T <: Real = b + a

# Subtraction
Base.:-(b::DistWithArgs{T}, a::Real) where T <: Real = b + (-1 * a)
Base.:-(b::DistWithArgs{T}, a::Arg) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, -, +, (x, a) -> (1.0, 1.0)), (a, b.arglist...))
Base.:-(a::Real, b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, x -> a - x, x -> a - x, x -> (-1.0,)), b.arglist)
Base.:-(a::Arg, b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, (x, a) -> a - x, (x, a) -> a - x, (x, a) -> (-1.0, 1.0)), (a, b.arglist...))

# Multiplication
Base.:*(b::DistWithArgs{T}, a::Real) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, x -> x * a, x -> x / a, x -> (1.0/a,)), b.arglist)
Base.:*(b::DistWithArgs{T}, a::Arg) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, *, /, (x, a) -> (1.0/a, -x/(a*a))), (a, b.arglist...))
Base.:*(a::Real, b::DistWithArgs{T}) where T <: Real = b * a
Base.:*(a::Arg, b::DistWithArgs{T}) where T <: Real = b * a

# Division
Base.:/(b::DistWithArgs{T}, a::Real) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, x -> x / a, x -> x * a, x -> (a,)), b.arglist)
Base.:/(b::DistWithArgs{T}, a::Arg) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, /, *, (x, a) -> (a, x)), (a, b.arglist...))
Base.:/(a::Real, b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, x -> a / x, x -> a / x, x -> (-a / (x*x),)), b.arglist)
Base.:/(a::Arg, b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 1, (x, a) -> a/x, (x, a) -> a/x, (x, a) -> (-a/(x*x), 1.0/x)), (a, b.arglist...))

# Exponentiation
Base.exp(b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, exp, log, x -> (1.0 / x,)), b.arglist)
Base.log(b::DistWithArgs{T}) where T <: Real = DistWithArgs(TransformedDistribution{T, T}(b.base, 0, log, exp, x -> (exp(x),)), b.arglist)

# Indexing
Base.getindex(collection::Arg, d::DistWithArgs{T}) where T = DistWithArgs{Any}(WithLabelArg{Any, T}(d.base), (collection, d.arglist...))
Base.getindex(collection::AbstractArray{T}, d::DistWithArgs{U}) where {T, U} = DistWithArgs{T}(RelabeledDistribution{T, U}(d.base, collection), d.arglist)
Base.getindex(collection::AbstractDict{U, T}, d::DistWithArgs{U}) where {T, U} = DistWithArgs{T}(RelabeledDistribution{T, U}(d.base, collection), d.arglist)
# Ensure no ambiguity with `Base` implementation.
Base.getindex(collection::Dict{U, T}, d::DistWithArgs{U}) where {T, U} = DistWithArgs{T}(RelabeledDistribution{T, U}(d.base, collection), d.arglist)

# `Enum` constructors
function (::Type{T})(d::DistWithArgs{U}) where {T <: Enum, U}
    lookup = Dict(Int(i) => i for i in instances(T))
    DistWithArgs{T}(RelabeledDistribution{T, U}(d.base, lookup), d.arglist)
end

export @dist

max_stars_repo_name: cameronraysmith/Gen

max_issues_repo_name: cameronraysmith/Gen

max_forks_repo_name: cameronraysmith/Gen


text:
function solve()
    n, k = [parse(Int, x) for x in split(readline())]
    triked = ones(Bool, n)
    for _ in 1:k
        d = parse(Int, readline())
        a = [parse(Int, x) for x in split(readline())]
        for aa in a
            triked[aa] = false
        end
    end
    sum(triked)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
import ProgressMeter
#######################################################################################
# Original correlation sum
#######################################################################################
using Distances, Roots
export correlationsum, grassberger_dim, boxed_correlationsum,
estimate_r0_buenoorovio, data_boxing, autoprismdim, estimate_r0_theiler

"""
    correlationsum(X, ε::Real; w = 0, norm = Euclidean(), q = 2) → C_q(ε)
Calculate the `q`-order correlation sum of `X` (`Dataset` or timeseries)
for a given radius `ε` and `norm`. They keyword `show_progress = false` can be used
to display a progress bar for large `X`.

The function [`boxed_correlationsum`](@ref) is faster and should be preferred over this one.

## Description
The correlation sum is done using the formula:
```math
C_2(\\epsilon) = \\frac{2}{(N-w)(N-w-1)}\\sum_{i=1}^{N}\\sum_{j=1+w+i}^{N} B(||X_i - X_j|| < \\epsilon)
```
for `q=2` and
```math
C_q(\\epsilon) = \\left[\\frac{1}{\\alpha} \\sum_{i=w+1}^{N-w}\\left[\\sum_{j:|i-j| > w} B(||X_i - X_j|| < \\epsilon)\\right]^{q-1}\\right]^{1/(q-1)}
```
where
```math
\\alpha = (N-2w)(N-2w-1)^{(q-1)}
```
for `q≠2`, where ``N`` is its length and ``B`` gives 1 if the argument is
`true`. `w` is the [Theiler window](@ref). If `ε` is a vector its values have to be
ordered. See the article of Grassberger for the general definition [^Grassberger2007] and the book "Nonlinear Time Series Analysis" [^Kantz2003], Ch. 6, for
a discussion around `w` and choosing best values and Ch. 11.3 for the
explicit definition of the q-order correlationsum.

    correlationsum(X, εs::AbstractVector; w, norm, q) → C_q(ε)

If `εs` is a vector, `C_q` is calculated for each `ε ∈ εs`.
If also `q=2`, we attempt to do further optimizations are done, if the allocation
a matrix of size `N×N` is possible

See [`grassberger`](@ref) for more. See also [`takens_best_estimate`](@ref).

[^Grassberger]: Peter Grassberger (2007) [Grassberger-Procaccia algorithm. Scholarpedia, 2(5):3043.](http://dx.doi.org/10.4249/scholarpedia.3043)

[^Kantz]: Kantz, H., & Schreiber, T. (2003). [More about invariant quantities. In Nonlinear Time Series Analysis (pp. 197-233). Cambridge: Cambridge University Press.](https://doi:10.1017/CBO9780511755798.013)
"""
function correlationsum(X, ε; q = 2, norm = Euclidean(), w = 0, show_progress = false)
    if q == 2
        correlationsum_2(X, ε, norm, w, show_progress)
    else
        correlationsum_q(X, ε, q, norm, w, show_progress)
    end
end

function correlationsum_2(X, ε::Real, norm, w, show_progress)
    N = length(X)
    if show_progress
        progress = ProgressMeter.Progress(N; desc = "Correlation sum: ", dt = 1.0)
    end
    C = zero(eltype(X))
    @inbounds for (i, x) in enumerate(X)
        for j in i+1+w:N
            C += evaluate(norm, x, X[j]) < ε
        end
        show_progress && ProgressMeter.next!(progress)
    end
    return C * 2 / ((N-w-1)*(N-w))
end

function correlationsum_q(X, ε::Real, q, norm, w, show_progress)
    N, C = length(X), zero(eltype(X))
    normalisation = (N-2w)*(N-2w-one(eltype(X)))^(q-1)
    if show_progress
        progress = ProgressMeter.Progress(length(1+w:N-w); desc = "Correlation sum: ", dt = 1.0)
    end
    for i in 1+w:N-w
        x = X[i]
        C_current = zero(eltype(X))
        # computes all distances from 0 up to i-w
        for j in 1:i-w-1
            C_current += evaluate(norm, x, X[j]) < ε
        end
        # computes all distances after i+w till the end
        for j in i+w+1:N
            C_current += evaluate(norm, x, X[j]) < ε
        end
        C += C_current^(q - 1)
        show_progress && ProgressMeter.next!(progress)
    end
    return (C / normalisation) ^ (1 / (q-1))
end

# Optimized version
function correlationsum_2(X, εs::AbstractVector, norm, w, show_progress)
    @assert issorted(εs) "Sorted εs required for optimized version."
    d = try
        distancematrix(X, norm)
    catch err
        @warn "Couldn't create distance matrix ($(typeof(err))). Using slower algorithm..."
        return [correlationsum_2(X, ε, norm, w, show_progress) for ε in εs]
    end
    return correlationsum_2_fb(X, εs, d, w, show_progress) # function barrier
end
function correlationsum_2_fb(X, εs, d, w, show_progress)
    Cs = zeros(eltype(X), length(εs))
    N = length(X)
    factor = 2/((N-w)*(N-1-w))
    if show_progress
        K = length(length(εs)÷2:-1:1)
        M = K + length((length(εs)÷2 + 1):length(εs))
        progress = ProgressMeter.Progress(M; desc = "Correlation sum: ", dt = 1.0)
    end

    # First loop: mid-way ε until lower saturation point (C=0)
    for (ki, k) in enumerate(length(εs)÷2:-1:1)
        ε = εs[k]
        for i in 1:N
            @inbounds Cs[k] += count(d[j, i] < ε for j in i+1+w:N)
        end
        show_progress && ProgressMeter.update!(progress, ki)
        Cs[k] == 0 && break
    end
    # Second loop: mid-way ε until higher saturation point (C=max)
    for (ki, k) in enumerate((length(εs)÷2 + 1):length(εs))
        ε = εs[k]
        for i in 1:N
            @inbounds Cs[k] += count(d[j, i] < ε for j in i+1+w:N)
        end
        show_progress && ProgressMeter.update!(progress, ki+K)
        if Cs[k] ≈ 1/factor
            Cs[k:end] .= 1/factor
            break
        end
    end
    show_progress && ProgressMeter.finish!(progress)
    return Cs .* factor
end

function correlationsum_q(X, εs::AbstractVector, q, norm, w, show_progress)
    @assert issorted(εs) "Sorted εs required for optimized version."
    Nε, T, N = length(εs), eltype(X), length(X)
    Cs = zeros(T, Nε)
    normalisation = (N-2w)*(N-2w-one(T))^(q-1)
    if show_progress
        progress = ProgressMeter.Progress(length(1+w:N-w); desc = "Correlation sum: ", dt = 1.0)
    end
    for i in 1+w:N-w
        x = X[i]
        C_current = zeros(T, Nε)
        # Compute distances for j outside the Theiler window
        for j in Iterators.flatten((1:i-w-1, i+w+1:N))
            dist = evaluate(norm, x, X[j])
            for k in Nε:-1:1
                if dist < εs[k]
                    C_current[k] += 1
                else
                    break
                end
            end
        end
        Cs .+= C_current .^ (q-1)
        show_progress && ProgressMeter.next!(progress)
    end
    return (Cs ./ normalisation) .^ (1/(q-1))
end

function distancematrix(X, norm = Euclidean())
    N = length(X)
    d = zeros(eltype(X), N, N)
    @inbounds for i in 1:N
        for j in i+1:N
            d[j, i] = evaluate(norm, X[i], X[j])
        end
    end
    return d
end

"""
    grassberger_dim(data, εs = estimate_boxsizes(data); kwargs...) → D_C
Use the method of Grassberger and Proccacia[^Grassberger1983], and the correction by
Theiler[^Theiler1986], to estimate the correlation dimension `D_C` of the given `data`.

This function does something extremely simple:
```julia
cm = correlationsum(data, εs; kwargs...)
return linear_region(log.(sizes), log(cm))[2]
```
i.e. it calculates [`correlationsum`](@ref) for various radii and then tries to find
a linear region in the plot of the log of the correlation sum versus log(ε).
See [`generalized_dim`](@ref) for a more thorough explanation.

See also [`takens_best_estimate`](@ref).

[^Grassberger1983]: Grassberger and Proccacia, [Characterization of strange attractors, PRL 50 (1983)](https://journals-aps-org.e-bis.mpimet.mpg.de/prl/abstract/10.1103/PhysRevLett.50.346)

[^Theiler1986]: Theiler, [Spurious dimension from correlation algorithms applied to limited time-series data. Physical Review A, 34](https://doi.org/10.1103/PhysRevA.34.2427)
"""
function grassberger_dim(data::AbstractDataset, εs = estimate_boxsizes(data); kwargs...)
    cm = correlationsum(data, εs; kwargs...)
    return linear_region(log.(εs), log.(cm))[2]
end

################################################################################
# Boxed Correlation sum (we distribute data to boxes beforehand)
################################################################################
"""
    boxed_correlationsum(X::Dataset, εs, r0 = maximum(εs); kwargs...) → Cs

Estimate the box assisted q-order correlation sum[^Kantz2003] `Cs` out of a
dataset `X` for each radius in `εs`, by splitting the data into boxes of size `r0`
beforehand. This method is much faster than [`correlationsum`](@ref), **provided that** the
box size `r0` is significantly smaller than then the attractor length.
A good estimate for `r0` is [`estimate_r0_buenoorovio`](@ref).

    boxed_correlationsum(X; kwargs...) → εs, Cs

In this method the minimum inter-point distance and [`estimate_r0_buenoorovio`](@ref)
are used to estimate good `εs` for the calculation, which are also returned.

## Keywords
* `q = 2`
* `P = autoprismdim(data)` : The prism dimension.
* `w = 0` : The [Theiler window](@ref).
* `show_progress = false` : Whether to display a progress bar for the calculation.

## Description
`C_q(ε)` is calculated for every `ε ∈ εs` and each of the boxes to then be
summed up afterwards. The method of splitting the data into boxes was
implemented according to Theiler[^Theiler1987]. `w` is the [Theiler window](@ref).
`P` is the prism dimension. If `P` is unequal to the dimension of the data, only the
first `P` dimensions are considered for the box distribution (this is called the
prism-assisted version). By default `P` is choosen automatically.

The function is explicitly optimized for `q = 2` but becomes quite slow for `q ≠ 2`.

See [`correlationsum`](@ref) for the definition of `C_q`
and also [`data_boxing`](@ref) to use the algorithm that splits data into boxes.

[^Kantz]: Kantz, H., & Schreiber, T. (2003). [More about invariant quantities. In Nonlinear Time Series Analysis (pp. 197-233). Cambridge: Cambridge University Press.](https://doi:10.1017/CBO9780511755798.013)

[^Theiler1987]: Theiler, [Efficient algorithm for estimating the correlation dimension from a set of discrete points. Physical Review A, 36](https://doi.org/10.1103/PhysRevA.36.4456)
"""
function boxed_correlationsum(data; q = 2, P = autoprismdim(data), w = 0, show_progress=false)
    r0, ε0 = estimate_r0_buenoorovio(data, P)
    @assert  r0 < ε0 "The calculated box size was smaller than the minimum interpoint " *
    "distance. Please choose manually."
    εs = 10.0 .^ range(log10(ε0), log10(r0), length = 16)
    boxed_correlationsum(data, εs, r0; q, P, w, show_progress)
end

function boxed_correlationsum(
        data, εs, r0 = maximum(εs);
        q = 2, P = autoprismdim(data), w = 0,
        show_progress = false,
    )
    @assert P ≤ size(data, 2) "Prism dimension has to be lower or equal than " *
    "data dimension."
    boxes, contents = data_boxing(data, r0, P)
    if q == 2
        boxed_correlationsum_2(boxes, contents, data, εs; w, show_progress)
    else
        boxed_correlationsum_q(boxes, contents, data, εs, q; w, show_progress)
    end
end

"""
    autoprismdim(data, version = :bueno)

An algorithm to find the ideal choice of a prism dimension for [`boxed_correlationsum`](@ref).
`version = :bueno` uses `P=2`, while `version = :theiler` uses Theiler's original suggestion.
"""
function autoprismdim(data, version = :bueno)
    D = dimension(data)
    N = length(data)
    if version == :bueno
        return min(D, 2)
    elseif version == :theiler
        if D > 0.75 * log2(N)
            return max(2, ceil(0.5 * log2(N)))
        else
            return D
        end
    end
end

"""
    data_boxing(data, r0, P) → boxes, contents
Distribute the `data` points into boxes of size `r0`. Return box positions
and the contents of each box as two separate vectors. Implemented according to
the paper by Theiler[^Theiler1987] improving the algorithm by Grassberger and
Procaccia[^Grassberger1983]. If `P` is smaller than the dimension of the data,
only the first `P` dimensions are considered for the distribution into boxes.

See also: [`boxed_correlationsum`](@ref).

[^Theiler1987]: Theiler, [Efficient algorithm for estimating the correlation dimension from a set of discrete points. Physical Review A, 36](https://doi.org/10.1103/PhysRevA.36.4456)

[^Grassberger1983]: Grassberger and Proccacia, [Characterization of strange attractors, PRL 50 (1983)](https://journals-aps-org.e-bis.mpimet.mpg.de/prl/abstract/10.1103/PhysRevLett.50.346)
"""
function data_boxing(data, r0, P)
    @assert P ≤ size(data, 2) "Prism dimension has to be ≤ than data dimension."
    mini = minima(data)[1:P]

    # Map each datapoint to its bin edge and sort the resulting list:
    bins = map(point -> floor.(Int, (point[1:P] - mini)/r0), data)
    permutations = sortperm(bins, alg=QuickSort)

    boxes = unique(bins[permutations])
    contents = Vector{Vector{Int}}()
    sizehint!(contents, length(boxes))

    prior, prior_perm = 1, permutations[1]
    # distributes all permutation indices into boxes
    for (index, perm) in enumerate(permutations)
        if bins[perm] ≠ bins[prior_perm]
            push!(contents, permutations[prior:index-1])
            prior, prior_perm = index, perm
        end
    end
    push!(contents, permutations[prior:end])

    Dataset(boxes), contents
end

"""
    boxed_correlationsum_2(boxes, contents, data, εs; w = 0)
For a vector of `boxes` and the indices of their `contents` inside of `data`,
calculate the classic correlationsum of a radius or multiple radii `εs`.
`w` is the Theiler window, for explanation see [`boxed_correlationsum`](@ref).
"""
function boxed_correlationsum_2(boxes, contents, data, εs; w = 0, show_progress = false)
    Cs = zeros(eltype(data), length(εs))
    N = length(data)
    M = length(boxes)
    if show_progress
        progress = ProgressMeter.Progress(M; desc = "Boxed correlation sum: ", dt = 1.0)
    end
    for index in 1:M
        indices_neighbors = find_neighborboxes_2(index, boxes, contents)
        indices_box = contents[index]
        Cs .+= inner_correlationsum_2(indices_box, indices_neighbors, data, εs; w)
        show_progress && ProgressMeter.update!(progress, index)
    end
    Cs .* (2 / ((N - w) * (N - w - 1)))
end

"""
    find_neighborboxes_2(index, boxes, contents) → indices
For an `index` into `boxes` all neighbouring boxes beginning from the current
one are searched. If the found box is indeed a neighbour, the `contents` of
that box are added to `indices`.
"""
function find_neighborboxes_2(index, boxes, contents)
    indices = Int[]
    box = boxes[index]
    N_box = length(boxes)
    for index2 in index:N_box
        if evaluate(Chebyshev(), box, boxes[index2]) < 2
            append!(indices, contents[index2])
        end
    end
    indices
end

"""
    inner_correlationsum_2(indices_X, indices_Y, data, εs; norm = Euclidean(), w = 0)
Calculates the classic correlation sum for values `X` inside a box, considering
`Y` consisting of all values in that box and the ones in neighbouring boxes for
all distances `ε ∈ εs` calculated by `norm`. To obtain the position of the
values in the original time series `data`, they are passed as `indices_X` and
`indices_Y`.

`w` is the Theiler window. Each index to the original array is checked for the
distance of the compared index. If this absolute value is not higher than `w`
its element is not used in the calculation of the correlationsum.

See also: [`correlationsum`](@ref)
"""
function inner_correlationsum_2(indices_X, indices_Y, data, εs; norm = Euclidean(), w = 0)
    @assert issorted(εs) "Sorted εs required for optimized version."
    Cs, Ny, Nε = zeros(length(εs)), length(indices_Y), length(εs)
    for (i, index_X) in enumerate(indices_X)
    	x = data[index_X]
        for j in i+1:Ny
            index_Y = indices_Y[j]
            # Check for Theiler window.
            if abs(index_Y - index_X) > w
                # Calculate distance.
		        dist = evaluate(norm, data[index_Y], x)
		        for k in Nε:-1:1
		            if dist < εs[k]
		                Cs[k] += 1
		            else
		                break
		            end
		        end
		    end
        end
    end
    return Cs
end

"""
    boxed_correlationsum_q(boxes, contents, data, εs, q; w = 0)
For a vector of `boxes` and the indices of their `contents` inside of `data`,
calculate the `q`-order correlationsum of a radius or radii `εs`.
`w` is the Theiler window, for explanation see [`boxed_correlationsum`](@ref).
"""
function boxed_correlationsum_q(boxes, contents, data, εs, q; w = 0, show_progress = false)
    q <= 1 && @warn "This function is currently not specialized for q <= 1" *
    " and may show unexpected behaviour for these values."
    Cs = zeros(eltype(data), length(εs))
    N = length(data)
    M = length(boxes)
    if show_progress
        progress = ProgressMeter.Progress(M; desc = "Boxed correlation sum: ", dt = 1.0)
    end
    for index in 1:M
        indices_neighbors = find_neighborboxes_q(index, boxes, contents, q)
        indices_box = contents[index]
        Cs .+= inner_correlationsum_q(indices_box, indices_neighbors, data, εs, q; w)
        show_progress && ProgressMeter.update!(progress, index)
    end
    clamp.((Cs ./ ((N - 2w) * (N - 2w - 1) ^ (q-1))), 0, Inf) .^ (1 / (q-1))
end

"""
    find_neighborboxes_q(index, boxes, contents, q) → indices
For an `index` into `boxes` all neighbouring boxes are searched. If the found
box is indeed a neighbour, the `contents` of that box are added to `indices`.
"""
function find_neighborboxes_q(index, boxes, contents, q)
    indices = Int[]
    box = boxes[index]
    for (index2, box2) in enumerate(boxes)
        if evaluate(Chebyshev(), box, box2) < 2
            append!(indices, contents[index2])
        end
    end
    indices
end

"""
    inner_correlationsum_q(indices_X, indices_Y, data, εs, q::Real; norm = Euclidean(), w = 0)
Calculates the `q`-order correlation sum for values `X` inside a box,
considering `Y` consisting of all values in that box and the ones in
neighbouring boxes for all distances `ε ∈ εs` calculated by `norm`. To obtain
the position of the values in the original time series `data`, they are passed
as `indices_X` and `indices_Y`.

`w` is the Theiler window. The first and last `w` points of this data set are
not used by themselves to calculate the correlationsum.

See also: [`correlationsum`](@ref)
"""
function inner_correlationsum_q(indices_X, indices_Y, data, εs, q::Real; norm = Euclidean(), w = 0)
    @assert issorted(εs) "Sorted εs required for optimized version."
    Cs = zeros(length(εs))
    N, Ny, Nε = length(data), length(indices_Y), length(εs)
    for i in indices_X
        # Check that this index is not within Theiler window of the boundary
        # This step is neccessary for easy normalisation.
        (i < w + 1 || i > N - w) && continue
        C_current = zeros(Nε)
        x = data[i]
        for j in indices_Y
            # Check that this index is not whithin the Theiler window
        	if abs(i - j) > w
                # Calculate the distance for the correlationsum
		        dist = evaluate(norm, x, data[j])
		        for k in Nε:-1:1
		            if dist < εs[k]
		                C_current[k] += 1
		            else
		                break
		            end
		        end
		    end
        end
        Cs .+= C_current .^ (q-1)
    end
    return Cs
end

#######################################################################################
# Good boxsize estimates for boxed correlation sum
#######################################################################################
"""
    estimate_r0_theiler(X::Dataset) → r0, ε0
Estimate a reasonable size for boxing the data `X` before calculating the
[`boxed_correlationsum`](@ref) proposed by Theiler[^Theiler1987].
Return the boxing size `r0` and minimum inter-point distance in `X`, `ε0`.

To do so the dimension is estimated by running the algorithm by Grassberger and
Procaccia[^Grassberger1983] with `√N` points where `N` is the number of total
data points. Then the optimal boxsize ``r_0`` computes as
```math
r_0 = R (2/N)^{1/\\nu}
```
where ``R`` is the size of the chaotic attractor and ``\\nu`` is the estimated dimension.

[^Theiler1987]: Theiler, [Efficient algorithm for estimating the correlation dimension from a set of discrete points. Physical Review A, 36](https://doi.org/10.1103/PhysRevA.36.4456)

[^Grassberger1983]: Grassberger and Proccacia, [Characterization of strange attractors, PRL 50 (1983)](https://journals-aps-org.e-bis.mpimet.mpg.de/prl/abstract/10.1103/PhysRevLett.50.346)
"""
function estimate_r0_theiler(data)
    N = length(data)
    mini, maxi = minmaxima(data)
    R = mean(maxi .- mini)
    # Sample √N datapoints for a rough estimate of the dimension.
    data_sample = data[unique(rand(1:N, ceil(Int, sqrt(N))))] |> Dataset
    # Define radii for the rough dimension estimate
    min_d, _ = minimum_pairwise_distance(data)
    if min_d == 0
        @warn(
        "Minimum distance in the dataset is zero! Probably because of having data "*
        "with low resolution, or duplicate data points. Setting to `d₊/1000` for now.")
        min_d = R/(10^3)
    end
    lower = log10(min_d)
    εs = 10 .^ range(lower, stop = log10(R), length = 12)
    # Actually estimate the dimension.
    cm = correlationsum(data_sample, εs)
    ν = linear_region(log.(εs), log.(cm), tol = 0.5, warning = false)[2]
    # The combination yields the optimal box size
    r0 = R * (2/N)^(1/ν)
    return r0, min_d
end

"""
    estimate_r0_buenoorovio(X::Dataset, P = autoprismdim(X))
Estimates a reasonable size for boxing the time series `X` proposed by
Bueno-Orovio and Pérez-García[^Bueno2007] before calculating the correlation
dimension as presented by Theiler[^Theiler1983]. If instead of boxes, prisms
are chosen everything stays the same but `P` is the dimension of the prism.
To do so the dimension `ν` is estimated by running the algorithm by Grassberger
and Procaccia[^Grassberger1983] with `√N` points where `N` is the number of
total data points.
An effective size `ℓ` of the attractor is calculated by boxing a small subset
of size `N/10` into boxes of sidelength `r_ℓ` and counting the number of filled
boxes `η_ℓ`.
```math
\\ell = r_\\ell \\eta_\\ell ^{1/\\nu}
```
The optimal number of filled boxes `η_opt` is calculated by minimising the number
of calculations.
```math
\\eta_\\textrm{opt} = N^{2/3}\\cdot \\frac{3^\\nu - 1}{3^P - 1}^{1/2}.
```
`P` is the dimension of the data or the number of edges on the prism that don't
span the whole dataset.

Then the optimal boxsize ``r_0`` computes as
```math
r_0 = \\ell / \\eta_\\textrm{opt}^{1/\\nu}.
```

[^Bueno2007]: Bueno-Orovio and Pérez-García, [Enhanced box and prism assisted algorithms for computing the correlation dimension. Chaos Solitons & Fractrals, 34(5)](https://doi.org/10.1016/j.chaos.2006.03.043)

[^Theiler1987]: Theiler, [Efficient algorithm for estimating the correlation dimension from a set of discrete points. Physical Review A, 36](https://doi.org/10.1103/PhysRevA.36.4456)

[^Grassberger1983]: Grassberger and Proccacia, [Characterization of strange attractors, PRL 50 (1983)](https://journals-aps-org.e-bis.mpimet.mpg.de/prl/abstract/10.1103/PhysRevLett.50.346)
"""
function estimate_r0_buenoorovio(X, P = autoprismdim(X))
    mini, maxi = minmaxima(X)
    N = length(X)
    R = mean(maxi .- mini)
    # The possibility of a bad pick exists, if so, the calculation is repeated.
    ν = zero(eltype(X))
    min_d, _ = minimum_pairwise_distance(X)
    if min_d == 0
        @warn(
        "Minimum distance in the dataset is zero! Probably because of having data "*
        "with low resolution, or duplicate data points. Setting to `d₊/1000` for now.")
        min_d = R/(10^3)
    end

    # Sample N/10 datapoints out of data for rough estimate of effective size.
    sample1 = X[unique(rand(1:N, N÷10))] |> Dataset
    r_ℓ = R / 10
    η_ℓ = length(data_boxing(sample1, r_ℓ, P)[1])
    r0 = zero(eltype(X))
    while true
        # Sample √N datapoints for rough dimension estimate
        sample2 = X[unique(rand(1:N, ceil(Int, sqrt(N))))] |> Dataset
        # Define logarithmic series of radii.
        εs = 10.0 .^ range(log10(min_d), log10(R); length = 16)
        # Estimate ν from a sample using the Grassberger Procaccia algorithm.
        cm = correlationsum(sample2, εs)
        ν = linear_region(log.(εs), log.(cm); tol = 0.5, warning = false)[2]
        # Estimate the effictive size of the chaotic attractor.
        ℓ = r_ℓ * η_ℓ^(1/ν)
        # Calculate the optimal number of filled boxes according to Bueno-Orovio
        η_opt = N^(2/3) * ((3^ν - 1/2) / (3^P - 1))^(1/2)
        # The optimal box size is the effictive size divided by the box number # to the power of the inverse dimension.
        r0 = ℓ / η_opt^(1/ν)
        !isnan(r0) && break
    end
    return r0, min_d
end


text:
using Base.Iterators: repeated

test = """
    L.LL.LL.LL
    LLLLLLL.LL
    L.L.L..L..
    LLLL.LL.LL
    L.LL.LL.LL
    L.LLLLL.LL
    ..L.L.....
    LLLLLLLLLL
    L.LLLLLL.L
    L.LLLLL.LL
    """

is_floor(spot) = spot == '.'
is_empty(spot) = spot == 'L'
is_occupied(spot) = spot == '#'

function seats(string::AbstractString)::Matrix
    rows = split(string, r"\s+"; keepempty = false)
    seats = map(only, hcat(split.(rows, "")...)) |> permutedims
end

function neighbors(seats::M, ci::CartesianIndex)::M where {M <: AbstractMatrix}
    return neighbors(seats, ci.I...)
end

function neighbors(seats::M, r::Int, c::Int)::M where {M <: AbstractMatrix}
    return seats[(r - 1:r + 1) ∩ (1:end),
                 (c - 1:c + 1) ∩ (1:end)]
end

function visible(seats::AbstractMatrix, ci::CartesianIndex)::Vector
    return visible(seats, ci.I...)
end

function visible(seats::AbstractMatrix, r::Int, c::Int)::Vector
    (R, C) = size(seats)
    (cs, rs) = repeated.((c, r))
    rays = [(r - 1:-1:1, cs), (r + 1:R, cs), # vertical
            (rs, c - 1:-1:1), (rs, c + 1:C), # horizontal
            (r - 1:-1:1, c - 1:-1:1), (r + 1:R, c + 1:C), # diagonal
            (r - 1:-1:1, c + 1:C), (r + 1:R, c - 1:-1:1)]
    return map(rays) do ray
        ray = seats[map(CartesianIndex, ray...)]
        i = findfirst(!is_floor, ray)
        return i === nothing ? '.' : ray[i]
    end
end

function step(seats::M, look_at::Function)::M where {M <: AbstractMatrix}
    return map(keys(seats)) do ci
        n_occ = count(is_occupied, look_at(seats, ci))
        return seats[ci] == 'L' && n_occ == 0 ? '#' :
            seats[ci] == '#' && n_occ >= 5 ? 'L' :
            seats[ci]
    end
end

function run(seats::M, look_at::Function)::M where {M <: AbstractMatrix}
    prev_seats = nothing
    while prev_seats != seats
        prev_seats = seats
        seats = step(seats, look_at)
    end
    return seats
end

function q1(seats::AbstractMatrix)::Integer
    final = run(seats, neighbors)
    return count(is_occupied, final)
end

function q2(seats::AbstractMatrix)::Integer
    final = run(seats, visible)
    return count(is_occupied, final)
end

let lounge = seats(test)
    @assert q1(lounge) == 37
    @assert q2(lounge) == 26
end

lounge = read("day11.in", String) |> seats
println("Q1: ", q1(lounge))
println("Q2: ", q2(lounge))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using ParamLevelSet
using LinearAlgebra
using Test
using jInv.Mesh

plotting = false;
if plotting
	using jInvVisPyPlot
	using PyPlot
	close("all");
end

function getDiamondModel(n::Array{Int64})
n_tup = tuple(n...);
n = collect(n_tup);
u = zeros(Float32,n_tup);
X,Y,Z = meshgrid(-div(n[1],2):div(n[1],2),-div(n[2],2):div(n[2],2),-div(n[3],2):div(n[3],2))
I = (abs.(X)+abs.(Y)+abs.(Z)) .<= div(n[1],4);
u[I] .= 1.0;
return u;
end

println("Rotation operator transpose test");
n = [64,64,64].+1;

u = getDiamondModel(n);
u = convert(Array{Float64}, u)

if plotting
	figure()
	plotModel(u)
end

theta_phi = deg2rad.([37.0;24.0]);
display(theta_phi)
b = [1.0;2.0;3.0];

# RT = generateSamplingMatrix(u,theta_phi,1,1)
# v = RT'*u[:]
# println("The norm is:")
# println(norm(RT*v - u[:]))

v, = rotateAndMove3D(u,theta_phi,b,false);

if plotting
	figure()
	plotModel(v)
end
#uu = rotateAndMove3DTranspose(v,theta_phi,b);
uu, = rotateAndMove3D(v,theta_phi,b,true);
if plotting
	figure()
	plotModel(uu)
end

@test abs(dot(v,v) - dot(uu,u)) < 1e-5
@test norm(uu[:]-u[:])/norm(u[:]) < 1.0


text:
module TestSymbolic

using Test
using VanFoFy.Symbolic: VarLinForm, variables, add!, add_conjugated!, mul!

function test()
    @testset "VarLinForm" begin
        form1 = VarLinForm{ComplexF64}()
        @test length(variables(form1)) == 0

        @test form1[5] == 0.0im

        form1[5] = 2.0+3.0im
        @test form1[5] == 2.0+3.0im

        form2 = VarLinForm{ComplexF64}()
        form2[3] = 1.0-8.0im
        add!(form1, form2, 2)
        @test form1[5] == 2.0+3.0im
        @test form1[3] == 2.0 - 16.0im

        add_conjugated!(form1, form2, 1.5)
        @test form1[5] == 2.0+3.0im
        @test form1[3] == 3.5 - 4.0im

        add!(form1, form1)
        @test form1[5] == 4.0+6.0im
        @test form1[3] == 7.0 - 8.0im

        add_conjugated!(form1, form1)
        @test form1[5] == 8.0+0.0im
        @test form1[3] == 14.0 - 0.0im

        form1[5] += 4.0im
        mul!(form1, 0.5)
        @test form1[5] == 4.0+2.0im
        @test form1[3] == 7.0 - 0.0im

        mul!(form1, 2)
        @test form1[5] == 8.0+4.0im
        @test form1[3] == 14.0 - 0.0im

        mul!(form1, 1im)
        @test form1[5] == 8.0im-4.0
        @test form1[3] == 14.0im - 0.0
    end
end
end # module TestSymbolic
using .TestSymbolic
TestSymbolic.test()


text:
baremodule MultiFloatsCodeGen

######################################################################## EXPORTS

export one_pass_renorm_func, two_pass_renorm_func,
    multifloat_add_func, multifloat_float_add_func,
    multifloat_mul_func, multifloat_float_mul_func,
    multifloat_div_func, multifloat_sqrt_func,
    multifloat_eq_func, multifloat_ne_func, multifloat_lt_func,
    multifloat_gt_func, multifloat_le_func, multifloat_ge_func

######################################################################## IMPORTS

using Base: +, -, *, div, !, (:), ==, !=, <, <=, >, >=, min, max,
    Vector, isempty, length, push!, deleteat!, reverse!, Dict, haskey, @assert

###################################################### METAPROGRAMMING UTILITIES

const SymExpr = Union{Symbol,Expr}

meta_tuple(x...) = Expr(:tuple, x...)

inline_block() = [Expr(:meta, :inline)] # same as @inline

function_def(name, args, body) =
    Expr(:function,
        Expr(:where, Expr(:call, name, args...), :T),
        Expr(:block, body...))

function_def_typed(name, arg_type, args, body) =
    function_def(name, [Expr(:(::), arg, arg_type) for arg in args], body)

meta_two_sum(s::Symbol, e::Symbol, a::SymExpr, b::SymExpr) =
    Expr(:(=), meta_tuple(s, e), Expr(:call, :two_sum, a, b))

meta_two_prod(p::Symbol, e::Symbol, a::SymExpr, b::SymExpr) =
    Expr(:(=), meta_tuple(p, e), Expr(:call, :two_prod, a, b))

meta_quick_two_sum(s::Symbol, e::Symbol, a::Symbol, b::Symbol) =
    Expr(:(=), meta_tuple(s, e), Expr(:call, :quick_two_sum, a, b))

meta_sum(addends::Vector{Symbol}) =
    length(addends) == 1 ? addends[1] : Expr(:call, :+, addends...)

meta_sum(result::Symbol, addends::Vector{Symbol}) =
    Expr(:(=), result, meta_sum(addends))

meta_multifloat(N::Int) = :(MultiFloat{T,$N})

renorm_name(n::Int) = Symbol("renorm_", n)

mpadd_name(src_len::Int, dst_len::Int) = Symbol("mpadd_", src_len, '_', dst_len)

################################################################ RENORMALIZATION

function one_pass_renorm_func(n::Int; sloppy::Bool=false)
    args = [Symbol('a', i) for i = 0 : n - sloppy]
    sums = [Symbol('s', i) for i = 0 : n - 1]
    if (n == 2) && sloppy
        return function_def_typed(renorm_name(2), :T, args,
            [Expr(:call, :quick_two_sum, args...)])
    end
    code = inline_block()
    push!(code, meta_quick_two_sum(sums[1], sums[2], args[1], args[2]))
    for i = 1 : n-2
        push!(code, meta_quick_two_sum(
            sums[i+1], sums[i+2], sums[i+1], args[i+2]))
    end
    push!(code, sloppy ? meta_tuple(sums...) :
        meta_tuple(sums[1:n-1]..., meta_sum([sums[n], args[n+1]])))
    function_def_typed(renorm_name(n), :T, args, code)
end

function two_pass_renorm_func(n::Int; sloppy::Bool=false)
    args = [Symbol('a', i) for i = 0 : n - sloppy]
    sums = [Symbol('s', i) for i = 0 : n - 1]
    if (n == 2) && sloppy
        return function_def_typed(renorm_name(2), :T, args,
            [Expr(:call, :quick_two_sum, args...)])
    end
    code = inline_block()
    push!(code, meta_quick_two_sum(:t, args[end], args[end-1], args[end]))
    for i = n-2-sloppy : -1 : 1
        push!(code, meta_quick_two_sum(:t, args[i+2], args[i+1], :t))
    end
    push!(code, meta_quick_two_sum(sums[1], sums[2], args[1], :t))
    for i = 1 : n-2
        push!(code, meta_quick_two_sum(
            sums[i+1], sums[i+2], sums[i+1], args[i+2]))
    end
    push!(code, sloppy ? meta_tuple(sums...) :
        meta_tuple(sums[1:n-1]..., meta_sum([sums[n], args[n+1]])))
    function_def_typed(renorm_name(n), :T, args, code)
end

################################################################### ACCUMULATION

function add_var!(vars::Vector{Tuple{Symbol,Int}}, prefix::Char, i::Int)
    var = Symbol(prefix, i)
    push!(vars, (var, i))
    var
end

function add_var!(vars::Vector{Tuple{Symbol,Int}}, prefix::Char, i::Int, j::Int)
    var = Symbol(prefix, i, '_', j)
    push!(vars, (var, j))
    var
end

function mpadd_func(src_len::Int, dst_len::Int)
    @assert src_len >= dst_len
    args = [Symbol('a', i) for i = 1 : src_len]
    sums = [Symbol('s', i) for i = 0 : dst_len - 1]
    vars = [(a, 0) for a in args]
    code = inline_block()
    k = 0
    for i = 0 : dst_len-2
        while true
            addends = [v[1] for v in vars if v[2] == i]
            if length(addends) == 1
                push!(code, Expr(:(=), sums[i+1], addends[1]))
                deleteat!(vars, [v[1] == addends[1] for v in vars])
                break
            end
            for j = 1 : div(length(addends), 2)
                a, b = addends[2*j-1], addends[2*j]
                push!(code, meta_two_sum(add_var!(vars, 'm', k, i),
                    add_var!(vars, 'm', k+1, i+1), a, b))
                k += 2
                deleteat!(vars, [v[1] == a || v[1] == b for v in vars])
            end
        end
    end
    push!(code, meta_sum(sums[dst_len], [v[1] for v in vars]))
    push!(code, meta_tuple(sums...))
    function_def_typed(mpadd_name(src_len, dst_len), :T, args, code)
end

const MPADD_CACHE = Dict{Symbol,Expr}()

function _mpsum(results::Vector{Symbol}, addends::Vector{Symbol})
    src_len = length(addends)
    dst_len = length(results)
    @assert src_len >= dst_len
    if dst_len == 1
        meta_sum(results[1], addends)
    elseif src_len == dst_len == 2
        meta_two_sum(results[1], results[2], addends[1], addends[2])
    else
        func_name = mpadd_name(src_len, dst_len)
        if !haskey(MPADD_CACHE, func_name)
            MPADD_CACHE[func_name] = mpadd_func(src_len, dst_len)
        end
        Expr(:(=), meta_tuple(results...), Expr(:call, func_name, addends...))
    end
end

function generate_accumulation_code!(code::Vector{Expr},
        vars::Vector{Tuple{Symbol,Int}}, N::Int; sloppy::Bool=false)
    sums = [Symbol('s', i) for i = 0 : N - sloppy]
    for i = 0 : N - sloppy
        addends = [v[1] for v in vars if v[2] == i]
        results = [sums[i+1]]
        for j = i + 1 : min(i + length(addends) - 1, N - sloppy)
            push!(results, add_var!(vars, 'm', i, j))
        end
        push!(code, _mpsum(results, addends))
        deleteat!(vars, [v[2] <= i for v in vars])
    end
    push!(code, Expr(:call, meta_multifloat(N),
        Expr(:call, renorm_name(N), sums...)))
end

########################################################### ARITHMETIC FUNCTIONS

function multifloat_add_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    vars = Tuple{Symbol,Int}[]
    for i = 1 : N - sloppy
        tmp = add_var!(vars, 't', i-1)
        err = add_var!(vars, 'e', i)
        push!(code, meta_two_sum(tmp, err, :(a._limbs[$i]), :(b._limbs[$i])))
    end
    if sloppy
        tmp = add_var!(vars, 't', N-1)
        push!(code, :($tmp = a._limbs[$N] + b._limbs[$N]))
    end
    reverse!(vars)
    generate_accumulation_code!(code, vars, N, sloppy=sloppy)
    function_def_typed(:multifloat_add, meta_multifloat(N), [:a, :b], code)
end

function multifloat_float_add_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    vars = Tuple{Symbol,Int}[]
    push!(code, meta_two_sum(add_var!(vars, 't', 0), add_var!(vars, 'e', 1),
        :(a._limbs[1]), :b))
    for i = 2 : N
        push!(code, Expr(:(=), add_var!(vars, 't', i-1), :(a._limbs[$i])))
    end
    reverse!(vars)
    generate_accumulation_code!(code, vars, N, sloppy=sloppy)
    function_def(:multifloat_float_add,
        [Expr(:(::), :a, meta_multifloat(N)), Expr(:(::), :b, :T)], code)
end

function multifloat_mul_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    for i = 0 : N-1-sloppy, j = 0 : i
        push!(code, meta_two_prod(
            Symbol('t', j, '_', i), Symbol('e', j, '_', i+1),
            :(a._limbs[$(j+1)]), :(b._limbs[$(i-j+1)])))
    end
    for i = 0 : N-2+sloppy
        push!(code, Expr(:(=), Symbol('t', i, '_', N - sloppy),
            :(a._limbs[$(i + 2 - sloppy)] * b._limbs[$(N - i)])))
    end
    vars = Tuple{Symbol,Int}[]
    for i = 0 : N-1-sloppy, j = 0 : i; add_var!(vars, 't', j, i       ); end
    for i = 0 : N-2+sloppy           ; add_var!(vars, 't', i, N-sloppy); end
    for i = 0 : N-1-sloppy, j = 0 : i; add_var!(vars, 'e', j, i+1     ); end
    generate_accumulation_code!(code, vars, N, sloppy=sloppy)
    function_def_typed(:multifloat_mul, meta_multifloat(N), [:a, :b], code)
end

function multifloat_float_mul_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    for i = 0 : N-1-sloppy
        push!(code, meta_two_prod(Symbol('t', i), Symbol('e', i+1),
            :(a._limbs[$(i+1)]), :b))
    end
    if sloppy
        push!(code, Expr(:(=), Symbol('t', N-1), :(a._limbs[$N] * b)))
    end
    vars = Tuple{Symbol,Int}[]
    for i = 0 : N-1;      add_var!(vars, 't', i); end
    for i = 1 : N-sloppy; add_var!(vars, 'e', i); end
    generate_accumulation_code!(code, vars, N, sloppy=sloppy)
    function_def(:multifloat_float_mul,
        [Expr(:(::), :a, meta_multifloat(N)), Expr(:(::), :b, :T)], code)
end

function multifloat_div_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    quots = [Symbol('q', i) for i = 0 : N - sloppy]
    push!(code, :($(quots[1]) = a._limbs[1] / b._limbs[1]))
    push!(code, :(r = a - b * $(quots[1])))
    for i = 2 : N-sloppy
        push!(code, :($(quots[i]) = r._limbs[1] / b._limbs[1]))
        push!(code, :(r -= b * $(quots[i])))
    end
    push!(code, :($(quots[N+1-sloppy]) = r._limbs[1] / b._limbs[1]))
    push!(code, Expr(:call, meta_multifloat(N),
        Expr(:call, renorm_name(N), quots...)))
    function_def_typed(:multifloat_div, meta_multifloat(N), [:a, :b], code)
end

num_sqrt_iters(N::Int, sloppy::Bool) =
    (2 <= N <= 3) ? 2 :
    (4 <= N <= 5) ? 3 :
    sloppy ? div(N + 1, 2) : div(N + 2, 2)

function multifloat_sqrt_func(N::Int; sloppy::Bool=false)
    code = inline_block()
    push!(code, :(r = MultiFloat{T,$N}(inv(unsafe_sqrt(x._limbs[1])))))
    push!(code, :(h = scale(T(0.5), x)))
    for _ = 1 : num_sqrt_iters(N, sloppy)
        push!(code, :(r += r * (T(0.5) - h * (r * r))))
    end
    push!(code, :(r * x))
    function_def_typed(:multifloat_sqrt, meta_multifloat(N), [:x], code)
end

########################################################### RECURSIVE COMPARISON

eq_expr(n::Int) = (n == 1
    ? :(x._limbs[$n] == y._limbs[$n])
    : :($(eq_expr(n-1)) & (x._limbs[$n] == y._limbs[$n])))

ne_expr(n::Int) = (n == 1
    ? :(x._limbs[$n] != y._limbs[$n])
    : :($(ne_expr(n-1)) | (x._limbs[$n] != y._limbs[$n])))

lt_expr(m::Int, n::Int) = (m == n
    ? :(x._limbs[$m] < y._limbs[$m])
    : :((x._limbs[$m] < y._limbs[$m]) |
        ((x._limbs[$m] == y._limbs[$m]) & $(lt_expr(m+1, n)))))

gt_expr(m::Int, n::Int) = (m == n
    ? :(x._limbs[$m] > y._limbs[$m])
    : :((x._limbs[$m] > y._limbs[$m]) |
        ((x._limbs[$m] == y._limbs[$m]) & $(gt_expr(m+1, n)))))

le_expr(m::Int, n::Int) = (m == n
    ? :(x._limbs[$m] <= y._limbs[$m])
    : :((x._limbs[$m] < y._limbs[$m]) |
        ((x._limbs[$m] == y._limbs[$m]) & $(le_expr(m+1, n)))))

ge_expr(m::Int, n::Int) = (m == n
    ? :(x._limbs[$m] >= y._limbs[$m])
    : :((x._limbs[$m] > y._limbs[$m]) |
        ((x._limbs[$m] == y._limbs[$m]) & $(ge_expr(m+1, n)))))

########################################################### COMPARISON FUNCTIONS

cmp_func(name, N, expr) = function_def_typed(name,
    meta_multifloat(N), [:x, :y], push!(inline_block(), expr))

multifloat_eq_func(N::Int) = cmp_func(:multifloat_eq, N, eq_expr(N))
multifloat_ne_func(N::Int) = cmp_func(:multifloat_ne, N, ne_expr(N))
multifloat_lt_func(N::Int) = cmp_func(:multifloat_lt, N, lt_expr(1, N))
multifloat_gt_func(N::Int) = cmp_func(:multifloat_gt, N, gt_expr(1, N))
multifloat_le_func(N::Int) = cmp_func(:multifloat_le, N, le_expr(1, N))
multifloat_ge_func(N::Int) = cmp_func(:multifloat_ge, N, ge_expr(1, N))

end # baremodule MultiFloatsCodeGen


text:
@testset "Testing the Dynamic" begin

# We start by testing the coordinate change
import ResponseManneville.InducedLSVMapDefinition: CoordinateChange, InvCoordinateChange

@test CoordinateChange(0.5) == 0
@test CoordinateChange(1) == 1
@test InvCoordinateChange(0) == 0.5
@test CoordinateChange(1) == 1

# We start now the map
import ResponseManneville.InducedLSVMapDefinition: ApproxInducedLSV, derleft, derderleft
D = ApproxInducedLSV(0.5, 10)
@test D.α == 0.5
@test derleft(D, 0) == 1
@test derderleft(D, 0) == Inf

using TaylorSeries, IntervalArithmetic

T(x) = ResponseManneville.DynamicDefinition.plottable(D, x)
x = Taylor1([0.6, 1], 2)
y = T(x)

# we test the Taylor Series, the last branch is linear with derivative 2
@test y[1] == 2
@test y[2] == 0

# we test now the preimwithder_derder function
test_x = Base.rand(Float64, 10)

for t in test_x
    for k in 1:ResponseManneville.DynamicDefinition.nbranches(D)
        # preimwithder_derder computes the preimage through the k-th branch
        # with tight enclosures for the derivative and the second derivative
        val = ResponseManneville.InducedLSVMapDefinition.preimwithder_derder(D, k, t, 0.000001)
        val_z = T(Taylor1([mid(val[1]), 1], 2)) #this is computed through TaylorSeries
        @test abs(val_z[0]-t)<2^-15 && abs(val_z[1]-mid(Interval(val[2])))<2^-15 && abs(2*val_z[2]-mid(Interval(val[3])))<2^-15
    end
end

end;

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
"""
    artefact_description(x)
    artefact_description!(x, val)

Freeform description of the observed subject artefact and its possible cause (e.g.,
"door open", "nurse walked into room at 2 min", "seizure at 10 min"). If this
field is left empty, it will be interpreted as absence of artifacts.
"""
@defprop ArtefactDescription{:artefact_description}::String

"""
    intended_for(x)
    intended_for!(x, val)

Path or list of path relative to the subject subfolder pointing to the structural
MRI, possibly of different types if a list is specified, to be used with the MEG
recording. The path(s) need(s) to use forward slashes instead of backward slashes
(e.g. ses-/anat/sub-01_T1w.nii.gz).
"""
@defprop IntendedFor{:intended_for}::Vector{String}

"""
    instructions(x)
    instructions!(x, val)

Text of the instructions given to participants before the scan. This is especially
important in context of resting state fMRI and distinguishing between eyes open and
eyes closed paradigms.
"""
@defprop Instructions{:instructions}::String

"""
    cog_atlas_id(x)
    cog_atlas_id!(x, val)

URL of the corresponding [Cognitive Atlas Task](https://www.cognitiveatlas.org/) term.
"""
@defprop CogAtlasID{:cog_atlas_id}::String

"""
    cog_poid(x)
    cog_poid!(x, val)

URL of the corresponding [CogPO](http://www.cogpo.org/) term.
"""
@defprop CogPOID{:cog_poid}::String