text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

# Functions
for T in (:Taylor1, :TaylorN)
    @eval begin
        function exp(a::$T)
            order = a.order
            aux = exp(constant_term(a))
            aa = one(aux) * a
            c = $T( aux, order )
            for k in eachindex(a)
                exp!(c, aa, k)
            end
            return c
        end

        function log(a::$T)
            iszero(constant_term(a)) && throw(DomainError(a,
                    """The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))

            order = a.order
            aux = log(constant_term(a))
            aa = one(aux) * a
            c = $T( aux, order )
            for k in eachindex(a)
                log!(c, aa, k)
            end
            return c
        end

        sin(a::$T) = sincos(a)[1]
        cos(a::$T) = sincos(a)[2]
        function sincos(a::$T)
            order = a.order
            aux = sin(constant_term(a))
            aa = one(aux) * a
            s = $T( aux, order )
            c = $T( cos(constant_term(a)), order )
            for k in eachindex(a)
                sincos!(s, c, aa, k)
            end
            return s, c
        end

        function tan(a::$T)
            order = a.order
            aux = tan(constant_term(a))
            aa = one(aux) * a
            c = $T(aux, order)
            c2 = $T(aux^2, order)
            for k in eachindex(a)
                tan!(c, aa, c2, k)
            end
            return c
        end

        function asin(a::$T)
            a0 = constant_term(a)
            a0^2 == one(a0) && throw(DomainError(a,
                """Series expansion of asin(x) diverges at x = ±1."""))

            order = a.order
            aux = asin(a0)
            aa = one(aux) * a
            c = $T( aux, order )
            r = $T( sqrt(1 - a0^2), order )
            for k in eachindex(a)
                asin!(c, aa, r, k)
            end
            return c
        end

        function acos(a::$T)
            a0 = constant_term(a)
            a0^2 == one(a0) && throw(DomainError(a,
            """Series expansion of asin(x) diverges at x = ±1."""))

            order = a.order
            aux = acos(a0)
            aa = one(aux) * a
            c = $T( aux, order )
            r = $T( sqrt(1 - a0^2), order )
            for k in eachindex(a)
                acos!(c, aa, r, k)
            end
            return c
        end

        function atan(a::$T)
            order = a.order
            a0 = constant_term(a)
            aux = atan(a0)
            aa = one(aux) * a
            c = $T( aux, order)
            r = $T(one(aux) + a0^2, order)
            iszero(constant_term(r)) && throw(DomainError(a,
                """Series expansion of atan(x) diverges at x = ±im."""))

            for k in eachindex(a)
                atan!(c, aa, r, k)
            end
            return c
        end

        function atan(a::$T, b::$T)
            c = atan(a/b)
            c[0] = atan(constant_term(a), constant_term(b))
            return c
        end

        sinh(a::$T) = sinhcosh(a)[1]
        cosh(a::$T) = sinhcosh(a)[2]
        function sinhcosh(a::$T)
            order = a.order
            aux = sinh(constant_term(a))
            aa = one(aux) * a
            s = $T( aux, order)
            c = $T( cosh(constant_term(a)), order)
            for k in eachindex(a)
                sinhcosh!(s, c, aa, k)
            end
            return s, c
        end

        function tanh(a::$T)
            order = a.order
            aux = tanh( constant_term(a) )
            aa = one(aux) * a
            c = $T( aux, order)
            c2 = $T( aux^2, order)
            for k in eachindex(a)
                tanh!(c, aa, c2, k)
            end
            return c
        end

        function asinh(a::$T)
            order = a.order
            a0 = constant_term(a)
            aux = asinh(a0)
            aa = one(aux) * a
            c = $T( aux, order )
            r = $T( sqrt(a0^2 + 1), order )
            iszero(constant_term(r)) && throw(DomainError(a,
                """Series expansion of asinh(x) diverges at x = ±im."""))
            for k in eachindex(a)
                asinh!(c, aa, r, k)
            end
            return c
        end

        function acosh(a::$T)
            a0 = constant_term(a)
            a0^2 == one(a0) && throw(DomainError(a,
                """Series expansion of acosh(x) diverges at x = ±1."""))

            order = a.order
            aux = acosh(a0)
            aa = one(aux) * a
            c = $T( aux, order )
            r = $T( sqrt(a0^2 - 1), order )
            for k in eachindex(a)
                acosh!(c, aa, r, k)
            end
            return c
        end

        function atanh(a::$T)
            order = a.order
            a0 = constant_term(a)
            aux = atanh(a0)
            aa = one(aux) * a
            c = $T( aux, order)
            r = $T(one(aux) - a0^2, order)
            iszero(constant_term(r)) && throw(DomainError(a,
                """Series expansion of atanh(x) diverges at x = ±1."""))

            for k in eachindex(a)
                atanh!(c, aa, r, k)
            end
            return c
        end

    end
end

# Recursive functions (homogeneous coefficients)
for T in (:Taylor1, :TaylorN)
    @eval begin
        @inline function identity!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            @inbounds c[k] = identity(a[k])
            return nothing
        end

        @inline function zero!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            @inbounds c[k] = zero(a[k])
            return nothing
        end

        @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds c[0] = one(a[0])
            else
                @inbounds c[k] = zero(a[k])
            end
            return nothing
        end

        @inline function abs!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            z = zero(constant_term(a))
            if constant_term(constant_term(a)) > constant_term(z)
                return add!(c, a, k)
            elseif constant_term(constant_term(a)) < constant_term(z)
                return subst!(c, a, k)
            else
                throw(DomainError(a,
                """The 0th order coefficient must be non-zero
                (abs(x) is not differentiable at x=0)."""))
            end
            return nothing
        end

        @inline abs2!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} = sqr!(c, a, k)

        @inline function exp!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds c[0] = exp(constant_term(a))
                return nothing
            end
            if $T == Taylor1
                @inbounds c[k] = k * a[k] * c[0]
            else
                @inbounds mul!(c[k], k * a[k], c[0])
            end
            @inbounds for i = 1:k-1
                if $T == Taylor1
                    c[k] += (k-i) * a[k-i] * c[i]
                else
                    mul!(c[k], (k-i) * a[k-i], c[i])
                end
            end
            @inbounds c[k] = c[k] / k
            return nothing
        end

        @inline function log!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds c[0] = log(constant_term(a))
                return nothing
            elseif k == 1
                @inbounds c[1] = a[1] / constant_term(a)
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * a[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1)*a[1], c[k-1])
            end
            @inbounds for i = 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * a[i] * c[k-i]
                else
                    mul!(c[k], (k-i)*a[i], c[k-i])
                end
            end
            @inbounds c[k] = (a[k] - c[k]/k) / constant_term(a)
            return nothing
        end

        @inline function sincos!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds s[0], c[0] = sin( a0 ), cos( a0 )
                return nothing
            end

            x = a[1]
            if $T == Taylor1
                @inbounds s[k] = x * c[k-1]
                @inbounds c[k] = -x * s[k-1]
            else
                mul!(s[k], x, c[k-1])
                mul!(c[k], -x, s[k-1])
            end
            @inbounds for i = 2:k
                x = i * a[i]
                if $T == Taylor1
                    s[k] += x * c[k-i]
                    c[k] -= x * s[k-i]
                else
                    mul!(s[k], x, c[k-i])
                    mul!(c[k], -x, s[k-i])
                end
            end

            @inbounds s[k] = s[k] / k
            @inbounds c[k] = c[k] / k
            return nothing
        end

        @inline function tan!(c::$T{T}, a::$T{T}, c2::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds aux = tan( constant_term(a) )
                @inbounds c[0] = aux
                @inbounds c2[0] = aux^2
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * a[k-1] * c2[1]
            else
                @inbounds mul!(c[k], (k-1) * a[k-1], c2[1])
            end
            @inbounds for i = 1:k-1
                if $T == Taylor1
                    c[k] += (k-i) * a[k-i] * c2[i]
                else
                    mul!(c[k], (k-i) * a[k-i], c2[i])
                end
            end
            @inbounds c[k] = a[k] + c[k]/k
            sqr!(c2, c, k)

            return nothing
        end

        @inline function asin!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = asin( a0 )
                @inbounds r[0] = sqrt( 1 - a0^2 )
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            sqrt!(r, 1-a^2, k)
            @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r)
            return nothing
        end

        @inline function acos!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = acos( a0 )
                @inbounds r[0] = sqrt( 1 - a0^2 )
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            sqrt!(r, 1-a^2, k)
            @inbounds c[k] = -(a[k] + c[k]/k) / constant_term(r)
            return nothing
        end

        @inline function atan!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = atan( a0 )
                @inbounds r[0] = 1 + a0^2
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            @inbounds sqr!(r, a, k)
            @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r)
            return nothing
        end

        @inline function sinhcosh!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds s[0] = sinh( constant_term(a) )
                @inbounds c[0] = cosh( constant_term(a) )
                return nothing
            end

            x = a[1]
            if $T == Taylor1
                @inbounds s[k] = x * c[k-1]
                @inbounds c[k] = x * s[k-1]
            else
                @inbounds mul!(s[k], x, c[k-1])
                @inbounds mul!(c[k], x, s[k-1])
            end
            @inbounds for i = 2:k
                x = i * a[i]
                if $T == Taylor1
                    s[k] += x * c[k-i]
                    c[k] += x * s[k-i]
                else
                    mul!(s[k], x, c[k-i])
                    mul!(c[k], x, s[k-i])
                end
            end
            s[k] = s[k] / k
            c[k] = c[k] / k
            return nothing
        end

        @inline function tanh!(c::$T{T}, a::$T{T}, c2::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds aux = tanh( constant_term(a) )
                @inbounds c[0] = aux
                @inbounds c2[0] = aux^2
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = k * a[k] * c2[0]
            else
                @inbounds mul!(c[k], k * a[k], c2[0])
            end
            @inbounds for i = 1:k-1
                if $T == Taylor1
                    c[k] += (k-i) * a[k-i] * c2[i]
                else
                    mul!(c[k], (k-i) * a[k-i], c2[i])
                end
            end
            @inbounds c[k] = a[k] - c[k]/k
            sqr!(c2, c, k)

            return nothing
        end

        @inline function asinh!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = asinh( a0 )
                @inbounds r[0] = sqrt( a0^2 + 1 )
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            sqrt!(r, a^2+1, k)
            @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r)
            return nothing
        end

        @inline function acosh!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = acosh( a0 )
                @inbounds r[0] = sqrt( a0^2 - 1 )
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            sqrt!(r, a^2-1, k)
            @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r)
            return nothing
        end

        @inline function atanh!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T<:Number}
            if k == 0
                a0 = constant_term(a)
                @inbounds c[0] = atanh( a0 )
                @inbounds r[0] = 1 - a0^2
                return nothing
            end

            if $T == Taylor1
                @inbounds c[k] = (k-1) * r[1] * c[k-1]
            else
                @inbounds mul!(c[k], (k-1) * r[1], c[k-1])
            end
            @inbounds for i in 2:k-1
                if $T == Taylor1
                    c[k] += (k-i) * r[i] * c[k-i]
                else
                    mul!(c[k], (k-i) * r[i], c[k-i])
                end
            end
            @inbounds sqr!(r, a, k)
            @inbounds c[k] = (a[k] + c[k]/k) / constant_term(r)
            return nothing
        end

    end
end

@doc doc"""
    inverse(f)

Return the Taylor expansion of ``f^{-1}(t)``, of order `N = f.order`,
for `f::Taylor1` polynomial if the first coefficient of `f` is zero.
Otherwise, a `DomainError` is thrown.

The algorithm implements Lagrange inversion at ``t=0`` if ``f(0)=0``:
```math
\begin{equation*}
f^{-1}(t) = \sum_{n=1}^{N} \frac{t^n}{n!} \left.
    \frac{{\rm d}^{n-1}}{{\rm d} z^{n-1}}\left(\frac{z}{f(z)}\right)^n
    \right\vert_{z=0}.
\end{equation*}
```

""" inverse

function inverse(f::Taylor1{T}) where {T<:Number}
    if f[0] != zero(T)
        throw(DomainError(f,
        """
        Evaluation of Taylor1 series at 0 is non-zero. For high accuracy, revert
        a Taylor1 series with first coefficient 0 and re-expand about f(0).
        """))
    end
    z = Taylor1(T,f.order)
    zdivf = z/f
    zdivfpown = zdivf
    S = eltype(zdivf)
    coeffs = zeros(S,f.order+1)

    @inbounds for n in 1:f.order
        coeffs[n+1] = zdivfpown[n-1]/n
        zdivfpown *= zdivf
    end
    Taylor1(coeffs, f.order)
end

# Documentation for the recursion relations
@doc doc"""
    exp!(c, a, k) --> nothing

Update the `k-th` expansion coefficient `c[k+1]` of `c = exp(a)`
for both `c` and `a` either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{equation*}
c_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j.
\end{equation*}
```

""" exp!

@doc doc"""
    log!(c, a, k) --> nothing

Update the `k-th` expansion coefficient `c[k+1]` of `c = log(a)`
for both `c` and `a` either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{equation*}
c_k = \frac{1}{a_0} \big(a_k - \frac{1}{k} \sum_{j=0}^{k-1} j a_{k-j} c_j \big).
\end{equation*}
```

""" log!

@doc doc"""
    sincos!(s, c, a, k) --> nothing

Update the `k-th` expansion coefficients `s[k+1]` and `c[k+1]`
of `s = sin(a)` and `c = cos(a)` simultaneously, for `s`, `c` and `a`
either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{eqnarray*}
s_k &=&  \frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} c_j ,\\
c_k &=& -\frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} s_j.
\end{eqnarray*}
```

""" sincos!

@doc doc"""
    tan!(c, a, p, k::Int) --> nothing

Update the `k-th` expansion coefficients `c[k+1]` of `c = tan(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `p = c^2` and
is passed as an argument for efficiency.

The coefficients are given by

```math
\begin{equation*}
c_k = a_k + \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j.
\end{equation*}
```

""" tan!

@doc doc"""
    asin!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = asin(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = sqrt(1-c^2)` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = \frac{1}{ \sqrt{r_0} }
    \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big).
\end{equation*}
```

""" asin!

@doc doc"""
    acos!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = acos(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = sqrt(1-c^2)` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = - \frac{1}{ r_0 }
    \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big).
\end{equation*}
```

""" acos!

@doc doc"""
    atan!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = atan(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = 1+a^2` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = \frac{1}{r_0}\big(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j\big).
\end{equation*}
```

""" atan!

@doc doc"""
    sinhcosh!(s, c, a, k)

Update the `k-th` expansion coefficients `s[k+1]` and `c[k+1]`
of `s = sinh(a)` and `c = cosh(a)` simultaneously, for `s`, `c` and `a`
either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{eqnarray*}
s_k &=& \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j, \\
c_k &=& \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} s_j.
\end{eqnarray*}
```

""" sinhcosh!

@doc doc"""
    tanh!(c, a, p, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = tanh(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `p = a^2` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = a_k - \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j.
\end{equation*}
```

""" tanh!

@doc doc"""
    asinh!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = asinh(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = sqrt(1-c^2)` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = \frac{1}{ \sqrt{r_0} }
    \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big).
\end{equation*}
```

""" asinh!

@doc doc"""
    acosh!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = acosh(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = sqrt(c^2-1)` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = \frac{1}{ r_0 }
    \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big).
\end{equation*}
```

""" acosh!

@doc doc"""
    atanh!(c, a, r, k)

Update the `k-th` expansion coefficients `c[k+1]` of `c = atanh(a)`,
for `c` and `a` either `Taylor1` or `TaylorN`; `r = 1-a^2` and
is passed as an argument for efficiency.

```math
\begin{equation*}
c_k = \frac{1}{r_0}\big(a_k + \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j\big).
\end{equation*}
```

""" atanh!


text:
# make sure to include this file
# AFTER BarPhysics.jl has been included
module BarPhysicsInversion
using HiQGA.transD_GP.AbstractOperator
# using AbstractOperator (above) says that transD_GP can
# use this kind of physics, through a misfit function

import HiQGA.transD_GP.AbstractOperator.get_misfit
# get_misfit (above) extends all the other misfit functions
# in transD_GP. Because we are extending a method
# we must use import instead of using.

import HiQGA.transD_GP.Model, HiQGA.transD_GP.Options
# the above line now recognizes GP models and McMC options defined in src/TransD_GP_MCMC.jl

# The above three lines are the bare minimum we need to use transD_GP

# But wait! we need to specialize this inversion module
# for BarPhysics ... so let's do that. This is basically the
# geophysics "forward call"
using ..BarPhysics
import ..BarPhysics.returnphysics!

# Now let's actually define our BarInversion operator and for simplicity
# say it is a subtype of 1D operator (operators are defined in src/AbstractOperator.jl)
mutable struct BarInversion<:Operator1D
    d # some observations data
    F # some forward physics struct
end

# some constructor for BarInversion
function BarInversion()
    BarInversion([0, 0.], BarPhysics.Bar(20, rand(1:100))) # return value
end

# Now let's define (extend) a misfit function for
# a 1D model parameter vector m, the physics for which
# can be computed through methods in BarPhysics
# using an operator of type Bar, using McMC prior bounds
# and options stored in a struct of type Options
function get_misfit(m::Model, opt::Options, Bar::BarInversion)
    get_misfit(m.fstar, opt, Bar) # dispatches to helper function below
end
# above defined function and type signature MUST be defined

function get_misfit(x::AbstractArray, opt::Options, Bar::BarInversion)
    # m.fstar is the 1D array interpolated by the GP mean
    # at locations opt.xall
    returnphysics!(Bar.F, x) # think of this as F(x)
    chi2by2 = 0*Bar.d[1] - 0*Bar.F.somefield[1]
end
# this misfit function returns the same value for very call
# so we're sampling the prior
# this type signature is easier to debug than feeding in a GP model

end

ext: jl

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
# test_maps.jl

function suitable_point_to_map(m)
    T = domaintype(m)
    if T <: SVector
        randvec(eltype(T),length(T))
    else
        x = T(rand())
    end
end

suitable_point_to_map(m::DomainSets.EmbeddingMap{S,T}) where {S,T} = one(S)

function suitable_point_to_map(m::DomainSets.EmbeddingMap{SVector{N,S},T}) where {S,T,N}
    x = @SVector ones(N)
end

function suitable_point_to_map(m::DomainSets.ProductMap)
    x = ()
    for map in elements(m)
        x = (x..., suitable_point_to_map(map))
    end
    x
end

suitable_point_to_map(::CartToPolarMap{T}) where {T} = randvec(T,2)
suitable_point_to_map(::PolarToCartMap{T}) where {T} = randvec(T,2)
# Test a map m with dimensions n
function test_generic_map(T, m)
    # Try to map a random vector
    isreal(T(0)) && (@test isreal(m))
    x = suitable_point_to_map(m)
    y1 = applymap(m, x)
    y2 = m * x
    @test y1 == y2
    y3 = m(x)
    @test y1 == y3

    mi = inv(m)
    xi1 = applymap(mi, y1)
    @test xi1 ≈ x
    xi2 = mi * y1
    @test xi2 ≈ x
    xi3 = m\y1
    @test xi3 ≈ x
    xi4 = apply_inverse(m, y1)
    @test xi4 ≈ x

    # if is_linear(m)
    #     x = suitable_point_to_map(m, n)
    #     y1 = applymap(m, x)
    #     x0 = @SVector zeros(T,n)
    #     a,b = linearize(m, x0)
    #     y2 = a*x+b
    #     @test y1 ≈ y2
    # end
end

randvec(T,n) = SVector{n,T}(rand(n))
randvec(T,m,n) = SMatrix{m,n,T}(rand(m,n))

function test_maps(T)
    a = T(0)
    b = T(1)
    c = T(2)
    d = T(3)
    m = interval_map(a, b, c, d)
    @test m(a) ≈ c
    @test m(b) ≈ d

    test_generic_map(T, m)

    m2 = AffineMap(randvec(T, 2, 2), randvec(T, 2))
    test_generic_map(T, m2)

    m3 = LinearMap(randvec(T, 2, 2))
    test_generic_map(T, m3)
    r = randvec(T,2,1)
    @test jacobian(m3, r)*r ≈ m3*r

    # Test an affine map with a a scalar and b a vector
    m4 = AffineMap(T(1.2), randvec(T, 2))
    test_generic_map(T, m4)

    m5 = AffineMap(randvec(T, 3, 3), randvec(T, 3))
    test_generic_map(T, m5)

    m6 = m3∘m4
    @test typeof(m6) <: AffineMap
    test_generic_map(T, m6)

    # Test special maps
    test_composite_map(T)

    test_product_map(T)

    test_scaling_maps(T)

    test_identity_map(T)

    test_rotation_map(T)

    test_translation_map(T)

    test_cart_polar_map(T)

    # Diagonal maps

end

function test_scaling_maps(T)
    test_generic_map(T, scaling_map(T(2)))

    test_generic_map(T, scaling_map(T(2)))
    test_generic_map(T, scaling_map(T(2), T(3)))
    test_generic_map(T, scaling_map(T(2), T(3), T(4)))
    test_generic_map(T, scaling_map(T(2), T(3), T(4), T(5)))
end

function test_identity_map(T)
    test_generic_map(T, IdentityMap{T}())
    test_generic_map(T, IdentityMap{SVector{2,T}}())
    @test islinear(IdentityMap{T}())
end

function test_rotation_map(T)
    ϕ = T(pi)/4
    m = rotation_map(ϕ)
    x = [one(T), zero(T)]
    y = m*x
    @test y[1] ≈ sqrt(T(2))/2
    @test y[2] ≈ sqrt(T(2))/2

    ϕ = T(pi)/4
    m = rotation_map(ϕ, 0, 0)
    x = [zero(T), one(T), zero(T)]
    y = m*x
    @test y[1] ≈ 0
    @test y[2] ≈ sqrt(T(2))/2
    @test y[3] ≈ sqrt(T(2))/2

    # TODO: add more tests for a 3D rotation

    theta = T(rand())
    phi = T(rand())
    psi = T(rand())
    m2 = rotation_map(theta)
    test_generic_map(T, m2)
    m3 = rotation_map(phi, theta, psi)
    test_generic_map(T, m3)

    r = suitable_point_to_map(m2)
    @test norm(m2*r)≈norm(r)

    r = suitable_point_to_map(m3)
    @test norm(m3*r)≈norm(r)
    @test islinear(m3)
end

function test_translation_map(T)
    v = randvec(T,3)
    m = translation_map(v)
    test_generic_map(T, m)
    @test islinear(m)
end

function test_cart_polar_map(T)
    m1 = CartToPolarMap{T}()
    test_generic_map(T, m1)
    @test !islinear(m1)

    m2 = PolarToCartMap{T}()
    test_generic_map(T, m2)
    @test !islinear(m2)
end

function test_embedding_map(T1, T2)
    x = nonzero_element(T1)
    m = embedding_map(T1, T2)

    x = suitable_point_to_map(m)
    y1 = applymap(m, x)
    y2 = m * x
    @test y1 == y2
    y3 = m(x)
    @test y1 == y3

    @test applymap(m, x) == convert_space(spacetype(T2), x)
    m2 = restriction_map(T2, T1)
    y = nonzero_element(T2)
    @test applymap(m2, y) == restrict_space(spacetype(T1), y)
end

function test_isomorphism_map(T1, T2)
    test_embedding_map(T1, T2)
    test_embedding_map(T2, T1)
    @test inv(isomorphism_map(T1,T2)) == isomorphism_map(T2, T1)
    x = nonzero_element(T1)
    m = isomorphism_map(T1, T2)
    y = applymap(m, x)
    @test apply_inverse(m, y) == x
    m2 = isomorphism_map(T2, T1)
    y2 = nonzero_element(T2)
    x2 = applymap(m2, y2)
    @test apply_inverse(m2, x2) == y2
end

function test_composite_map(T)
    a = T(0)
    b = T(1)
    c = T(2)
    d = T(3)
    ma = IdentityMap{T}()
    mb = interval_map(a, b, c, d)

    r = suitable_point_to_map(ma)
    m1 = ma∘mb
    test_generic_map(T, m1)
    @test m1(r) ≈ ma(mb(r))
    m2 = m1∘mb
    test_generic_map(T, m2)
    @test m2(r) ≈ m1(mb(r))
    m3 = mb∘m2
    test_generic_map(T, m3)
    @test m3(r) ≈ mb(m2(r))
    m = m2∘m3
    test_generic_map(T, m)
    @test m(r) ≈ m2(m3(r))

    @test !(typeof(CompositeMap(ma)) <: CompositeMap)
end

function test_product_map(T)
    a = T(0)
    b = T(1)
    c = T(2)
    d = T(3)
    ma = IdentityMap{T}()
    mb = interval_map(a, b, c, d)

    r1 = suitable_point_to_map(ma)
    r2 = suitable_point_to_map(ma)
    r3 = suitable_point_to_map(ma)
    r4 = suitable_point_to_map(ma)
    r5 = suitable_point_to_map(ma)

    m1 = tensorproduct(ma,mb)
    test_generic_map(T, m1)
    @test compare_tuple(m1((r1,r2)), (ma(r1),mb(r2)))
    m2 = tensorproduct(m1,mb)
    test_generic_map(T, m2)
    @test compare_tuple(m2((r1,r2,r3)), (ma(r1),mb(r2),mb(r3)) )
    m3 = tensorproduct(mb,m2)
    test_generic_map(T, m3)
    @test compare_tuple(m3((r1,r2,r3,r4)),(mb(r1),ma(r2),mb(r3),mb(r4)))
    m = tensorproduct(m1,m2)
    test_generic_map(T, m)
    @test compare_tuple(m((r1,r2,r3,r4,r5)),(m1((r1,r2))...,m2((r3,r4,r5))...))
end

Base.isapprox(a::NTuple{L,SVector{N,T}}, b::NTuple{L,SVector{N,T}}) where {L,N,T} = compare_tuple(a,b)
Base.isapprox(a::NTuple{L,T}, b::NTuple{L,T}) where {L,T} = compare_tuple(a,b)

# Compare two iterable sequences for element-wise equality
compare_tuple(a, b) = reduce(&, map(isapprox, a, b))

@testset "maps" begin
    test_maps(Float64)
    test_maps(BigFloat)

    @testset "embedding map" begin
        T1 = Float64
        T2 = Complex{Float64}
        T3 = SVector{1,Float64}
        T4 = SVector{2,Float64}
        T5 = SVector{2,Complex{Float64}}

        test_embedding_map(T1, T2)
        test_embedding_map(T3, T4)
        test_isomorphism_map(T2, T4)
    end
end


text:
println("# test.jl ...")
using Random
using Test
include("keys0.jl")

egall() = begin eglib(); egsome() end

function eglib()
  println("#-- lib")
  @assert same(1)==1 "not same"
  @assert int(3.2)==3  "bad int convert"
  @assert per([1,2,3,4], .25) == 1 "bad per"
  @assert thing(".32") == .32 "bad thing"
  Random.seed!(1)
  @assert few([10,20,30,40],4) == [10,20,20,10]
  @assert o([10,20,Dict(1=>2,3=>4)]) ==
            "[10, 20, {3=4, 1=2}]" "bad string"
  lst=[x for x in  csv("data/auto93.csv")]
  @assert typeof(lst[end][1]) == Float64 "badtype"
  @assert 399 == length(lst) "rows missed"
end

function egsym()
  println("#-- sym")
  s=inc!(Sym(), ["a","a","a","a","b","b","c"])
  @assert 1.378 <= var(s) <= 1.379
end

function egsome()
  println("#-- some")
  s=Some()
  Random.seed!(1)
  s=inc!(Some(),[Random.rand() for _ in 1:1000])
  @assert it.some.max == length(s._all)
  @assert 0.47 <= mid(s) <= 0.53
  @assert .3 <= var(s) <=.31
end

function egsomes()
  println("#-- somes")
  s=Some()
  Random.seed!(1)
  r(z) =  round(z,digits=3)
  for i in 1:10
    inc!(s,[Random.rand()^2 for _ in 1:100])
    @assert .3  <=  var(s)  <= .34
    @assert .21 <=  mid(s)  <= .26
  end
end

egsAll= Dict("$x"[3:end] => getfield(Main,x)
             for x in names(Main)
             if length("$x")>2 && "$x"[1:2]=="eg")

for (n,item) in enumerate(ARGS)
 if item=="-t" && haskey(egsAll,ARGS[n+1])
   egsAll[ARGS[n+1]]()
   exit()
 end
 if  item=="-T"
   for (_,f) in egs f() end
 end
end

max_stars_repo_name: timm/peek

max_issues_repo_name: timm/peek

max_forks_repo_name: timm/peek


text:
"""
    check_feasibility(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point) where T<:AbstractNodeLabel

Evaluate the constraints defined by `elemid` of type `elem` at `bus` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if there is not violated constraints.

"""

###Genbounds constraints
function check_feasibility(element::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:GOCGenerator
    if elem_formulation == :NewVar
      lb = element.power_min
      ub = element.power_max
      cstrname = get_cstrname(scenario, bus, elemid, "Genbounds")
      varname = variable_name("Sgen", bus, elemid, scenario)
      value = point[Variable(varname,Complex)]
      #println(lb, " <= ", value, " <= ", ub)
      slack_left = value - lb
      slack_right = ub - value
      min_slack = min(real(slack_left),real(slack_right),imag(slack_left),imag(slack_right))
      if min_slack < -epsilon
          #infeasibility
          message = "Power generated not in diagram PQ."
          if real(slack_left) < - epsilon
              message = message*" Pgen < Pmin"
          elseif imag(slack_left) < - epsilon
              message = message*" Qgen < Qmin"
          elseif real(slack_right) < - epsilon
              message = message*" Pgen > Pmax"
          elseif imag(slack_left) < - epsilon
              message = message*" Qgen > Qmax"
          end
          return Dict(cstrname => message)
      end
    end
end

# ##test
# scenario = "BaseCase"
# bus = "BUS_1"
# ds_bus = OPFpbs[scenario].ds.bus[bus]
# elemid = "Generator_1"
# elem = ds_bus[elemid]
# elem_formulation = OPFpbs[scenario].mp.node_formulations[bus][elemid]
# point = global_point
# point_infeasible = copy(global_point)
# point_infeasible[Variable(variable_name("Sgen", bus, elemid, scenario), Complex)] += 1000
#
# check_feasibility(elem, bus, elemid, elem_formulation, scenario, point,1e-6)
# check_feasibility(elem, bus, elemid, elem_formulation, scenario, point_infeasible,1e-6)

###VOLTM constraints
function check_feasibility(element::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:GOCVolt
    cstrname = get_cstrname(scenario, bus, elemid, "VOLTM")
    Vmin = element.voltage_magnitude_min
    Vmax = element.voltage_magnitude_max
    varname = variable_name("VOLT", bus, "", scenario)
    value = point[Variable(varname,Complex)]
    module_value = value * conj(value)
    if abs(imag(module_value)) > 1e-15
        warn("Imaginary part of $(varname) module is not zero")
    end
    #println(Vmin^2, " <= ", module_value, " <= ", Vmax^2)
    slack_left = module_value - Vmin^2
    slack_right = Vmax^2 - module_value
    if real(slack_left) > - epsilon
        if real(slack_right) > - epsilon
            ##feasibility at epsilon OK
            return
        else
            return Dict(cstrname => "Voltage magnitude above Vmax")
        end
    else
        return Dict( cstrname => "Voltage magnitude under Vmin")
    end
end

##test
# scenario = "BaseCase"
# bus = "BUS_1"
# ds_bus = OPFpbs[scenario].ds.bus[bus]
# elemid = "VoltGOC"
# elem = ds_bus[elemid]
# elem_formulation = OPFpbs[scenario].mp.node_formulations[bus][elemid]
# point = global_point
#
# check_feasibility(elem, bus, elemid, elem_formulation, scenario, point, 1e-6)

# function check_feasibility(elem::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point) where T<:AbstractLinkLabel
#     return
# end
#
"""
    check_feasibility(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point, epsilon::Float64) where T<:AbstractLinkLabel

Evaluate the constraints defined by `elemid` of type `elem` at `link` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if there is not violated constraints.

"""

##Smax constraints TODO: define for all types of abstract link ???
function check_feasibility(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point, epsilon::Float64) where T<:GOCLineπ_notransformer
    origid, destid = link.orig, link.dest
    ctr_name_orig = get_cstrname(scenario, link, elemid, "Smaxorig")
    ctr_name_dest = get_cstrname(scenario, link, elemid, "Smaxdest")
    Sor = Sorig(element, link, elemid, elem_formulation, link_vars)
    Sde = Sdest(element, link, elemid, elem_formulation, link_vars)
    Smax = element.power_magnitude_max
    slack_Sorig = Smax^2 - evaluate(abs2(Sor), point)
    slack_Sdest = Smax^2 - evaluate(abs2(Sde), point)
    if abs(imag(slack_Sorig)) > 1e-15 || abs(imag(slack_Sdest)) > 1e-15
        warn("Imaginary part of Sorig or Sdest module is not zero")
    end
    ctr_slacks = Dict(ctr_name_orig => real(slack_Sorig), ctr_name_dest => real(slack_Sdest))
    min_violation_ctr = minimum(ctr_slacks)
    if min_violation_ctr[2] > - epsilon
        return
    else
        if real(slack_Sorig) < - epsilon
            if real(slack_Sdest) < - epsilon
                return Dict(ctr_name_orig => "Power magnitude at origin above Smax", ctr_name_dest => "Power magnitude at destination")
            else
                return Dict(ctr_name_orig => "Power magnitude at origin is above Smax")
            end
        else
            return Dict(ctr_name_dest => "Power magnitude at destination Smax")
        end
    end
end

function check_feasibility(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, point::Point, epsilon::Float64) where T<:GOCLineπ_withtransformer
    origid, destid = link.orig, link.dest
    ctr_name_orig = get_cstrname(scenario, link, elemid, "Smaxorig")
    ctr_name_dest = get_cstrname(scenario, link, elemid, "Smaxdest")
    Sor = Sorig(element, link, elemid, elem_formulation, link_vars)
    Sde = Sdest(element, link, elemid, elem_formulation, link_vars)
    Smax = element.power_magnitude_max
    slack_Sorig = Smax^2 - evaluate(abs2(Sor), point)
    slack_Sdest = Smax^2 - evaluate(abs2(Sde), point)
    if abs(imag(slack_Sorig)) > 1e-15 || abs(imag(slack_Sdest)) > 1e-15
        warn("Imaginary part of Sorig or Sdest module is not zero")
    end
    ctr_slacks = Dict(ctr_name_orig => real(slack_Sorig), ctr_name_dest => real(slack_Sdest))
    min_violation_ctr = minimum(ctr_slacks)
    if min_violation_ctr[2] > - epsilon
        return
    else
        if real(slack_Sorig) < - epsilon
            if real(slack_Sdest) < - epsilon
                return Dict(ctr_name_orig => "Power magnitude at origin above Smax", ctr_name_dest => "Power magnitude at destination")
            else
                return Dict(ctr_name_orig => "Power magnitude at origin is above Smax")
            end
        else
            return Dict(ctr_name_dest => "Power magnitude at destination Smax")
        end
    end
end

#test
# scenario = "BaseCase"
# link = collect(keys(OPFpbs[scenario].ds.link))[1]
# ds_link = OPFpbs[scenario].ds.link[link]
# elemid = linepi_notransformer_name("BL")
# elem = ds_link[elemid]
# link_vars = OPFpbs[scenario].mp.link_vars[link]
# elem_formulation = OPFpbs[scenario].mp.link_formulations[link][elemid]
# point = global_point
# check_feasibility(elem, link, elemid, elem_formulation, link_vars, scenario, point, 1e-6)

"""
    check_feasibility_cc_active_power(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:AbstractNodeLabel

Evaluate the coupling constraints about active power defined by `elemid` of type `elem` at `bus` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if the elem is not a generator or if there is not violated constraints.

"""
function check_feasibility_cc_active_power(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:AbstractNodeLabel
    return
end

function check_feasibility_cc_active_power(elem::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:GOCGenerator
    if !is_basecase_scenario(scenario)
        cc_name = cc_active_power(scenario, bus, elemid)
        part_factor = elem.participation_factor
        delta_name = delta_varname(scenario)
        sgen_bc_varname = variable_name("Sgen", bus, elemid, basecase_scenario_name())
        sgen_scen_varname = variable_name("Sgen", bus, elemid, scenario)
        delta = point[Variable(delta_name, Real)]
        Sgen_basecase = point[Variable(sgen_bc_varname,Complex)]
        Sgen_scenario = point[Variable(sgen_scen_varname,Complex)]
        diff = real(Sgen_scenario - Sgen_basecase - delta*part_factor)
        if abs(diff) < epsilon
            #feasibility
            return
        else
            return Dict(cc_name => "Active power generated in contingency not equal to the one in Basecase plus delta * participation factor")
        end
    else
        #warn("No active power coupling constraints for BaseCase scenario")
        return
    end
end

# scenario = "Scenario1"
# check_feasibility_cc_active_power(elem, bus, elemid, elem_formulation, scenario, point,1e-15)
"""
    check_feasibility_cc_reactive_power(element::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:AbstractNodeLabel

Evaluate the coupling constraints about reactive power defined by `elemid` of type `element` at `bus` for `scenario` (depending on `elem_formulation`) on `point` and return a dictionary of infeasible constraints Dict(ctrname => message with information)
Return nothing if the elem is not a generator or if there is not violated constraints.

"""
function check_feasibility_cc_reactive_power(element::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:AbstractNodeLabel
    return
end

function check_feasibility_cc_reactive_power(element::T, bus::String, elemid::String, elem_formulation::Symbol, scenario::String, point::Point,epsilon::Float64) where T<:GOCGenerator
    if !is_basecase_scenario(scenario)
        cc_name = cc_reactive_power(scenario, bus, elemid)
        sgen_scen_varname = variable_name("Sgen", bus, elemid, scenario)
        volt_bc_name = variable_name("VOLT", bus, "", basecase_scenario_name())
        volt_scen_name = variable_name("VOLT", bus, "", scenario)
        Sgen_scenario = point[Variable(sgen_scen_varname,Complex)]
        V_basecase = point[Variable(volt_bc_name,Complex)]
        V_scenario = point[Variable(volt_scen_name,Complex)]
        Qmin = imag(element.power_min)
        Qmax = imag(element.power_max)
        #println(abs(V_basecase)," ", abs(V_scenario))

        if  -epsilon <= abs(V_scenario) - abs(V_basecase) <= epsilon
            ##Vscen = Vbc => no constraint
            #println("Vscenario = Vbasecase")
            return
        elseif abs(V_scenario) - abs(V_basecase) < - epsilon
            #Vscen < Vbc => Qscen = Qmax
            println(bus, elemid, " Vscenario < Vbasecase")
            slack = Qmax - imag(Sgen_scenario)
            println(Qmax, imag(Sgen_scenario))
            println(slack)
            if abs(slack) > epsilon
                return Dict(cc_name => "Vscenario < Vbasecase, reactive power generated should be equal to Qmax")
            end

        elseif abs(V_scenario) - abs(V_basecase) > epsilon
            #Vscn > Vbc => Qscen = Qmin
            println(bus, elemid, " Vscenario > Vbasecase")
            slack = imag(Sgen_scenario) - Qmin
            println(Qmin, imag(Sgen_scenario))
            println(slack)
            if abs(slack) > epsilon
                return Dict(cc_name => "Vscenario > Vbasecase, reactive power generated should be equal to Qmin")
            end
        end

    else
        #warn("No reactive power coupling constraints for BaseCase scenario")
        return
    end
end

# scenario = "Scenario1"
# bus = "BUS_1"
# ds_bus = OPFpbs[scenario].ds.bus[bus]
# elemid = "Generator_1"
# elem = ds_bus[elemid]
# elem_formulation = OPFpbs[scenario].mp.node_formulations[bus][elemid]
# point = global_point
# check_feasibility_cc_reactive_power(elem, bus, elemid, elem_formulation, scenario, point,1e-7)

## TODO::
# function check_feasibility_power_balance(bus::String, ds::DataSource, gs::GridStructure, mp::MathematicalProgramming, scenario::String, point::Point, epsilon::Float64)
#
#     bus_elems_formulations = mp.node_formulations[bus]
#     bus_elems_var = mp.node_vars[bus]
#     cstrnames, Snode, lb, ub = get_Snodal(busid, ds.bus[bus], bus_elems_formulations, bus_elems_var)
#     S_balance = Snode + get_Slinks_in(busid, ds, gs, mp) + get_Slinks_out(busid, ds, gs, mp)
#
#
# end
"""
    test_feasibility_GOC(filenames, epsilon::Float64)

Check feasibility at `epsilon` for instance defined by `filenames`
Return a dictonary of violated constraints by scenario : dict( scenario => dict(ctrname => message))

"""

function test_feasibility_GOC(instance_path, epsilon::Float64)
    OPFpbs = load_OPFproblems(GOCInput, instance_path)
    introduce_Sgenvariables!(OPFpbs)
    pb_global = build_globalpb(OPFpbs)
    global_point = solution_point(instance_path)
    dict_scenario_infeasible_ctr = Dict{String, Dict{String,String}}()
    for (scenario, OPFpb) in OPFpbs
        dict_scenario_infeasible_ctr[scenario] = check_feasibility(scenario, OPFpb, global_point,epsilon)
    end
    return dict_scenario_infeasible_ctr
end

"""
    check_feasibility(scenario::String, OPFpb::Scenario, point::Point, epsilon::Float64)

Checks feasibility at `epsilon` for each constraint generated for in the scenario `scenario` constructed from `OPFpb` evaluated at `point`.
Returns a dictionary of constraints violated Dict: ctrname => message of error
"""
function check_feasibility(scenario::String, OPFpb::Scenario, point::Point, epsilon::Float64)
    ds = OPFpb.ds
    mp = OPFpb.mp

    dict_infeasible_ctr = Dict{String, String}()
        for (bus, elems) in ds.bus
            bus_elems_formulations = mp.node_formulations[bus]
            for (elemid, elem) in elems
                elem_formulation = bus_elems_formulations[elemid]
                dict2 = check_feasibility(elem, bus, elemid, elem_formulation, scenario, point, epsilon)
                if typeof(dict2) != Void
                    merge!(dict_infeasible_ctr, dict2)
                end
                dict3 = check_feasibility_cc_active_power(elem, bus, elemid, elem_formulation, scenario, point, epsilon)
                if typeof(dict3) != Void
                    merge!(dict_infeasible_ctr, dict3)
                end
                dict4 = check_feasibility_cc_reactive_power(elem, bus, elemid, elem_formulation, scenario, point, epsilon)
                if typeof(dict4) != Void
                    merge!(dict_infeasible_ctr, dict4)
                end
            end
        end

        for (link, elems) in ds.link
            link_elems_formulations = mp.link_formulations[link]
            link_vars = mp.link_vars[link]
            for (elemid, elem) in elems
                elem_formulation = link_elems_formulations[elemid]
                dict2 = check_feasibility(elem, link, elemid, elem_formulation, link_vars, scenario, point, epsilon)
                if typeof(dict2) != Void
                    merge!(dict_infeasible_ctr, dict2)
                end
            end
        end
    return dict_infeasible_ctr
end

function is_basecase_scenario(scenario::String)
    if scenario == basecase_scenario_name()
        return true
    else
        return false
    end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Asciicast
using Test

@testset "Asciicast.jl" begin
    # Write your tests here.
end


text:
# This is a demo for solving shock structure problem with neural networks

using Kinetic
using KitML.DiffEqFlux
using OrdinaryDiffEq
using Plots

set = Setup(
    "gas",
    "shock",
    "1d1f1v",
    "kfvs",
    "shakhov",
    1,
    1,
    "vanleer",
    "fix",
    0.5,
    10.0,
)

ps = PSpace1D(-25, 25, 50, 0)
gas = Gas(1.0, 2.0, 2/3, 0.0, 3.0, 0.81, 1.0, 0.5)
vs = VSpace1D(-10, 10, 72)
wL, primL, fL, bcL, wR, primR, fR, bcR = ib_rh(gas.Ma, gas.γ, vs.u)
ib = IB1F(wL, primL, fL, bcL, wR, primR, fR, bcR)
ks = SolverSet(set, ps, vs, gas, ib, pwd())

ctr = Array{ControlVolume1D1F}(undef, ks.pSpace.nx)
face = Array{Interface1D1F}(undef, ks.pSpace.nx + 1)
for i in eachindex(ctr)
    if i <= ks.pSpace.nx ÷ 2
        ctr[i] = ControlVolume1D1F(
            ks.pSpace.x[i],
            ks.pSpace.dx[i],
            Float32.(ks.ib.wL),
            Float32.(ks.ib.primL),
            Float32.(ks.ib.fL),
        )
    else
        ctr[i] = ControlVolume1D1F(
            ks.pSpace.x[i],
            ks.pSpace.dx[i],
            Float32.(ks.ib.wR),
            Float32.(ks.ib.primR),
            Float32.(ks.ib.fR),
        )
    end
end
for i = 1:ks.pSpace.nx+1
    face[i] = Interface1D1F(ks.ib.wL, ks.ib.fL)
end

dt = timestep(ks, ctr, 0.0)
sumRes = zeros(3)
for iter = 1:123
    Kinetic.evolve!(ks, ctr, face, dt)
    Kinetic.update!(ks, ctr, face, dt, sumRes)
end

X = Array{Float32}(undef, ks.vSpace.nu, ks.pSpace.nx)
for i = 1:ks.pSpace.nx
    X .= ctr[i].f
end

function shakhov!(df, f, p, t)
    M, tau = p
    df .= (M .- f) ./ tau
end

M = Array{Float32}(undef, vs.nu, size(X, 2))
SM = Array{Float32}(undef, vs.nu, size(X, 2))
τ = Array{Float32}(undef, 1, size(X, 2))
for i in axes(X, 2)
    M[:, i] .= maxwellian(ks.vSpace.u, ctr[i].prim)
    q = heat_flux(ctr[i].f, ctr[i].prim, ks.vSpace.u, ks.vSpace.weights)
    S = shakhov(ks.vSpace.u, M[:,i], q, ctr[i].prim, ks.gas.Pr)
    SM[:,i] .= M[:,i] .+ S
    τ[1, i] = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
end
P = [SM, τ]

tspan = (0, dt)
prob = ODEProblem(shakhov!, X, tspan, P)
Y = solve(prob, Euler(), dt=dt) |> Array

model_univ = FastChain(FastDense(vs.nu, vs.nu * 2, tanh), FastDense(vs.nu * 2, vs.nu))

p_model = initial_params(model_univ)

function dfdt(f, p, t)
    df = (M .- f) ./ τ .+ model_univ(f .- M, p)
end

prob_ube = ODEProblem(dfdt, X, tspan, p_model)

function loss(p)
    sol_ube = solve(prob_ube, Euler(), u0=X, p=p, dt=dt)
    loss = sum(abs2, Array(sol_ube) .- Y)
    return loss
end

cb = function (p, l)
    display(l)
    return false
end

res = sci_train(loss, p_model, ADAM(), cb=Flux.throttle(cb, 1), maxiters=200)

sol = solve(prob_ube, Euler(), u0=X, p=res.u, dt=dt)

contour(ks.ps.x, ks.vs.u, sol.u[end])


text:
using AbstractTrees

struct TrajectoryLog
    x
    v
    chart

    g
    ginv
    Γ

    JftWJf
    JftWginvℱ
    A
    B
    vmΓmvm

    TrajectoryLog() = new([],[],[],[],[],[],[],[],[],[],[])
end

function clear!(traj_log::TrajectoryLog)
    empty!(traj_log.x)
    empty!(traj_log.v)
    empty!(traj_log.chart)
    empty!(traj_log.g)
    empty!(traj_log.ginv)
    empty!(traj_log.Γ)
    empty!(traj_log.JftWJf)
    empty!(traj_log.JftWginvℱ)
    empty!(traj_log.A)
    empty!(traj_log.B)
    empty!(traj_log.vmΓmvm)
end

mutable struct TreeNode{T}
    data::T # Manifold for root, Task for leaves, TaskMap otherwise
    children::Vector{TreeNode}
    chart::Chart
    coord_rep::CoordinateRep
    traj_log::TrajectoryLog
    parent::TreeNode

    function TreeNode(data::M; CM::Chart{I,M}=default_chart(M),
            robot_coord_rep=ChartRep()) where {M<:Manifold,I}
        new{M}(data, [], CM, robot_coord_rep, TrajectoryLog())
    end
    function TreeNode(data::T; CN=codomain_manifold(data),
            task_coord_rep=ChartRep()) where T<:Union{TaskMapX,TaskX}
        new{T}(data, [], CN, task_coord_rep, TrajectoryLog())
    end
end

AbstractTrees.children(node::TreeNode) = node.children
AbstractTrees.printnode(io::IO, node::TreeNode{T}) where T = print(io, T)
function AbstractTrees.printnode(io::IO, node::TreeNode{T}, inds) where T
    if !isempty(inds)
        print(io, '(')
        show(IOContext(io, :compact => true), inds[end])
        print(io, ") ")
    end
    AbstractTrees.printnode(io, T)
end

function add_child!(node::TreeNode{M}, taskx::Union{TaskMapX,TaskX};
        CN=default_chart(codomain_manifold(taskx)), task_coord_rep=ChartRep()) where M <: Manifold
    if domain_manifold(taskx) == M
        new_node = TreeNode(taskx; CN, task_coord_rep)
        push!(node.children, new_node)
        return new_node
    else
        throw(ArgumentError("Domain manifold mismatch"))
    end
    nothing
end

function add_child!(node::TreeNode{<:TaskMapX}, taskx::Union{TaskMapX,TaskX};
        CN=default_chart(codomain_manifold(taskx)), task_coord_rep=ChartRep()) where M <: Manifold
    if codomain_manifold(node.data) == domain_manifold(taskx)
        new_node = TreeNode(taskx; CN, task_coord_rep)
        push!(node.children, new_node)
        return new_node
    else
        throw(ArgumentError("Manifold mismatch"))
    end
    nothing
end

add_child!(node::TreeNode{<:TaskX}, args...) = throw(ArgumentError("All tasks must be leaf nodes"))

function check_tree(root::TreeNode)
    for node in Leaves(root)
        if typeof(node.data) <: Union{Manifold, TaskMapX}
            error("Check failed: All leaf nodes must be tasks")
        end
    end
    println("Check passed")
end

function set_charts_emb!(node::TreeNode{M}, xme, CM) where M <: Manifold
    !isglobal(node.chart) ? node.chart = choose_chart_emb(xme, node.chart) : nothing
    for child in children(node)
        set_charts_emb!(child, xme, node.chart)
    end
end
function set_charts_chart!(node::TreeNode{M}, xm, CM) where M <: Manifold
    !isglobal(CM) ? node.chart = choose_chart_chart(xm, CM) : nothing
    for child in children(node)
        set_charts_chart!(child, xm, node.chart)
    end
end
function set_charts_emb!(node::TreeNode, xme, CM)
    taskx = node.data
    !isglobal(node.chart) ? node.chart = choose_chart_emb(xme, taskx, CM, node.chart) : nothing
    for child in children(node)
        set_charts_emb!(child, task_map_emb(xme, taskx, CM, node.chart), node.chart)
    end
end
function set_charts_chart!(node::TreeNode, xm, CM)
    taskx = node.data
    !isglobal(node.chart) ? node.chart = choose_chart_chart(xm, taskx, CM, node.chart) : nothing
    for child in children(node)
        set_charts_chart!(child, task_map_chart(xm, taskx, CM, node.chart), node.chart)
    end
end

function clear_logs!(node::TreeNode)
    clear!(node.traj_log)
    for child in node.children
        clear_logs!(child)
    end
    nothing
end


text:
mutable struct Controller1 <: Controller
    g::Real
    m::Real
    M::Real
    l::Real

    c::Real
    k::Real
    Fs::Real

    F::Function
    control!::Function

    function Controller1(F,g,m,M,l,c,k,Fs)

        new(g,m,M,l,c,k,Fs,F,control!)
    end
end

function control!(mechanism,controller::Controller1,kstep)
    g = controller.g
    m = controller.m
    M = controller.M
    l = controller.l
    c = controller.c
    k = controller.k
    Fs = controller.Fs
    F = controller.F
    cart = mechanism.bodies[1]
    pend = mechanism.bodies[2]

    v = cart.state.vc[2]
    q = pend.state.qc
    θ = (rotation_angle(q) * rotation_axis(q))[1]
    ω = pend.state.ωc[1]
    V = 2/3*m*l^2*ω^2 + m*g*l*(1+cos(θ))

    Fcart = F(θ,ω,V)
    Fcart, Fpend = friction(Fcart, 0, v, ω, k, c, Fs)

    setForce!(mechanism, geteqconstraint(mechanism, 3), [Fcart])
    setForce!(mechanism, geteqconstraint(mechanism, 4), [Fpend])

    return
end

function run1!(str; vis::Bool=true)
    open("Files/E1.jl",create=true,write=true) do file
        write(file,"
        Δt = 0.01

        # Parameters from Table 1.1
        g = 9.81
        L1 = 0.28125
        l1 = L1/2
        l = l1 #0.14344
        M = 4.5
        m = 0.2585
        # I = 1/3*m*l^2
        c = 0.002
        k = 0.3
        Fs = 8

        # Joint axes
        ex = [1.0;0.0;0.0]
        ey = [0.0;1.0;0.0]

        p = [0.0;0.0;l1] # joint connection point

        # Desired orientation
        θ1 = 0

        # Links
        origin = Origin{Float64}()
        cart = Box(0.1, 0.5, 0.1, M)
        pend = Box(0.055, 0.055, L1, m)

        # Constraints
        joint1 = EqualityConstraint(Prismatic(origin, cart, ey))
        joint2 = EqualityConstraint(Revolute(cart, pend, ex; p2 = -p))

        links = [cart;pend]
        constraints = [joint1;joint2]

        mech = Mechanism(origin, links, constraints, Δt = Δt, g = -g)

        xinit = 0.0
        theta1init = pi
        F(theta,dtheta,Vis) = 0
        sgn(x) = x==0 ? 1 : sign(x)

        ### Begin Student Input

        ")
        write(file, str)
        write(file, "

        ### End Student Input

        setPosition!(origin,cart,Δx = [0;xinit;0])
        setPosition!(cart,pend,p2 = -p, Δq = UnitQuaternion(RotX(theta1init)))

        controller = Controller1(F,g,m,M,l,c,k,Fs)

        steps = Base.OneTo(Int(10/Δt))
        storage = Storage{Float64}(steps,2)

        try
            simulate!(mech,storage,controller,record = true)
        catch
            println(\"Error. Unstable behavior.\")
        end")

        if vis
            write(file, "
            visualize(mech,storage)")
        end
    end

    try
        include("Files/E1.jl")
    catch
        println("Error. Bad code.")
    end

    return
end


text:
using BenchmarkTools;
using Revise;

#########################################################################################
# HERE DEFINING THE NAMED TUPLE ON THE FLY FIXES THE ISSUE
#########################################################################################
# function L0(arg1; kwarg1, kwarg2, kwarg3, kwarg4)
#     arg1 + kwarg1 + kwarg2 + kwarg3 + kwarg4
# end

# allkwargs = (kwarg1=1, kwarg2=2, kwarg3=3, kwarg4=4)
# firstkwargs = (kwarg1=1, kwarg2=2)
# lastkwargs = (kwarg3=3, kwarg4=4)

# L1_fromglob(arg1; kwargs...) = L0(arg1; firstkwargs..., kwargs...)
# L1(arg1; kwargs...) = L0(arg1; (kwarg1=1, kwarg2=2)..., kwargs...)

# @btime L0(0; $allkwargs...) # 1.506 ns (0 allocations: 0 bytes)
# @btime L1_fromglob(0; $lastkwargs...) # 108.332 ns (2 allocations: 80 bytes)
# @btime L1(0; $lastkwargs...) # 1.205 ns (0 allocations: 0 bytes)

########################################################################################
# TESTING ALL THE DIFFERENT CASES...
########################################################################################
# function L0(arg1; kwarg1, kwarg2, kwarg3, kwarg4)
#     arg1 + kwarg1 + kwarg2 + kwarg3 + kwarg4
# end

# function freezeparams(func_to_freeze, params_to_freeze)
#     (arg1; kwargz...) -> func_to_freeze(arg1; kwargz..., params_to_freeze...)
# end

# function freezeparams_let(func_to_freeze::T, params_to_freeze::NamedTuple) where T<:Function
#     f = let func_to_freeze=func_to_freeze, params_to_freeze=params_to_freeze, arg1, kwargs, kwarg1, kwarg2, kwarg3, kwarg4
#         (arg1; kwargs...) -> func_to_freeze(arg1; kwargs..., params_to_freeze...)
#     end
#     f
# end

# allkwargs = (kwarg1=1, kwarg2=2, kwarg3=3, kwarg4=4)
# firstkwargs = (kwarg1=1, kwarg2=2)
# lastkwargs = (kwarg3=3, kwarg4=4)

# L1_FP_WITHARG_GLOBAL(arg1; kwargs...) = freezeparams(L0, firstkwargs)
# L1_FP_WITHARG_INTERP(arg1; kwargs...) = freezeparams(L0, (kwarg1=1, kwarg2=2))

# L1_FPL_WITHARG_GLOBAL(arg1; kwargs...) = freezeparams_let(L0, firstkwargs)
# L1_FPL_WITHARG_INTERP(arg1; kwargs...) = freezeparams_let(L0, (kwarg1=1, kwarg2=2))

# println("\nReference case with no freezing")
# @btime L0(0; $allkwargs...) # 1.205 ns (0 allocations: 0 bytes)

# println("\nL1_FP_WITHARG_GLOBAL")
# @btime L1_FP_WITHARG_GLOBAL(0; $lastkwargs...) # 271.585 ns (2 allocations: 64 bytes)

# println("\nL1_FP_WITHARG_INTERP")
# @btime L1_FP_WITHARG_INTERP(0; $lastkwargs...) # 1.506 ns (0 allocations: 0 bytes)

# println("\nL1_FPL_WITHARG_GLOBAL")
# @btime L1_FPL_WITHARG_GLOBAL(0; $lastkwargs...) # 269.493 ns (2 allocations: 64 bytes)

# println("\nL1_FPL_WITHARG_INTERP")
# @btime L1_FPL_WITHARG_INTERP(0; $lastkwargs...) # 1.506 ns (0 allocations: 0 bytes)

# println("\ntest done\n")

#########################################################################################
# Trying succesful approach above but within struct... Still doesn't work nicely.
#########################################################################################
# struct thang{T<:Function}
#     a::T
# end

# function L0(arg1; kwarg1, kwarg2, kwarg3, kwarg4)
#     arg1 + kwarg1 + kwarg2 + kwarg3 + kwarg4
# end

# function freeze_params(struct_to_freeze, params_to_freeze)
#     newF(arg1; kwargs...) = let params_to_freeze = params_to_freeze
#                                 struct_to_freeze.a(arg1; kwargs..., params_to_freeze...)
#                             end
#     thang(newF)
# end

# thang1 = thang(L0)

# allkwargs = (kwarg1=1, kwarg2=2, kwarg3=3, kwarg4=4)
# firstkwargs = (kwarg1=1, kwarg2=2)
# lastkwargs = (kwarg3=3, kwarg4=4)

# thang2 = freeze_params(thang1, firstkwargs)
# thang2_explicit = freeze_params(thang1, (kwarg1=1, kwarg2=2))

# @btime thang1.a(0; $allkwargs...) # 112.925 ns (1 allocation: 48 bytes)
# @btime thang2.a(0; $lastkwargs...) # 178.169 ns (7 allocations: 224 bytes)
# @btime thang2_explicit.a(0; $lastkwargs...) # 179.551 ns (7 allocations: 224 bytes)

#########################################################################################
# Try overload get property/field operators for struct.
#########################################################################################
# import Base.getproperty

# struct thingy{T<:Function}
#     free_params::Tuple
#     fixed_params::NamedTuple
#     a::T
# end

# function Base.getproperty(inst::thingy, sym::Symbol)
#     if sym==:a
#         f = let inst::thingy=inst, sym::Symbol=sym, getfield=getfield, fixed::NamedTuple=inst.fixed_params, t::Float64, kwargs::NamedTuple
#                 (t; kwargs...) -> getfield(inst, sym)(t; kwargs..., fixed...)
#             end
#         return f
#     elseif sym==:free_params || sym==:fixed_params
#         return getfield(inst, sym)
#     end
# end

# function freezeparams(orig::thingy, tofix::NamedTuple)
#     newfree =  Tuple([i for i in orig.free_params if !(i in keys(tofix))])
#     newfixed = merge(orig.fixed_params, tofix)
#     return thingy(newfree, newfixed, orig.a)
# end

# function funk1(t; a, b, c)
#     t+a+b+c
# end

# t1 = thingy((:a, :b, :c), NamedTuple{}(), funk1)
# t2 = freezeparams(t1, (a=1,))

# tarray = collect(0.0:0.01:1000);

# @btime t1.a.($tarray; a=1, b=2, c=3);
# @btime t2.a.($tarray; b=2, c=3);

# println("\nDone")

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using BinDeps
using Compat

@BinDeps.setup

version = "1.18"
url = "https://github.com/google/leveldb/archive/v$(version).tar.gz"

libleveldbjl = library_dependency("libleveldbjl")

provides(Sources, URI(url), libleveldbjl, unpacked_dir="leveldb-$(version)")

leveldbbuilddir = BinDeps.builddir(libleveldbjl)
leveldbsrcdir = joinpath(BinDeps.depsdir(libleveldbjl),"src", "leveldb-$version")
leveldblibdir = BinDeps.libdir(libleveldbjl)
leveldblibfile = joinpath(leveldblibdir,libleveldbjl.name*".$(Libdl.dlext)")

provides(BuildProcess,
    (@build_steps begin
        GetSources(libleveldbjl)
        CreateDirectory(leveldblibdir)
        @build_steps begin
            ChangeDirectory(leveldbsrcdir)
            FileRule(leveldblibfile, @build_steps begin
                `make`
                `cp libleveldb.$(Libdl.dlext) $(leveldblibfile)`
            end)
        end
    end), libleveldbjl, os = :Unix)

@BinDeps.install [ :libleveldbjl => :libleveldbjl]


text:
function solve_a(data::String)::Int
	total = 0
	len_data = length(data)

	for i = 1:len_data
		if i != len_data
			target_index = i+1
		else
			target_index = 1
		end

		if data[i] == data[target_index]
			total += parse(Int, data[i])
		 end
	end

	total
end

function solve_b(data::String)::Int
	total = 0
	len_data = length(data)
	half_len = len_data ÷ 2

	for i = 1:len_data
		target_index = i + half_len
		if target_index > len_data
			target_index %= len_data
		end

		if data[i] == data[target_index]
			total += parse(Int, data[i])
		end
	end

	total
end

data = readline("input.txt")
sol_a = solve_a(data)
sol_b = solve_b(data)

println("sol a: $(sol_a)")
println("sol b: $(sol_b)")


text:
# Download Pfam
# =============

"""
Download a gzipped stockholm full alignment for the `pfamcode`.
The extension of the downloaded file is `.stockholm.gz` by default. The `filename` can be changed, but the `.gz` at the end is mandatory.
"""
function downloadpfam(pfamcode::ASCIIString; filename::ASCIIString="$pfamcode.stockholm.gz")
  if ismatch(r"^PF\d{5}$"i, pfamcode) # length(pfamcode)== 7 && ( pfamcode[1:2] == "PF" || pfamcode[1:2] == "pf" )
    #namegz = string(filename, ".gz")
    download(string("http://pfam.xfam.org/family/PF", pfamcode[3:end], "/alignment/full/gzipped"), filename) # namegz)
    #run(`gzip -d $namegz`)
  else
    throw( ErrorException( string(pfamcode, " is not a correct Pfam code") ) )
  end
end


text:
include("sosdemo7.jl")
include("sosdemo8.jl")
include("BPT12e399.jl")


text:
# Definition of StrBootTest "class" for holding intermediate results, with associated utilities and get functions

struct StrClust{T<:Real}
	N::Int; multiplier::T; even::Bool
	order::Vector{Int64}
	info::Vector{UnitRange{Int64}}
end

struct StrFE{T<:Real}
	is::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}
  wt::Vector{T}
  sqrtwt::Vector{T}
end

mutable struct StrEstimator{T<:AbstractFloat}
  isDGP::Bool; liml::Bool; fuller::T; κ::T
  R₁perp::Matrix{T}; Rpar::Matrix{T}

  kZ::Int64
  y₁::Vector{T}; ü₁::Vector{T}; u⃛₁::Vector{T}; β̈::Vector{T}; γ̈::Vector{T}; β̈₀::Vector{T}; invXXXy₁par::Vector{T}
  Yendog::Matrix{Bool}
  invZperpZperp::Symmetric{T,Matrix{T}}; invZperpZperpZperpX::Matrix{T}; XZ::Matrix{T}; YPXY::Symmetric{T,Matrix{T}}; R₁invR₁R₁::Matrix{T}
	restricted::Bool; RperpX::Matrix{T}; RperpXperp::Matrix{T}; RRpar::Matrix{T}; RparX::Matrix{T}; RparY::Matrix{T}; RR₁invR₁R₁::Matrix{T}
	∂β̈∂r::Matrix{T}; YY::Symmetric{T,Matrix{T}}; AR::Matrix{T}; XAR::Matrix{T}; R₁invR₁R₁Y::Matrix{T}; invXXXZ::Matrix{T}; Ü₂::Matrix{T}; Rt₁::Vector{T}
	invXX::Symmetric{T,Matrix{T}}; Y₂::Matrix{T}; X₂::Matrix{T}; invH::Symmetric{T,Matrix{T}}
	y₁par::Vector{T}; Xy₁par::Vector{T}
	A::Symmetric{T,Matrix{T}}; Z::Matrix{T}; Zperp::Matrix{T}; X₁::Matrix{T}
	WXAR::Matrix{T}; CT_XAR::Vector{Matrix{T}}

  # IV/GMM only
  ZZ::Symmetric{T,Matrix{T}}; XY₂::Matrix{T}; XX::Symmetric{T,Matrix{T}}; H_2SLS::Symmetric{T,Matrix{T}}; V::Matrix{T}; ZY₂::Matrix{T}; ZR₁ZR₁::Symmetric{T,Matrix{T}}; X₂ZR₁::Matrix{T}; ZR₁Y₂::Matrix{T}; X₁ZR₁::Matrix{T}
  ZR₁Z::Matrix{T}; X₂y₁::Vector{T}; X₁y₁::Vector{T}; Zy₁::Vector{T}; ZXinvXXXZ::Matrix{T}; H_2SLSmZZ::Symmetric{T,Matrix{T}}
  ZXinvXXXy₁par::Vector{T}; t₁Y::Vector{T}
  Y₂y₁::Vector{T}; twoZR₁y₁::Vector{T}; y₁y₁::T; y₁pary₁par::T; Y₂Y₂::Matrix{T}
  X₂y₁par::Vector{T}; X₁y₁par::Vector{T}; Zy₁par::Vector{T}; Y₂y₁par::Vector{T}
  Rperp::Matrix{T}; ZR₁::Matrix{T}
  kX::Int64; kX₁::Int64; kZperp::Int64
	Π̂::Matrix{T}
	S✻⋂ZperpZpar::Array{T,3}; S✻⋂ZperpY₂::Array{T,3}; S⋂y₁X₁::Matrix{T}; S⋂y₁X₂::Matrix{T}; S⋂Xy₁::Matrix{T}; S✻⋂ZperpZperp::Array{T,3}; ZperpX::Matrix{T}; S✻⋂Zperpy₁::Matrix{T}; S✻⋂XZR₁::Array{T,3}; S✻⋂XY₂::Array{T,3}; S✻⋂XZperp::Array{T,3}; S✻⋂Xy₁::Matrix{T}; S✻⋂XX::Array{T,3}; S✻⋂XZpar::Array{T,3}
	S✻⋂X₁Y₂::Array{T,3}; S✻⋂X₂Y₂::Array{T,3}; S✻ZparY₂::Array{T,3}; S✻y₁y₁::Vector{T}; S✻Zpary₁::Matrix{T}; S✻ZR₁y₁::Matrix{T}; S✻ZR₁Y₂::Array{T,3}; S✻ZR₁ZR₁::Array{T,3}; S✻ZR₁Z::Array{T,3}
	ZperpZR₁::Matrix{T}; S✻⋂ZperpZR₁::Array{T,3}
	S✻Y₂Y₂::Array{T,3}; S✻ZparZpar::Array{T,3}; S✻Y₂y₁::Matrix{T}
	Xpar₁::Matrix{T}
	invZperpZperpZperpX₁::Matrix{T}; invZperpZperpZperpX₂::Matrix{T}; invZperpZperpZperpy₁::Vector{T}; invZperpZperpZperpY₂::Matrix{T}; S✻UY₂::Matrix{T}; invZperpZperpZperpZpar::Matrix{T}; invZperpZperpZperpZR₁::Matrix{T}

  StrEstimator{T}(isDGP, liml, fuller, κ) where T<:AbstractFloat = new(isDGP, liml, fuller, κ, Matrix{T}(undef,0,0))
end

mutable struct StrBootTest{T<:AbstractFloat}
  R::Matrix{T}; r::Vector{T}; R₁::Matrix{T}; r₁::Vector{T}
  y₁::Vector{T}; X₁::Matrix{T}; Y₂::Matrix{T}; X₂::Matrix{T}
  wt::Vector{T}; fweights::Bool
  liml::Bool; fuller::T; κ::T; arubin::Bool
  B::Int64; auxtwtype::Symbol; rng::AbstractRNG; maxmatsize::Float16
  ptype::Symbol; null::Bool; bootstrapt::Bool
	ID::Matrix{Int64}; NBootClustVar::Int8; NErrClustVar::Int8; issorted::Bool; small::Bool; clusteradj::Bool; clustermin::Bool
  FEID::Vector{Int64}; FEdfadj::Int64
  level::T; rtol::T
  madjtype::Symbol; NH₀::Int16
  ml::Bool; β̈::Vector{T}; A::Symmetric{T,Matrix{T}}; sc::Matrix{T}
  willplot::Bool; gridmin::Vector{T}; gridmax::Vector{T}; gridpoints::Vector{Float32}

  q::Int16; twotailed::Bool; scorebs::Bool; robust::Bool

  WRE::Bool; initialized::Bool; NFE::Int64; FEboot::Bool; granular::Bool; NErrClustCombs::Int16; subcluster::Int8; BFeas::Int64; interpolating::Bool
  v_sd::T
  confpeak::Vector{T}
  ID✻::Vector{Int64}; ID✻_✻⋂::Vector{Int64}
  anchor::Vector{T}; poles::Vector{T}; numer::Matrix{T}
  ci::Matrix{T}
  peak::NamedTuple{(:X, :p), Tuple{Vector{T}, T}}

	Nobs::Int64; NClustVar::Int8; kX₁::Int64; kX₂::Int64; kY₂::Int64; WREnonARubin::Bool; boottest!::Function
	coldotplus!::Function; colquadformminus!::Function; matmulplus!::Function; panelsum!::Function

  sqrt::Bool; _Nobs::T; kZ::Int64; sumwt::T; haswt::Bool; sqrtwt::Vector{T}; REst::Bool; multiplier::T; smallsample::T
		dof::Int64; dof_r::T; p::T; BootClust::Int8
		purerobust::Bool; N✻::Int64; N⋂::Int64; N✻⋂::Int64; Nw::Int64; enumerate::Bool; interpolable::Bool; interpolate_u::Bool; kX::Int64
  _FEID::Vector{Int64}; AR::Matrix{T}; v::Matrix{T}; u✻::Matrix{T}
  info✻::Vector{UnitRange{Int64}}; info✻_✻⋂::Vector{UnitRange{Int64}}; infoBootAll::Vector{UnitRange{Int64}}; info⋂_✻⋂::Vector{UnitRange{Int64}}
  JN⋂N✻::Matrix{T}; statDenom::Matrix{T}; uXAR::Matrix{T}; SuwtXA::Matrix{T}; numer₀::Matrix{T}; β̈dev::Matrix{T}
	YY✻_b::Matrix{T}; YPXY✻_b::Matrix{T}; numerw::Matrix{T}; numer_b::Vector{T}; dist::Matrix{T}

	distCDR::Matrix{T}; plotX::Vector{Vector{T}}; plotY::Vector{T}; ClustShare::Vector{T}; WeightGrp::Vector{UnitRange{Int64}}
  numersum::Vector{T}; u✻₀::Matrix{T}; invFEwt::Vector{T}
	β̈s::Matrix{T}; As::Matrix{T}
	info✻⋂::Vector{UnitRange{Int64}}; info⋂::Vector{UnitRange{Int64}}; ID✻⋂::Matrix{T}
	ü::Vector{T}
	DGP::StrEstimator{T}; Repl::StrEstimator{T}; M::StrEstimator{T}
	clust::Vector{StrClust{T}}
	denom::Matrix{Matrix{T}}; Kcd::Matrix{Matrix{T}}; Jcd::Matrix{Matrix{T}}; denom₀::Matrix{Matrix{T}}; Jcd₀::Matrix{Matrix{T}}; S✻UU::Matrix{Vector{T}}; CTUX::Matrix{Matrix{T}}
	∂u∂r::Vector{Matrix{T}}; ∂numer∂r::Vector{Matrix{T}}; S✻XU::Vector{Matrix{T}}; invXXS✻XU::Vector{Matrix{T}}; invZperpZperpS✻ZperpU::Vector{Matrix{T}}; S✻YU::Matrix{Vector{T}}; S✻UMZperp::Vector{Matrix{T}}; S✻UPX::Vector{Matrix{T}}; S✻ZperpU::Vector{Matrix{T}}; CTFEU::Vector{Matrix{T}}
  ∂denom∂r::Array{Matrix{T},3}; ∂Jcd∂r::Array{Matrix{T},3}; ∂²denom∂r²::Array{Matrix{T},4}
	FEs::Vector{StrFE{T}}
  T1L::Vector{Matrix{T}}; T1R::Vector{Matrix{T}}; T2::Matrix{T}; J⋂s::Vector{Array{T,3}}; Q::Array{T,3}; β̈v::Vector{Matrix{T}}
	crosstab⋂✻ind::Vector{Int64}; crosstab✻ind::Vector{Int64}
  seed::UInt64

	S✻XY₂::Array{T,3}; S✻XX::Array{T,3}; S✻XDGPZ::Array{T,3}; S✻Xy₁::Matrix{T}; S✻XZR₁::Array{T,3}
	invXXS✻XY₂::Array{T,3}; invXXS✻XX::Array{T,3}; invXXS✻XDGPZ::Array{T,3}; invXXS✻Xy₁::Matrix{T}; invXXS✻XDGPZR₁::Array{T,3}
	S✻⋂XY₂::Array{T,3}; S✻⋂XX::Array{T,3}; S✻⋂XDGPZ::Array{T,3}; S✻⋂Xy₁::Matrix{T}; S✻⋂X_DGPZR₁::Array{T,3}
	invXXS✻⋂XY₂::Array{T,3}; invXXS✻⋂XX::Array{T,3}; invXXS✻⋂XDGPZ::Array{T,3}; invXXS✻⋂Xy₁::Matrix{T}; invXXS✻⋂XDGPZR₁::Array{T,3}
	S✻ZperpY₂::Array{T,3}; S✻ZperpX::Array{T,3}; S✻ZperpDGPZ::Array{T,3}; S✻Zperpy₁::Matrix{T}; S✻ZperpDGPZR₁::Array{T,3}
	invZperpZperpS✻ZperpY₂::Array{T,3}; invZperpZperpS✻ZperpX::Array{T,3}; invZperpZperpS✻ZperpDGPZ::Array{T,3}; invZperpZperpS✻Zperpy₁::Matrix{T}; invZperpZperpS✻ZperpDGPZR₁::Array{T,3}
	_ID✻⋂::Vector{Int}
	S✻Y₂y₁::Matrix{T}; S✻DGPZy₁::Matrix{T}; S✻y₁y₁::Vector{T}; S✻DGPZR₁y₁::Matrix{T}; r₁S✻ReplZR₁Y₂::Array{T,3}; r₁S✻ReplZR₁X::Array{T,3}; r₁S✻ReplZR₁DGPZ::Array{T,3}; r₁S✻ReplZR₁y₁::Matrix{T}; r₁S✻ReplZR₁DGPZR₁::Array{T,3}; S✻ReplZY₂::Array{T,3}; S✻ReplZX::Array{T,3}; S✻ReplZDGPZ::Array{T,3}; S✻ReplZy₁::Matrix{T}; S✻ReplZDGPZR₁::Array{T,3}
	S✻U_S✻UMZperpX::Matrix{Array{T,3}}; negS✻UMZperpX::Vector{Array{T,3}}; S⋂XZperpinvZperpZperp::Array{T,3}; CT✻FEX::Array{T,3}; CT⋂FEX::Array{T,3}; CT✻FEY₂::Array{T,3}; CT✻FEZ::Array{T,3}; CT✻FEy₁::Matrix{T}; CT✻FEZR₁::Array{T,3}
	S✻Y₂Y₂::Array{T,3}; S✻ZparY₂Z_DGPZ::Array{T,3}; S✻DGPZY₂::Array{T,3}; S✻DGPZR₁Y₂::Array{T,3}; S✻DGPZDGPZ::Array{T,3};  S✻DGPZR₁DGPZR₁::Array{T,3}; S✻DGPZR₁DGPZ::Array{T,3}; S✻DGPZR₁X::Array{T,3}
	XinvXX::Matrix{T}; PXZ::Matrix{T}; FillingT₀::Matrix{Matrix{T}}
	S⋂ReplZX::Array{T,3}; S⋂Xy₁::Matrix{T}
	S✻⋂XU₂::Array{T,3}; S✻⋂XU₂RparY::Array{T,3}; S✻XU₂::Array{T,3}; S✻XU₂RparY::Array{T,3}; S✻ZperpU₂::Array{T,3}; S✻ZperpU₂RparY::Array{T,3}; invZperpZperpS✻ZperpU₂::Array{T,3}; invZperpZperpS✻ZperpU₂RparY::Array{T,3}; invXXS✻XU₂::Array{T,3}; invXXS✻XU₂RparY::Array{T,3}

	StrBootTest{T}(R, r, R₁, r₁, y₁, X₁, Y₂, X₂, wt, fweights, liml,
	               fuller, κ, arubin, B, auxtwtype, rng, maxmatsize, ptype, null, scorebs, bootstrapt, ID, NBootClustVar, NErrClustVar, issorted, robust, small, clusteradj, clustermin,
								 NFE, FEID, FEdfadj, level, rtol, madjtype, NH₀, ml,
								 β̈, A, sc, willplot, gridmin, gridmax, gridpoints) where T<:Real =
		begin
			kX₂ = ncols(X₂)
			scorebs = scorebs || iszero(B) || ml
			WREnonARubin = !(iszero(kX₂) || scorebs) && !arubin

			new(R, r, R₁, r₁, y₁, X₁, Y₂, X₂, wt, fweights, liml || !iszero(fuller),
					fuller, κ, arubin, B, auxtwtype, rng, maxmatsize, ptype, null, bootstrapt, ID, NBootClustVar, NErrClustVar, issorted, small, clusteradj, clustermin,
					FEID, FEdfadj, level, rtol, madjtype, NH₀, ml,
					β̈, A, sc, willplot, gridmin, gridmax, gridpoints,
				nrows(R), ptype == :symmetric || ptype == :equaltail, scorebs, robust || NErrClustVar>0,
				false, false, NFE, false, false, 0, 0, 0, false,
				one(T),
				[zero(T)],
				Vector{Int64}(undef,0), Vector{Int64}(undef,0),
				Vector{T}(undef,0), Vector{T}(undef,0), Matrix{T}(undef,0,0),
				Matrix{T}(undef,0,0),
				(X = Vector{T}(undef,0), p = T(NaN)),
				nrows(X₁), ncols(ID), ncols(X₁), kX₂, ncols(Y₂), WREnonARubin, WREnonARubin ? boottestWRE! : boottestOLSARubin!,
				coldotplus_nonturbo!,
				colquadformminus_nonturbo!,
				matmulplus_nonturbo!,
				panelsum_nonturbo!)
		end
end

function getdist(o::StrBootTest, diststat::Symbol=:none)
  if diststat == :numer
	  _numer = isone(o.v_sd) ? o.numer : o.numer / o.v_sd
	  o.distCDR = (@view _numer[:,2:end]) .+ o.r
	  # sort!(o.distCDR, dims=1)
  elseif nrows(o.distCDR)==0  # return test stats
    if length(o.dist) > 1
	    o.distCDR = (@view o.dist[1,2:end])' * o.multiplier
	    # sort!(o.distCDR, dims=1)
	  else
	    o.distCDR = zeros(0,1)
	  end
  end
  o.distCDR
end

function sumgreater(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x > i && (dest += 1)
  end
  dest
end
function sumless(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x < i && (dest += 1)
  end
  dest
end
function sumlessabs(x, v)
  dest = zero(Int64)
  @inbounds @simd for i in v
	  x < abs(i) && (dest += 1)
  end
  dest
end

# store p value in o.p. Return optionally-multiple-hypothesis-adjusted p value. Robust to missing bootstrapped values interpreted as +infinity.
function getp(o::StrBootTest{T}) where T
  o.boottest!(o)
  tmp = o.dist[1]
  isnan(tmp) && return tmp
  if o.B>0
  	if o.sqrt && o.ptype ≠ :upper
  	  if o.ptype == :symmetric
  	    n = sumlessabs(abs(tmp), o.dist)   # symmetric p value; do so as not to count missing entries in *dist
  	  elseif o.ptype == :equaltail
  	    n = 2min(sumgreater(tmp, o.dist) , sumless(tmp, o.dist))
  	  else
  		  n = sumgreater(tmp,  o.dist)  # lower-tailed p value
      end
  	else
  	  n = sumless(tmp, o.dist)  # upper-tailed p value or p value based on squared stats
    end
    o.p = n / o.BFeas |> T
  else
		tmp *= o.multiplier
		o.p = ccdf(o.small ? FDist{T}(T(o.dof), o.dof_r) : Chisq{T}(T(o.dof)), o.sqrt ? tmp^2 : tmp)
		if o.sqrt && !o.twotailed
			o.p /= 2
			(o.ptype == :upper) == (tmp<0) && (o.p = 1 - o.p)
		end
  end

	if o.madjtype == :bonferroni min(one(T), o.NH₀ * o.p)
  elseif o.madjtype == :sidak  one(T) - (one(T) - o.p) ^ o.NH₀
  else o.p
  end
end

getb(o::StrBootTest) = isone(o.v_sd) ? o.numer[:,1] : o.numer[:,1] / o.v_sd  # numerator for full-sample test stat
getV(o::StrBootTest) = o.statDenom / ((isone(o.v_sd) ? o.smallsample : o.v_sd^2 * o.smallsample) * (o.sqrt ? o.multiplier^2 : o.multiplier) * o.dof)  # denominator for full-sample test stat
getv(o::StrBootTest) = @views isone(o.v_sd) ? o.v[:,2:end] : o.v[:,2:end] / o.v_sd  # wild weights
@inline getstat(o::StrBootTest) = o.multiplier * o.dist[1]


text:
module chapters_faq_gotchas

using Dash, DashHtmlComponents, DashCoreComponents

app =  dash()

app.layout = html_div() do
    html_h1("Interactive Visualizations"),
    html_blockquote(dcc_markdown("This is the 6th and final chapter of the [Dash Tutorial](/).
    The [previous chapter](/sharing-data-between-callbacks) covered basic callback usage and the [next chapter](/state)
    describes how to share data between callbacks. The [rest of the Dash documentation](/) covers other
    topics like multi-page apps and component libraries.
    Just getting started? Make sure to [install the necessary
    dependencies](/installation)")),
    dcc_markdown("
    ### Frequently Asked Questions

    Q: How can I customize the appearance of my Dash app?
    A: Dash apps are rendered in the browser as modern standards-compliiant web apps. This means
    you can use CSS to style your Dash app as you would standard HTML.

    All `DashHtmlComponents` support inline CSS styling through a `style` attribute. An external
    CSS stylesheet can also be used to style `DashHtmlComponents` and `DashCoreComponents` by
    targeting the IDs or classnames of your components. Both `DashHtmlComponents` and `DashCoreComponents`
    accept the attribute `className`, which corresponds to the HTML element attribute `class`.

    To use an external stylesheet, simply declare the `external_stylesheets` attribute when initializing your Dash app:

    `app = dash(external_stylesheets=[\"URL TO YOUR STYLESHEET\"])`
    "),
    html_hr(),
    dcc_markdown("

    Q: How can I add JavaScript to my Dash app?
    A: To use an external script, simply declare the `external_scripts` attribute when initializing your Dash app:

    `app = dash(external_scripts=[\"URL TO YOUR SCRIPT\"])`

    ### Gotchas

    There are some aspects of how Dash works that can be counter-intuitive. This
    can be especially true of how the callback system works. This section outlines
    some common Dash gotchas that you might encounter as you start building out more
    complex Dash apps. If you have read through the rest of [Dash Tutorial](/) and are
    encountering unexpected behavior, this is a good section to read through. If you still have
    residual questions, the [Dash Community forums](https://community.plotly.com/c/dash/julia/20)
    is a good place to ask them.

    #### Callbacks require their `inputs`, `states`, and `output` to be present in the layout

    By default, Dash applies validation to your callbacks, which performs checks such as
    validating the types of callback arguments and checking to see whether the specified `input`
    and `output` components actually have the specified properties. For full validation, all components
    withiin your callback must thererfore appear in the initial layout of your app, and you will see an
    error if they do not.

    However, in the case of more complex Dash apps that involve dynamic modification of the layout
    (such as multipage apps), not every component appearing in your callbacks will be included in the
    initial layout. You can remove this restrictioin by disabling callback validation like this:

    `app = dash(suppress_callback_exceptions = true)`

    #### Callbacks require *all* `inputs`, `states`, and `output` to be rendered on the page

    If you have disabled callback validation in order to support dynamic layouts, then you won't be
    automatically alerted to the situation where a component within a callback is not found within a
    layout. In this situation, where a component registered with a callback is missing from the layout,
    the callback will fail to fire. For example, if you define a callback with only a subset of the
    specified `inputs` present in the current page layout, the callback will simply not fire at all.

    #### A component/property pair can only be the `output` of one callback

    For a given component/property pair, (eg `my-graph`, `figure`), it can only be registered as the
    `output` of one callback. If you want to associate two logically separate sets of `inputs` with
    the one output component/property pair, you'll have to bundle them up into a larger callback and
    detect which of the relevant `inputs` triggered the callback inside the function. For `html_button`
    elements, detecting which one triggered the callback can be done using the `n_clicks_timestamp` property.

    ### All callbacks must be defined before the server starts

    All your callbacks must be defined before your Dash app's server starts running, which is to say, before
    you call `run_server(app)`. This means that while you can assemble changed layout fragments dynamically
    during the handling of a callback, you can't define dynamic callbacks in response to the user input during
    the handling of a callback. If you have a dynamic interface, where a callback changes the layout to include
    a different set of input controls, then you must have already defined the callbacks required to service these
    new controls in advance.

    For example, a common scenario is a `dcc_dropdown` component that udpdates the current
    layout to replace a dashboard with another logically distinct dashboard that has a different set
    of controls (the number and type of which might depend on other user input) and different logic
    for generating the underlying data. A sensible organization would be for each of these dashboards
    to have separate callbacks. In this scenario, each of these callbacks must then be defined before the app starts running.

    Generally speaking, if a feature of your Dash app is that the number of `inputs` or `states` is determined
    by a user's input, then you must pre-define up front every permutation of callback that a user can potentially
    trigger.
    "),
    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")
end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
"""
    ZDT1(D, n_solutions)

ZDT1 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `D` number of variables (dimension)
- `n_solutions` number of pareto solutions.

Main properties:
- convex
"""
function ZDT1(D = 30, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        return ( [x[1], gx*(1 - sqrt(x[1] / gx))  ], [0.0], [0.0] )
    end

    bounds = Array([zeros(D) ones(D)]')

    x = range(0, 1, length=n_solutions)

    X = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ generateChild(x, f(x)) for x in X ]

    return f, bounds, pareto_set
end

"""
    ZDT2(D, n_solutions)

ZDT2 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `D` number of variables (dimension)
- `n_solutions` number of pareto solutions.

Main properties:
- nonconvex
"""
function ZDT2(D = 30, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        return ( [x[1], gx*(1 - (x[1] / gx)^2)  ], [0.0], [0.0] )
    end
    bounds = Array([zeros(D) ones(D)]')

    x = range(0, 1, length=n_solutions)

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ generateChild(x, f(x)) for x in x ]

    return f, bounds, pareto_set
end

"""
    ZDT3(D, n_solutions)

ZDT3 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `D` number of variables (dimension)
- `n_solutions` number of pareto solutions.

Main properties:
- convex disconected
"""
function ZDT3(D = 30, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        a = x[1] / gx
        return ( [x[1], gx*(1 - sqrt(a) - a*sin(10π*x[1]) )  ], [0.0], [0.0] )
    end
    bounds = Array([zeros(D) ones(D)]')

    regions = [ 0 0.0830015349
                0.182228780 0.2577623634
                0.4093136748 0.4538821041
                0.6183967944 0.6525117038
                0.8233317983 0.8518328654]

    n = n_solutions ÷ size(regions, 1)
    x = Float64[]
    for i in 1:size(regions, 1)
        x = vcat(x, range(regions[i,1], regions[i,2], length=n))
    end

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ generateChild(x, f(x)) for x in x ]

    return f, bounds, get_non_dominated_solutions(pareto_set)

end

"""
    ZDT4(D, n_solutions)

ZDT4 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `D` number of variables (dimension)
- `n_solutions` number of pareto solutions.

Main properties:
- nonconvex
"""
function ZDT4(D = 10, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 10*(length(x)-1) + sum( x[2:end].^2 - 10cos.(4π*x[2:end]))
        return ( [x[1], gx*(1 - sqrt(x[1] / gx))  ], [0.0], [0.0] )
    end
    bounds = Array([-5zeros(D) 5ones(D)]')
    bounds[:,1] = [0, 1.0]

    x = range(0, 1, length=n_solutions)

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ generateChild(x, f(x)) for x in x ]

    return f, bounds, pareto_set
end

"""
    ZDT6(D, n_solutions)

ZDT6 returns `(f::function, bounds::Matrix{Float64}, pareto_set::Array{xFgh_indiv})`
where `f` is the objective function and `pareto_set` is an array with optimal Pareto solutions
with `n_solutions`.

### Parameters
- `D` number of variables (dimension)
- `n_solutions` number of Pareto solutions.

Main properties:
- nonconvex
- non-uniformly spaced
"""
function ZDT6(D = 10, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1.0) )^(0.25)
        ff1 = 1.0 - exp(-4.0x[1])*sin(6.0π*x[1])^6
        return ( [ ff1 , gx*(1.0 - (ff1 / gx)^2)  ], [0.0], [0.0] )
    end

    bounds = Array([zeros(D) ones(D)]')

    #x = range(0, 1, length=n_solutions)

    #x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    xx = range(0.2807753191, 1, length=100)
    yy = 1 .- (xx).^2
    pareto_set = [ generateChild(zeros(0), ([xx[i], yy[i]], [0.0], [0.0])) for i in 1:length(xx) ]

    return f, bounds, pareto_set
end


text:
include("units.jl")
import Base.show

struct BuildOp
    unitBuildFunc
    unitType # needed for rendering
    buildLength::Number
    timer::GameTimer
    BuildOp(unitBuildFunc, unitType, buildLength) = new(unitBuildFunc, unitType, buildLength, GameTimer())
end
time_remaining(bop::BuildOp) = bop.buildLength - elapsed(bop.timer)
Base.show(io::IO, bop::BuildOp) = print(io, "BuildOp($(bop.unitBuildFunc)() in $(round(time_remaining(bop), 2)) secs)")

mutable struct Player
    units::PlayerUnits
    money
    build_ops::Array{BuildOp}
    Player() = new(PlayerUnits(), 0, BuildOp[])
    Player(m) = new(PlayerUnits(), m, BuildOp[])

    #Player() = (self = new(0, BuildOp[]); self.units = PlayerUnits(self))
    #Player(m) = (self = new(m, BuildOp[]); self.units = PlayerUnits(self))
end

function update!(p::Player, dt)
    p.money += length(p.units.collectors) * kMoneyPersecPercollector * dt

    for bop in p.build_ops
        update!(bop.timer, dt)
        if (time_remaining(bop) < 0)
            bop.unitBuildFunc()
        end
    end
    # Remove the completed build ops.
    filter!(bop->time_remaining(bop)>0, p.build_ops)
end

function purchase_unit!(p::Player, unitType::Type, addUnitFunc)
    if p.money < build_cost(unitType)
        return
    end
    p.money -= build_cost(unitType)
    bop = BuildOp(addUnitFunc, unitType, build_time(unitType))
    push!(p.build_ops, bop)
    start!(bop.timer)
    sort!(p.build_ops; by=time_remaining)
    return (bop, p.money)
end
function purchase_collector!(p::Player, pos = WorldPos(0,0))
     purchase_unit!(p, Collector, ()->add_unit!(p.units, Collector(p.units, pos)))
 end
function purchase_fighter!(p::Player, pos = WorldPos(0,0))
     purchase_unit!(p, Fighter, ()->add_unit!(p.units, Fighter(p.units, pos)))
 end

# Tests
p = Player(20)
purchase_collector!(p)
@test 17 == p.money
purchase_fighter!(p)
purchase_collector!(p)
@test 3 == length(p.build_ops)
update!(p, 0.05)
@test 3 == length(p.build_ops)
update!(p, 100)
@test 0 == length(p.build_ops)
@test 3 == length(p.units.units)
@test 2 == length(p.units.collectors)
# kill unit
try
    while true
        attack!(p.units.units[1], p.units.fighters[1])
    end
catch end
@test 0 == length(p.units.units)


text:
# parameters
path_to_working_folder = "data"
mass_bin_name = "0520_0540"#ARGS[1]
tslice = "t1"
tmin = "0.1"
tmax = "0.112853"
t_list = ["0.100000-0.112853","0.112853-0.127471","0.127471-0.144385","0.144385-0.164401","0.164401-0.188816"
                ,"0.188816-0.219907","0.219907-0.262177","0.262177-0.326380","0.326380-0.448588","0.448588-0.724294","0.724294-1.000000"]
######################################################
@show ARGS
# set names
for app in ["rd", "mc", "fu"]
    pwf = path_to_working_folder
    @eval $(Symbol("kinvar_"*app))    = joinpath($pwf,"variables_$(mass_bin_name)_$(tslice)_"*$app*".bin")
end
# Check that the files exist
#
name_of_the_file_rd = joinpath(path_to_working_folder,"$(mass_bin_name)_rd.root")
isfile(name_of_the_file_rd) || error("Cannot find the file $name_of_the_file_rd")
#
name_of_the_file_mc = joinpath(path_to_working_folder,"$(mass_bin_name)_mc.root")
isfile(name_of_the_file_mc) || error("Cannot find the file $name_of_the_file_mc")

# Converting data.root to data.txt
let
    @time run(`src/Extract data/$(mass_bin_name)_rd.root D $(tmin) $(tmax) $(kinvar_rd)`)
    @time run(`src/Extract data/$(mass_bin_name)_mc.root M $(tmin) $(tmax) $(kinvar_mc)`)
    @time run(`src/Extract data/$(mass_bin_name)_mc.root F $(tmin) $(tmax) $(kinvar_fu)`)
end


text:
struct ResNetModel
    entry
    pooling
    layers
    head

    size::Int64
end
Flux.@functor ResNetModel

(m::ResNetModel)(x) = x |> m.entry |> m.pooling |> m.layers |> m.head

function (m::ResNetModel)(x, ::Val{:stages})
    stages = typeof(x)[]

    o = x |> m.entry
    push!(stages, o)
    o = o |> m.pooling

    for l in m.layers
        o = o |> l
        push!(stages, o)
    end
    stages
end

const Config = Dict(
    18=>((2, 2, 2, 2), BasicBlock, 1),
    34=>((3, 4, 6, 3), BasicBlock, 1),
    50=>((3, 4, 6, 3), Bottleneck, 4),
    101=>((3, 4, 23, 3), Bottleneck, 4),
    152=>((3, 8, 36, 3), Bottleneck, 4),
)

function ResNetModel(
    model_size::Int64 = 18;
    in_channels::Int64 = 3,
    classes::Union{Int64, Nothing} = 1000,
)
    if !(model_size in keys(Config))
        throw(
            "Invalid mode size [$model_size]. " *
            "Supported sizes are $(keys(Config))."
        )
    end
    repeats, block, expansion = Config[model_size]

    entry = Chain(
        Conv((7, 7), in_channels=>64, pad=(3, 3), stride=(2, 2), bias=false),
        BatchNorm(64, relu),
    )
    pooling = MaxPool((3, 3), pad=(1, 1), stride=(2, 2))

    head = nothing
    if classes ≢ nothing
        head = Chain(
            MeanPool((7, 7)), flatten, Dense(512 * expansion, classes),
        )
    end

    in_channels = 64
    channels = (64, 128, 256, 512)
    strides = (1, 2, 2, 2)

    layers = []
    for (out_channels, repeat, stride) in zip(channels, repeats, strides)
        push!(layers, make_layer(
            block, in_channels=>out_channels, repeat, expansion, stride,
        ))
        in_channels = out_channels * expansion
    end

    ResNetModel(entry, pooling, Chain(layers...), head, model_size)
end

function stages_channels(r::ResNetModel)
    e = Config[r.size][3]
    (64, 64 * e, 128 * e, 256 * e, 512 * e)
end


text:
using IsDef
using Test
using Documenter
import FunctionWrappers: FunctionWrapper

@test isempty(detect_ambiguities(IsDef))

macro test_everything()
  quote

    # test apply
    # ----------

    @test apply(reduce, +, [2,3], init=1) == reduce(+, [2,3], init=1)

    # test ∨ \vee
    # -----------

    @test ∨(1, 2.0) == ∨(typeof(1), typeof(2.0)) == promote_type(typeof(1), typeof(2.0))

    # test isdef/Out
    # --------------

    @test IsDef.Out(Base.map, typeof(x->2x), Vector{Int}) == Vector{Int}
    @test isdef(Base.reduce, typeof(+), Vector{<:Number})

    # Caution!! Core.Compiler.return_type is actually not always as good as expected:
    # this by now actually inferes correctly...
    # Base.promote_op(Base.reduce, typeof(+), Vector{String}) === Union{}

    # works transparent with wrappers (unlike Base.`which`)
    wrapper(args...; kwargs...) = original(args...; kwargs...)
    original(a::Int, b::String) = true
    @test isdef(wrapper, Int, String)
    @test !isdef(wrapper, Float64)

    # does even work on compiler level
    @test Out(Out, typeof(Base.map), typeof(x->2x), Vector{Int}) == Type{Vector{Int}}

    # works in a strict open sense only with Any (would have to work for a whole newtype)
    # for everything else the concrete leave-types are used
    f(a) = a + a
    # we decided to leave Any as is, not going to newtype, and hence Any says often yes instead of no
    @test isdef(f, Any)
    @test isdef(f, Number)
    @test isdef(f, Integer)
    @test !isdef(f, String)
    @test !isdef(f, AbstractString)
    @test !isdef(f, Vector{<:String})
    @test !isdef(f, Vector{<:AbstractString})

    @test isdef(Base.map, typeof(x -> x+4), Array{<:Number, 3})
    @test !isdef(Base.map, typeof(x -> x+4), Vector{String})
    @test !isdef(Base.map, typeof(x -> x+4), Vector{AbstractString})
    @test !isdef(Base.map, typeof(x -> x+4), Vector{<:AbstractString})

    # test types
    @test Out(Some, Int) == Some{Int}
    @test isdef(Some, Int)

    # test values
    @test Out(sin, 1) == Float64
    @test isdef(sin, 1)

    # test inference
    # --------------

    @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(sin), Int) == Val{true}

    using IsDef
    using Test
    mywrapper(args...) = myfunc(args...)
    myfunc(::BigFloat) = "big"
    myfunc(::Float16) = 16
    myfunc(::Float32) = 32
    @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(mywrapper), Float64) == Val{false}
    @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(mywrapper), BigFloat) == Val{true}
    @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(mywrapper), Float16) == Val{true}
    @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(mywrapper), Float32) == Val{true}

    # TODO unfortunately type inference does not work through Union, which is surprising
    # inspecting the issue further it turns out that type-inference is extremely unstable, putting things into functions
    # instead of doing them at the REPL may influence type-inference to the worse...
    # better not to do much with type-inference

    # @test Base.promote_op((args...) -> Val(isdef(args...)), typeof(mywrapper), AbstractFloat) == Val{false}

    # test documentation
    # ------------------

    doctest(IsDef)
    @test Out(Base.map, typeof(string), Vector{Int}) == Vector{String}
    @test Out(Base.map, typeof(isodd), Vector{Int}) == Vector{Bool}
    @test Out(Base.map, FunctionWrapper{Bool, Tuple{Any}}, Vector{Int}) == Vector{Bool}

  end # quote
end # macro

@test_everything

# check that if we redefine the generated functions of the module, nothing breaks
# ----------------------------

IsDef.@redefine_generated
@test_everything