text:
#
# EasyFITS.jl -
#
# Utilities to simplify reading/writing FITS files.
#
#-------------------------------------------------------------------------------
#
# This file if part of EasyFITS software (https://github.com/emmt/EasyFITS.jl)
# licensed under the MIT license.
#
# Copyright (C) 2018-2020, Éric Thiébaut.
#

module EasyFITS

export
    FitsArray,
    FitsBitpix,
    FitsComment,
    FitsFile,
    FitsHDU,
    FitsHeader,
    FitsIO,
    FitsImage,
    FitsImageHDU,
    FitsTableHDU,
    exists,
    readfits,
    writefits,
    writefits!,
    write!

using FITSIO
using CFITSIO
using CFITSIO: libcfitsio

using Base: getfield, elsize, tail, OneTo,
    throw_boundserror, @propagate_inbounds
import Base: get, getindex, setindex!, keys, haskey, getkey,
    read, write, open, close

# Apply patches.
include("patches.jl")

#
# Predefine all types.
#

struct FitsIO
    fh::FITS
end

struct FitsBitpix{N} end

struct FitsHDU{T<:HDU}
    hdu::T
end

const FitsImageHDU{T<:ImageHDU} = FitsHDU{T}
const FitsTableHDU = FitsHDU{TableHDU}

struct FitsHeader
    hdr::FITSHeader
end

struct FitsImage{T,N} <: DenseArray{T,N}
    arr::Array{T,N}
    hdr::FitsHeader
end

# Singleton type for writing a FITS file.
struct FitsFile end

# Abstract type for retrieving the array part of an Image HDU in a FITS file.
abstract type FitsArray{T,N} <: AbstractArray{T,N} end

# Singleton type for retrieving a keyword comment.
struct FitsComment end

# Singleton type for marking a missing keyword.
struct Missing end

# This decoration is used to specialize `tryparse` for FITS keyword values.
struct FitsUnparsedValue
    str::String
end

# Annotated objects have a FITS header and implements indexation by keywords,
# they can also be directly read from a FITS file.
const Annotated = Union{FitsImage,FitsHeader}

# Readable objects are those which can be directly read from a FITS file.
const Readable = Union{Annotated,FitsArray}

# Allowed types for keywords values.
const KeywordValues = Union{Bool,Int,Float64,String}

# Valid types to specify a FITS extension.
const Extension = Union{Integer,AbstractString}

# Valid types to specify a FITS header.
const HeaderLike = Union{FitsHeader,NamedTuple,
                         Tuple{Vararg{Pair{<:AbstractString}}}}

readfits(filename::AbstractString; ext::Extension = 1) =
    read(FitsImage, filename; ext=ext)

readfits(::Type{Array}, filename::AbstractString; ext::Extension = 1) =
    read(FitsArray, filename; ext=ext)

function readfits(::Type{Array{T}}, filename::AbstractString;
                  ext::Extension = 1) where {T}
    read(FitsArray{T}, filename; ext=ext)
end

function readfits(::Type{Array{T,N}}, filename::AbstractString;
                  ext::Extension = 1) where {T,N}
    read(FitsArray{T,N}, filename; ext=ext)
end

writefits(filename::AbstractString, args...; kwds...) =
    write(FitsFile, filename, args...; kwds...)

writefits!(filename::AbstractString, args...; kwds...) =
    write(FitsFile, filename, args...; kwds...)

"""
    exists(path) -> boolean

yields whether file `path` exists.  Argument can be a file name or an instance
of `Base.Filesystem.StatStruct` as returned by the `stat` method.

"""
exists(st::Base.Filesystem.StatStruct) = (st.nlink != 0)
exists(path) = exists(stat(path))

"""
    EasyFITS.hduname(T) -> name, vers=1

yields the name and revision number of FITS Header Data Unit for storing data
of type `T`.  Other packages are encouraged to extend this method with their
types.

"""
hduname(::T) where {T} = hduname(T)
@noinline hduname(::Type{T}) where {T} =
    error(string("method `EasyFITS.hduname` has not been specialized ",
                 "for type `", T, "`"))

# Make sure that a (name, vers) tuple is returned.
fixhduname(x) = _fixhduname(hduname(x))
_fixhduname(str::AbstractString) = _fixhduname(str, 1)
_fixhduname(arg::Tuple{AbstractString,Integer}) = _fixhduname(arg[1], Int(arg[2]))
_fixhduname(arg::Tuple{AbstractString,Int}) = _fixhduname(String(arg[1]), arg[2])
_fixhduname(arg::Tuple{String,Int}) = arg

"""
    FitsIO(path, mode="r") -> io

opens of creates a FITS file whose name is `path` and returns an instance of
`FitsIO` which is an opaque handle to a FITS file.  Argument `mode` can be:

- `"r"` or `"r+"` to read or append to the contents of and existing FITS file.
  If file `path` does not exist but `path` does not end with the `".gz"`
  extension and `"\$path.gz"` does exist, then the compressed file
  `"\$path.gz"` is open instead.

- `"w"` to create a new FITS file named `path` that must not already exist.  An
  error is thrown if the file already exists.

- `"w!"` to open FITS file named `path` for writing.  If the file already
  exists, it is (silently) overwritten.

Call `close(io)` to close the FITS file associated with the `FitsIO` instance
`io`.  Call `isopen(io)` to check whether the FITS file associated with the
`FitsIO` instance `io` is open.  Closing the FITS file is automatically done,
if needed, when the instance is garbage collected.

The do-block syntax is supported to automatically close the FITS file:

    FitsIO(filename, mode="r") do io
        # use FITS handle io
        ...
    end

An instance of `FitsIO` is a collection of *Header Data Units* (HDU) and
implements indexation and iteration.  Assuming `io` is a `FitsIO` object, then:

- `io[i]` yields the `i`-th HDU.

- `length(io)` yields the number of HDUs.

- `io[name]` or `io[name,vers]` yields the HDU whose `EXTNAME` (or `HDUNAME`)
  keyword is equal to `name` (a string) and, optionally, whose `EXTVER` (or
  `HDUVER`) keyword is equal to `vers` (an integer).

- You can do `for hdu in io; ...; end` to iterate through all HDU's of `io`.

- Methods `findfirst(p,io)`, `findlast(p,io)`, `findnext(p,io,i)` and
  `findprev(p,io,i)` can be used on `FitsIO` object `io` to search for a
  specific HDU.  These methods test each HDU (starting at initial index `i` for
  `findnext` and `findprev`) with the predicate function `p` (called with a
  `FitsHDU` argument) and return the index (an integer) of the first HDU for
  which the predicate yields `true` or `nothing` if this never occurs.

Also see: [`EasyFITS.find`](@ref).

""" FitsIO

FitsIO(path::AbstractString, mode::AbstractString="r") = begin
    local handle::FITS
    if mode == "r" || mode == "r+"
        handle = FITS((isfile(path) || endswith(path, ".gz") ||
                       !isfile(path*".gz") ? path : path*".gz"), mode)
    else
        if mode == "w"
            if exists(path)
                throw_file_already_exists(path, "use mode \"w!\" to overwrite")
            end
        elseif mode != "w!"
            throw(ArgumentError("mode must be \"r\", \"r+\", \"w\" or \"w!\""))
        end
        handle = FITS(path, "w")
    end
    return FitsIO(handle)
end

FitsIO(func::Function, path::AbstractString, mode::AbstractString="r") = begin
    io = FitsIO(path, mode)
    try
        func(io)
    finally
        close(io)
    end
end

Base.get(::Type{FITS}, obj::FitsIO) = getfield(obj, :fh)
Base.get(::Type{FitsIO}, obj::FitsIO) = obj

Base.convert(::Type{FitsIO}, obj::FITS) = FitsIO(obj)
Base.convert(::Type{FITS}, obj::FitsIO) = get(FITS, obj)

Base.length(io::FitsIO) = length(get(FITS, io))
Base.getindex(io::FitsIO, args...) = FitsHDU(getindex(get(FITS, io), args...))
Base.iterate(io::FitsIO, args...) = begin
    result = iterate(get(FITS, io), args...)
    if result === nothing
        return nothing
    end
    return (FitsHDU(result[1]), result[2])
end
Base.lastindex(io::FitsIO) = lastindex(get(FITS, io))
Base.show(io::IO, obj::FitsIO) = show(io, get(FITS, obj))

Base.findfirst(pred::Function, io::FitsIO) = findnext(pred, io, 1)
Base.findlast(pred::Function, io::FitsIO) = findprev(pred, io, length(io))
Base.findnext(pred::Function, io::FitsIO, i::Integer) = findnext(pred, io, Int(i))
Base.findnext(pred::Function, io::FitsIO, i::Int=1) = find(pred, io, i:length(io))
Base.findprev(pred::Function, io::FitsIO, i::Integer) = findprev(pred, io, Int(i))
Base.findprev(pred::Function, io::FitsIO, i::Int=1) = find(pred, io, i:-1:1)

"""
    EasyFITS.find(pred, io, I=1:length(io))

yields the first index HDU `i ∈ I` of FITS instance `io` for which the
predicate `pred(io[i])` returns `true`, or [`nothing`](@ref) if `pred(io[i])`
returns `false` for all `i ∈ I`.

For instance:

    i = EasyFITS.find(hdu -> get(String,hdu,"EXTNAME","") == "CALIBRATION", io)
    if i === nothing
        # not found
        ...
    else
        # found at index i
        hdu = io[i]
    end

"""
function find(pred::Function, io::FitsIO,
              I::AbstractVector{Int} = 1:length(io)) :: Union{Nothing,Int}
    for i in I
        if pred(io[i])
            return i
        end
    end
    return nothing
end

"""

An instance of type `FitsHDU{T}` is an opaque handle to a FITS *Header Data
Unit* (HDU).  Parameter `T` is to distinguish between the different kinds of
HDU (*Image* or *Binary Table*).

""" FitsHDU

Base.get(::Type{HDU}, obj::FitsHDU) = getfield(obj, :hdu)
Base.get(::Type{ImageHDU}, obj::FitsImageHDU) = get(HDU, obj)
Base.get(::Type{TableHDU}, obj::FitsTableHDU) = get(HDU, obj)
Base.get(::Type{T}, obj::T) where {T<:FitsHDU} = obj

Base.convert(::Type{FitsHDU}, obj::HDU) = FitsHDU(obj)
Base.convert(::Type{FitsImageHDU}, obj::ImageHDU) = FitsHDU(obj)
Base.convert(::Type{FitsTableHDU}, obj::TableHDU) = FitsHDU(obj)

Base.convert(::Type{HDU}, obj::FitsHDU) = get(HDU, obj)
Base.convert(::Type{ImageHDU}, obj::FitsImageHDU) = get(HDU, obj)
Base.convert(::Type{TableHDU}, obj::FitsTableHDU) = get(HDU, obj)

Base.show(io::IO, obj::FitsHDU) = show(io, get(HDU, obj))
Base.size(obj::FitsImageHDU, args...) = size(get(HDU, obj), args...)
Base.ndims(obj::FitsImageHDU) = ndims(get(HDU, obj))
Base.length(obj::FitsImageHDU) = length(get(HDU, obj))
Base.lastindex(obj::FitsImageHDU) = lastindex(get(HDU, obj))

"""
    isprimary(hdu)

yields whether `hdu` is a primary FIT header data units (HDU).

"""
isprimary(obj::FitsHDU) = isprimary(get(HDU, obj))
isprimary(obj::HDU) = (getfield(obj, :ext) == 1)

"""

An object of type `FitsHeader` stores the header part of a FITS *Header Data
Unit* (HDU).  There are several ways to build an instance of `FitsHeader`.

    FitsHeader()

yields an empty instance of `FitsHeader`.

    FitsHeader(io, ext=1)

yields an instance of `FitsHeader` built over FITS header read from extension
`ext` in FITS file `io`.

    FitsHeader(; key1=val1, key2=val2, ...)

yields an instance of `FitsHeader` whose contents is set by keywords (values
can be a tuple with a value and a comment).  This is can also be done by
specifying key-value pairs:

    FitsHeader("key1" => val1, "key2" => val2, ...)

To avoid ambiguities the two styles cannot be mixed.

An object of type `FitsHeader`, say `obj`, implements indexation by keywords,
as with `obj[key]`, and `obj.key` syntax.

## Implementation

Type `FitsHeader` is a simple wrapper over `FITSIO.FITSHeader` (beware of the
different spellings) to implement indexation by keywords and `obj.key` syntax.
We have to define a new type for that because implementing the extended
interface over the existing `FITSIO.FITSHeader` would constitute a *type
piracy*.  An instance of `FitsHeader` can be simply built over an object `hdr`
of type `FITSIO.FITSHeader` by calling:

    FitsHeader(hdr)

""" FitsHeader

FitsHeader(hdu::FitsHDU) = FitsHeader(read_header(get(HDU, hdu)))

FitsHeader(io::FitsIO, ext::Extension = 1) = FitsHeader(io[ext])

function FitsHeader(; kwds...)
    hdr = FitsHeader(FITSHeader(String[], [], String[]))
    for (key, val) in kwds
        setproperty!(hdr, key, val)
    end
    return hdr
end

function FitsHeader(args::Pair{<:AbstractString}...)
    hdr = FitsHeader(FITSHeader(String[], [], String[]))
    for (key, val) in args
        hdr[key] = val
    end
    return hdr
end

Base.get(::Type{Array}, obj::FitsHeader) = nothing # FIXME: ???
Base.get(::Type{FitsHeader}, obj::FitsHeader) = obj
Base.get(::Type{FITSHeader}, obj::FitsHeader) = getfield(obj, :hdr)

Base.convert(::Type{FitsHeader}, hdr::FITSHeader) = FitsHeader(hdr)
Base.convert(::Type{FITSHeader}, hdr::FitsHeader) = get(FITSHeader, hdr)

Base.length(obj::FitsHeader) = length(get(FITSHeader, obj))
Base.show(io::IO, obj::FitsHeader) = show(io, get(FITSHeader, obj))

"""

Type `FitsComment` is a singleton used as a marker to indicate that the comment
of a FITS keyword is to be returned by the `get` method.

""" FitsComment

"""
    FitsBitpix(arg)

yields a singleton type which encapsulates the FITS *bitpix* (for
*bits-per-pixel*) code identifying array element type according to `arg`.
Argument `arg` can be a Julia type, an integer (interpreted as a FITS bitpix
code), an array, a FITS HDU/image/header.

Call `eltype(bpx)` to convert FITS bitpix `bpx` into a Julia type.

"""
FitsBitpix(n::Integer) = FitsBitpix{Int(n)}()
FitsBitpix(::Type{T}) where {T} = FitsBitpix(bitpix_from_type(T))
FitsBitpix(A::AbstractArray) = FitsBitpix(typeof(A))
FitsBitpix(::Type{<:AbstractArray{T}}) where {T} = FitsBitpix(T)
FitsBitpix(hdr::Union{<:FitsHeader,<:FITSHeader}) = FitsBitpix(hdr["BITPIX"])
FitsBitpix(hdu::Union{<:FitsImageHDU,<:ImageHDU}) =
    FitsBitpix(fits_get_img_equivtype(getfile(hdu)))

Base.eltype(::FitsBitpix{N}) where {N} = type_from_bitpix(Val(Cint(N)))

# FIXME: Base.eltype could also be extended to FITSHeader and ImageHDU.
Base.eltype(hdr::FitsHeader) =
    type_from_bitpix(Val(Cint(hdr["BITPIX"])))
Base.eltype(hdu::FitsImageHDU) =
    type_from_bitpix(Val(Cint(fits_get_img_equivtype(getfile(hdu)))))

"""
    FitsImage(arr, hdr)

yields a FITS Image that behaves like a Julia array when indexed by integers of
Cartesian indices and like a dictionary of FITS keywords when indexed by
strings.  Argument `arr` specifies the array contents of the object and `hdr`
specifies the keywords both are shared by the returned object.

The initial contents of the header may be specified by keywords.  For instance:

    FitsImage(arr; KEY1 = val1, KEY2 = (val2, com2), ...)

The embedded array may be creted by the constructor:

    FitsImage{T}(init, dims...; kwds...)

yields a FITS Image with an empty header and array-like data of dimensions
`dims` and with elements of type `T` initially set by `init`.  Specifying the
number of dimensions is also supported:

    FitsImage{T,N}(init, dims...; kwds...)

""" FitsImage

FitsImage(arr::AbstractArray; kwds...) = FitsImage(arr, FitsHeader(; kwds...))
FitsImage(arr::AbstractArray{T,N}, hdr::FitsHeader) where {T,N} =
    FitsImage{T,N}(convert(Array{T,N}, arr), hdr)
FitsImage(arr::AbstractArray, hdr::FITSHeader) =
    FitsImage(arr, FitsHeader(hdr))
FitsImage{T}(init, dims::Integer...; kwds...) where {T} =
    FitsImage{T}(init, dims; kwds...)
FitsImage{T}(init, dims::NTuple{N,Integer}; kwds...) where {T,N} =
    FitsImage{T,N}(init, dims; kwds...)
FitsImage{T,N}(init, dims::Integer...; kwds...) where {T,N} =
    FitsImage{T,N}(init, dims; kwds...)
FitsImage{T,N}(init, dims::NTuple{N,Integer}; kwds...) where {T,N} =
    FitsImage{T,N}(Array{T,N}(init, dims), FitsHeader(; kwds...))

Base.get(::Type{Array}, obj::FitsImage) = getfield(obj, :arr)
Base.get(::Type{FitsHeader}, obj::FitsImage) = getfield(obj, :hdr)
Base.get(::Type{FITSHeader}, obj::FitsImage) = get(FITSHeader, get(FitsHeader, obj))
Base.get(::Type{FitsImage}, obj::FitsImage) = obj

Base.convert(::Type{FitsImage}, arr::AbstractArray) = FitsImage(arr)
Base.convert(::Type{FitsHeader}, obj::FitsImage) = get(FitsHeader, obj)
Base.convert(::Type{FITSHeader}, obj::FitsImage) = get(FITSHeader, obj)
Base.convert(::Type{Array}, obj::FitsImage) = get(Array, obj)
Base.convert(T::Type{<:AbstractArray}, obj::FitsImage) =
    convert(T, get(Array, obj))

# Make FitsImage instances behave like arrays (indexing is considered later).
#Base.eltype(::FitsImage{T,N}) where {T,N} = T # FIXME: not needed
#Base.ndims(::FitsImage{T,N}) where {T,N} = N # FIXME: not needed
Base.length(A::FitsImage) = length(get(Array, A))
Base.size(A::FitsImage) = size(get(Array, A))
Base.size(A::FitsImage, d) = size(get(Array, A), d)
Base.axes(A::FitsImage) = axes(get(Array, A))
Base.axes(A::FitsImage, d) = axes(get(Array, A), d)
@inline Base.axes1(A::FitsImage) = axes1(get(Array, A))
Base.IndexStyle(::Type{<:FitsImage}) = IndexLinear()
Base.similar(A::FitsImage, T::Type=eltype(A), dims::NTuple{N,Int}=size(A)) where {N} =
    FitsImage{T,N}(undef, dims)
Base.elsize(::Type{FitsImage{T,N}}) where {T,N} = elsize(Array{T,N})
Base.sizeof(A::FitsImage) = sizeof(get(Array, A))

# For AbstractArray types, a more specific version of `Base.show` must be extended.
Base.show(io::IO, ::MIME"text/plain", A::FitsImage) = show(io, A)
function Base.show(io::IO, A::FitsImage{T,N}) where {T,N}
    println(io, "FitsImage{$T,$N} of size $(size(A)) and keywords:")
    show(io, FitsHeader(A))
end

# Make FitsImage's efficient iterators.
@inline Base.iterate(A::FitsImage, i=1) =
    ((i % UInt) - 1 < length(A) ? (@inbounds A[i], i + 1) : nothing)

# FIXME: copyto!, convert, unsafe_convert, pointer, etc.

# Maximum size of the value field including the final '\0'.
const VALUE_SIZE = 72

"""
    findkeys(obj, keys [, buf])

finds a FITS card in `obj` with one of the keywords in `keys` and returns a
string with the unparsed value of the keyword if found and `nothing` otherwise.

Optional argument `buf` is a small scratch buffer of bytes (`UInt8`) used to
temporarily store the unparsed value before it is converted into a Julia
string, its size is augmented if it is smaller than `VALUE_SIZE`.

"""
function findkeys(obj::Union{FitsIO,FitsHDU},
                  keys::NTuple{N,AbstractString},
                  buf::Vector{UInt8}=Vector{UInt8}(undef, VALUE_SIZE)) where {N}
    file = getfile(obj)
    status = Ref{Cint}()
    length(buf) ≥ VALUE_SIZE || resize!(buf, VALUE_SIZE)
    for key in keys
        status[] = 0
        ccall((:ffgkey, libcfitsio), Cint,
              (Ptr{Cvoid},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{Cint}),
              file.ptr, key, buf, C_NULL, status)

        # If the key is found, return it. If there was some other error
        # besides key not found, throw an error.
        if status[] == 0
            return unsafe_string(pointer(buf))
        elseif status[] != 202
            error(FITSIO.fits_get_errstatus(status[1]))
        end
    end
    return nothing
end

@doc @doc(trygetkeys) VALUE_SIZE

"""
    get([T,] obj, key[, def])

yields the value of the FITS keyword `key` in object `obj`.  Argument `T` may
be used to specify the type of the result (i.e. `Int`, `Bool`, `Float64` or
`String`) or `FitsComment` to retrieve the comment associated with the FITS
keyword `key`.  Argument `def` may be used to specify the value to return when
the keyword is not found (`def` is not converted to type `T` if this argument
is specified).  A `KeyError` exception is thrown if the keyword is not found
and no default value `def` is specifed.  An error is thrown if the keyword is
found but its value cannot be converted to `T` if it is specifed.

The syntax `get(T,obj)` is also extended to retrieve various things from object
`obj` of type `FitsHeader`, `FitsHDU` or `FitsImage`.  Examples:

```julia
get(Array, obj)            -> arr # get array associated to a `FitsImage`
get(FitsHeader, obj)       -> hdr # get FITS header of object `obj`
get(FitsComment, obj, key) -> str # yields comment for FITS keyword `key` in object `obj`
```

"""
get(obj::Annotated, key::AbstractString, def) = begin
    hdr = get(FITSHeader, obj)
    return (haskey(hdr, key) ? getindex(hdr, key) : def)
end

get(obj::Annotated, key::AbstractString) = begin
    hdr = get(FITSHeader, obj)
    haskey(hdr, key) || missing_keyword(key)
    return getindex(hdr, key)
end

get(::Type{FitsComment}, obj::Annotated, key::AbstractString, def) = begin
    hdr = get(FITSHeader, obj)
    return (haskey(hdr, key) ? FITSIO.get_comment(hdr, key) : def)
end

get(::Type{FitsComment}, obj::Annotated, key::AbstractString) = begin
    hdr = get(FITSHeader, obj)
    haskey(hdr, key) || missing_keyword(key)
    return get_comment(hdr, key)
end

for S in (AbstractFloat, Integer, AbstractString, Bool)
    @eval begin
        function get(::Type{T}, obj::Annotated,
                          key::AbstractString, def) where {T<:$S}
            hdr = get(FITSHeader, obj)
            return (haskey(hdr, key) ?
                    checkvalue(T, $S, getindex(hdr, key), key) : def)
        end

        function get(::Type{T}, obj::Annotated,
                          key::AbstractString) where {T<:$S}
            hdr = get(FITSHeader, obj)
            haskey(hdr, key) || missing_keyword(key)
            return checkvalue(T, $S, getindex(hdr, key), key)
        end
    end
end

function get(::Type{T}, obj::Union{FitsIO,FitsHDU}, key::AbstractString,
             def = Missing()) where {T<:KeywordValues}
    raw = findkeys(obj, (key,))
    if raw === nothing
        def === Missing() && missing_keyword(key)
        return def
    end
    if iscomment(key)
        if T == String
            return raw
        end
    else
        val = tryparse(T, FitsUnparsedValue(raw))
        if isa(val, T)
            return val
        end
    end
    error("value of FITS keyword \"", key,
          "\" cannot be converted to type `", T, "`")
end

function get(obj::Union{FitsIO,FitsHDU}, key::AbstractString, def = Missing())
    raw = findkeys(obj, (key,))
    if raw === nothing
        def === Missing() && missing_keyword(key)
        return def
    end
    if iscomment(key)
        return raw
    end
    c = first(raw)
    if c == 'T' || c == 'F'
        if length(raw) == 1
            return (c == 'T')
        end
    elseif c == Char(39) # single quote, using '\'' break Emacs fontify
        sval = tryparse(String, FitsUnparsedValue(raw))
        if isa(sval, String)
            return sval
        end
    else
        ival = tryparse(Int, FitsUnparsedValue(raw))
        if isa(ival, Int)
            return ival
        end
        rval = tryparse(Float64, FitsUnparsedValue(raw))
        if isa(rval, Float64)
            return rval
        end
    end
    error("unknown type of value for FITS keyword \"", key, "\"")
end

iscomment(key::AbstractString) =
    (key == "HISTORY" || key == "COMMENT")

Base.tryparse(::Type{T}, val::FitsUnparsedValue) where {T<:Real} =
    tryparse(T, val.str)

function Base.tryparse(::Type{Bool}, val::FitsUnparsedValue)
    if length(val.str) == 1
        if val.str[1] == 'T'
            return true
        elseif val.str[1] == 'F'
            return false
        end
    end
    return nothing
end

function Base.tryparse(::Type{String}, val::FitsUnparsedValue)
    # Check that the string in enclosed by single quotes.
    QUOTE = Char(39) # a single quote, avoing breaking Emacs' indentation
    str = val.str
    i = firstindex(str)
    j = lastindex(str)
    if length(str) < 2 || str[i] != QUOTE || str[j] != QUOTE
        return nothing
    end

    # Replace pairs of quotes by a single quote and discard trailing spaces.
    buf = Vector{UInt8}(undef, length(str) + 1)
    spaces = 0  # number of unwritten trailing spaces
    esc = false # previous character was a quote?
    n = 0
    for k in (i+1):(j-1)
        c = str[k]
        if esc
            if c != QUOTE
                return nothing
            end
            esc = false
        elseif c == QUOTE
            esc = true
            continue
        end
        if c == ' ' # FIXME: use `isspace`?
            spaces += 1
        else
            # Write any unwritten spaces and then the new non-space character.
            while spaces > 0
                spaces -= 1
                n += 1
                buf[n] = ' '
            end
            n += 1
            buf[n] = c
        end
    end
    if esc
        return nothing
    end
    buf[n+1] = 0
    return unsafe_string(pointer(buf))
end

#
# Override `getindex` and `setindex!` for efficient indexation by array
# indices.
#

@inline getindex(A::FitsImage, i::Int) = begin
    @boundscheck checkbounds(A, i)
    @inbounds val = get(Array, A)[i]
    return val
end

@inline setindex!(A::FitsImage, val, i::Int) = begin
    @boundscheck checkbounds(A, i)
    @inbounds get(Array, A)[i] = val
    return A
end

@inline Base.checkbounds(::Type{Bool}, A::FitsImage, i::Int) =
    1 ≤ i ≤ length(A)

#
# Override `getindex` and `setindex!` to implement `obj[key]`, `obj[key] = val`
# and `obj[key] = (val, com)` syntaxes.
#

@inline getindex(obj::Annotated, key::AbstractString) =
   getindex(get(FITSHeader, obj), key)

setindex!(obj::Annotated, val, key::AbstractString) = begin
    setindex!(get(FITSHeader, obj), val, key)
    return obj
end

setindex!(obj::Annotated, val::Tuple{Any,AbstractString}, key::AbstractString) = begin
    hdr = get(FITSHeader, obj)
    setindex!(hdr, val[1], key)
    set_comment!(hdr, key, val[2])
    return obj
end

setindex!(obj::Annotated, val::Tuple{Any,Nothing}, key::AbstractString) =
    setindex!(obj, (val[1], ""), key)

setindex!(obj::Annotated, val::Tuple{Any}, key::AbstractString) =
    setindex!(obj, val[1], key)

# Capture other calls as an error.
setindex!(obj::Annotated, val::Tuple, key::AbstractString) =
    throw(ArgumentError(string("unexpected value type `", typeof(val),
                               "` for FITS keyword")))

#
# Override `getproperty` and `setproperty!` to implement `obj.key` syntax.
#

@inline Base.getproperty(obj::T, sym::Symbol) where {T<:Annotated} =
    getindex(obj, propertyname(T, sym))

@inline Base.setproperty!(obj::T, sym::Symbol, val) where {T<:Annotated} =
    setindex!(obj, val, propertyname(T, sym))

"""
    EasyFITS.propertyname(T, sym)

converts symbol `sym` to a suitable key for object of type `T` throwing an
error if this conversion is not supported.

"""
@inline propertyname(::Type{<:Annotated}, sym::Symbol) = String(sym)
@noinline propertyname(::Type{T}, sym::Symbol) where {T} =
    error(string("Converting symbolic key for objects of type `", T, "` ",
                 "is not implemented. As a result, syntax `obj.key` is not ",
                 "supported for such objects."))

keys(obj::Annotated) = keys(get(FITSHeader, obj))
nkeys(obj::Annotated) = length(get(FITSHeader, obj))
nkeys(obj::Union{FITSHeader,AbstractDict}) = length(obj)
haskey(obj::Annotated, key) = haskey(get(FITSHeader, obj), key)
getkey(obj::Annotated, key, def) = getkey(get(FITSHeader, obj), key, def)

Base.delete!(obj::Annotated, i::Integer) = delete!(obj, Int(i))
Base.delete!(obj::Annotated, key::String) = delete!(obj, first_card(obj, key))
Base.delete!(obj::Annotated, idx::Int) = begin
    delete_card!(get(FITSHeader, obj), idx)
    return obj
end

Base.pop!(obj::Annotated, idx::Integer) = pop!(obj, Int(idx))
Base.pop!(obj::Annotated, idx::Integer, def) = pop!(obj, Int(idx), def)
Base.pop!(obj::Annotated, key::String) = pop!(obj, first_card(obj, key))
Base.pop!(obj::Annotated, key::String, def) = begin
    idx = first_card(obj, key, -1)
    return (idx == 0 ? def : pop!(obj, idx))
end
Base.pop!(obj::Annotated, idx::Int) = begin
    1 ≤ idx ≤ nkeys(obj) || error("out of bound FITS card index $idx")
    return unsafe_pop!(obj, idx)
end
Base.pop!(obj::Annotated, idx::Int, def) = begin
    1 ≤ idx ≤ nkeys(obj) || return def
    return unsafe_pop!(obj, idx)
end

unsafe_pop!(obj::Annotated, idx::Int) = begin
    hdr = get(FITSHeader, obj)
    val = hdr.values[idx]
    delete_card!(hdr, idx)
    return val
end

# Delete FITS card in header.  NOTE: Do not extend `delete!` to avoid
# type-piracy.  FIXME: Deleting a FITS card is expensive.
delete_card!(hdr::FITSHeader, i::Int) = begin
    if 1 ≤ i ≤ length(hdr)
        # Remove the entry from hdr.keys, hdr.values and hdr.comments and
        # rebuild the keyword map.
        delete_entry!(hdr.keys,     i)
        delete_entry!(hdr.values,   i)
        delete_entry!(hdr.comments, i)
        empty!(hdr.map)
        for j in 1:length(hdr.keys)
            key = hdr.keys[j]
            if !haskey(hdr.map, key)
                hdr.map[key] = j
            end
        end
    end
    nothing
end

# Delete an entry in a vector.
delete_entry!(A::Vector, i::Int) = begin
    n = length(A)
    if 1 ≤ i ≤ n
        for j in i:n-1
            A[j] = A[j+1]
        end
        resize!(A, n - 1)
    end
    nothing
end

# Yields index of first FITS card with given keyword.
first_card(obj::Annotated, key::String, def) =
    get(get(FITSHeader, obj).map, key, def)
first_card(obj::Annotated, key::String) = begin
    idx = first_card(obj, key, 0)
    idx == 0 && throw(KeyError(key))
    return idx
end

function checkvalue(::Type{T}, ::Type{S}, val,
                    key::AbstractString)::T where {S,T}
    isa(val, S) || bad_type(key)
    return T(val)
end

@noinline missing_keyword(key::AbstractString) = throw(KeyError(key))

@noinline bad_type(key::AbstractString) =
    error("bad type for FITS keyword \"$key\"")

"""
    read(::Type{T}, src) -> A

read an object of type `T<:Union{FitsImage,FitsHeader,FitsArray}` from FITS
extension `ext` in source `src` which can be the name of a FITS file or an
instance of `FitsIO` or `FitsHDU`.

If `src` is an instance of `FitsIO` or the name of a FITS file, keyword `ext`
may be used to specify the number or the name of the HDU number to read (the
first one by default).

Argument `T` specifies the type of the result.  Use `T<:Array` or
`T<:FitsArray` to retrieve the data part of the FITS HDU as a regular array.
Use `T<:FitsImage` to get a pseudo-array combining the header and data parts of
the FITS HDU.  Finally, use `T<:FitsHeader` to only retrieve the header part.

For `T<:Array` `T<:FitsArray`, or `T<:FitsImage`, the element type and number
of dimensions of the result can be specified, e.g. with `FitsImage{Float32,2}`

To avoid *type-piracy*, `T<:FitsArray` must be used instead of `T<:Array` when
`src` is a file name.

The result is indexable.  For `T<:FitsImage` and `T<:FitsHeader`, using a
string index yields the value of the corresponding FITS keyword in the header
part of the HDU.  Any other indices are used to access the contents of data
part of the HDU (as a regular Julia array).

Examples:

```julia
using EasyFITS
A = read(FitsImage, "image.fits")  # load the first HDU
A[2,3]                             # get value of data at indices (2,3)
A["BITPIX"]                        # get FITS bits per pixel
A.BITPIX                           # idem
getfitscomment(A, "BITPIX")        # get the associated comment
A["STUFF"] = 1                     # set value of FITS keyword STUFF
A["STUFF"] = (1, "Blah")           # idem with a value-comment pair
A.STUFF = 2                        # set value of FITS keyword STUFF
A.STUFF = (3, "Blah")              # idem with a comment
arr = get(Array, A)                # get the data part (a regular Julia array)
hdr = get(FitsHeader, A)           # get the header part
EasyFITS.nkeys(A)                  # get the number of keywords
EasyFITS.nkeys(hdr)                # get the number of keywords
length(hdr)                        # get the number of keywords
keys(A)                            # get the list of keywords
keys(hdr)                          # get the list of keywords
```

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

"""
read(::Type{T}, path::AbstractString, args...; kwds...) where {T<:Readable} =
    FitsIO(path, "r") do io
        return read(T, io, args...; kwds...)
    end

function read(::Type{T}, io::FitsIO, args...;
              ext::Extension = 1) where {T<:Union{Readable,Array}}
    read(T, io[ext], args...)
end

# Read header from FITS ImageHDU.
read(::Type{FitsHeader}, hdu::FitsHDU) =
    FitsHeader(read_header(get(HDU,hdu)))

# Read header and data from FITS Image HDU.
read(::Type{FitsImage{T,N}}, hdu::FitsHDU) where {T,N} =
    FitsImage(read(Array{T,N}, hdu), read(FitsHeader, hdu))
read(::Type{FitsImage}, hdu::FitsHDU) =
    FitsImage(read(Array, hdu), read(FitsHeader, hdu))
read(::Type{FitsImage{T}}, hdu::FitsHDU) where {T} =
    FitsImage(read(Array{T}, hdu), read(FitsHeader, hdu))

# Read array from FITS Image HDU.
read(hdu::FitsImageHDU, args...) = read(get(HDU, hdu), args...)
read(::Type{Array}, hdu::FitsImageHDU, args...) = read(hdu, args...)
read(::Type{Array{T}}, hdu::FitsImageHDU, args...) where {T} =
    convert(Array{T}, read(hdu, args...))
read(::Type{Array{T,N}}, hdu::FitsImageHDU, args...) where {T,N} =
    convert(Array{T,N}, read(hdu, args...))

read(::Type{FitsArray}, hdu::FitsHDU, args...) =
    read(Array, hdu, args...)
read(::Type{FitsArray{T}}, hdu::FitsHDU, args...) where {T} =
    read(Array{T}, hdu, args...)
read(::Type{FitsArray{T,N}}, hdu::FitsHDU, args...) where {T,N} =
    read(Array{T,N}, hdu, args...)

function Base.read!(hdu::FitsHDU, dst, args...)
    read!(get(HDU, hdu), dst, args...)
    return dst
end

"""
    write(FitsFile, path, args...; overwrite=false, kwds...)

creates a new FITS file `path` with contents built from the provided arguments
`args...` and keywords `kwds...`.  Unless `overwrite` is `true`, the file must
not already exist.  Instead of setting the `overwrite` keyword, simply call:

    write!(FitsFile, path, args...; kwds...)

to silently overwrites file `path` if it already exits.

Typical examples are:

    write(FitsFile, path, arr; KEY1 = val1, KEY2 = (val2, com2), ...)
    write(FitsFile, path, hdr, arr)
    write(FitsFile, path, (hdr1, arr1), arr2, (hdr3, arr3), ...)

with `arr*` arrays, `KEY*` keywords, `val*` keyword values, `com*` keyword
comments and `hdr*` objects of type `FitsHeader`.  In the last example, a FITS
file with 3 (or more) HDU's is created.

"""
function write!(::Type{FitsFile}, path::AbstractString, args...; kwds...)
    FitsIO(path, "w!") do io
        write(io, args...; kwds...)
    end
    nothing
end

function write(::Type{FitsFile}, path::AbstractString, args...;
               overwrite::Bool=false, kwds...)
    (overwrite == false && exists(path)) &&
        throw_file_already_exists(path, "try with `overwrite=true`")
    FitsIO(path, (overwrite ? "w!" : "w")) do io
        write(io, args...; kwds...)
    end
    nothing
end

write!(path::AbstractString, A::FitsImage) =
    write!(FitsFile, path, A)

write(path::AbstractString, A::FitsImage; overwrite::Bool=false) =
    write(FitsFile, path, A; overwrite=overwrite)

"""
    write(io::FitsIO, args...; kwds...)

writes data `args...` in FITS output `io` with keywords `kwds...`.  This method
can be specialized by foreign packages based on the types of the arguments
`args...`.

"""
write(io::FitsIO, arr::AbstractArray{<:Real}; kwds...) =
    write(io, FitsHeader(; kwds...), arr)
write(io::FitsIO, arr::AbstractArray{<:Real}, args::Pair{<:AbstractString}...) =
    write(io, FitsHeader(args...), arr)
write(io::FitsIO, hdr::FitsHeader, arr::AbstractArray{<:Real}) =
    write(get(FITS, io), arr, header = get(FITSHeader, hdr))
write(io::FitsIO, obj::FitsImage) =
    write(io, get(FitsHeader, obj), get(Array, obj))
write(io::FitsIO, args...) =
    for arg in args
        write(io, arg)
    end

# The following versions can build the header on the fly.
write(io::FitsIO, arg::Tuple{HeaderLike,AbstractArray{<:Real}}) =
        write(io, arg...)
write(io::FitsIO, arg::Tuple{AbstractArray{<:Real},HeaderLike}) =
    write(io, arg...)
write(io::FitsIO, arr::AbstractArray{<:Real}, hdr::HeaderLike) =
    write(io, hdr, arr)
function write(io::FitsIO,
               hdr::Tuple{Vararg{Pair{<:AbstractString}}},
               arr::AbstractArray)
    write(io::FitsIO, FitsHeader(hdr...), arr)
end
function write(io::FitsIO,
               hdr::NamedTuple,
               arr::AbstractArray)
    write(io::FitsIO, FitsHeader(; hdr...), arr)
end

"""
    EasyFITS.getfile(obj) -> file

yields the FITS file associated with object `obj` if it is open; otherwise
throws an error.  Object `obj` can be an instance of `FitsIO`, `FitsHDU`,
`FITSIO.FITS`, `FITSIO.HDU` or `FITSIO.FITSFile`.

    EasyFITS.getfile(obj, ext) -> file

throws an error if the FITS file associated with `obj` is not open; otherwise
moves to HDU number `ext` and returns the FITS file.  Object `obj` can be an
instance of `FitsIO`, `FITSIO.FITS`, or `FITSIO.FITSFile`.

This method is considered as *low-level* and is not exported.

"""
function getfile(file::FITSFile)
    fits_assert_open(file)
    return file
end

function getfile(file::FITSFile, ext::Integer)
    fits_assert_open(file)
    fits_movabs_hdu(file, ext)
    return file
end

get(::Type{FITSFile}, io::FitsIO) = getfield(get(FITS, io), :fitsfile)
get(::Type{FITSFile}, hdu::FitsHDU) = getfield(get(HDU, hdu), :fitsfile)
Base.isopen(obj::Union{FitsIO,<:FitsHDU}) = isopen(get(FITSFile, obj))
Base.close(obj::Union{FitsIO,<:FitsHDU}) = close(get(FITSFile, obj))

getfile(fh::FITS) = getfile(fh.fitsfile)
getfile(fh::FITS, ext::Integer) = getfile(fh.fitsfile, ext)

getfile(hdu::HDU) = getfile(hdu.fitsfile, hdu.ext)

getfile(io::FitsIO) = getfile(get(FITS, io))
getfile(io::FitsIO, ext::Integer) = getfile(get(FITS, io), ext)

getfile(hdu::FitsHDU) = getfile(get(HDU, hdu))

@noinline throw_file_already_exists(args::AbstractString...) =
    error(file_already_exists(args...))

file_already_exists(path::AbstractString) =
    string("file \"", path, "\" already exists")

file_already_exists(path::AbstractString, usage::AbstractString) =
    string(file_already_exists(path), ", ", usage)

end # module


text:
using COSMO
using SketchyCGAL
using LinearAlgebra, SparseArrays
using BSON
using Printf

include("../utils.jl")

# Run experiments single threaded
blas_threads = BLAS.set_num_threads(1)

## Data Load
function setup_maxcut_scgal(G)
    # Parameters
    n = size(G, 1)
    d = n

    # Data
    # objective = min -1/4( ∑_ij w_ij*(1-yi*yj) ) = -1/4⟨diag(∑ᵢ w_ij) - W), Y⟩
    # note that this reformulation uses the fact that Tr(Y) = 1
    C = -0.25*(Diagonal(G*ones(n)) - G)
    b = ones(n)

    # Scaling variables -- so trace is bounded by 1
    scale_C = 1 / norm(C)
    scale_X = 1 / n
    # b .= b .* scale_X
    # const C_const = C * scale_C

    # Linear map
    # AX = diag(X)
    function A!(y, X)
        n = size(X, 1)
        for i in 1:n
            y[i] = X[i,i]
        end
        return nothing
    end

    # Adjoint
    # A*z = Diagonal(z)
    function A_adj!(S::SparseMatrixCSC, z)
        for (i, j, v) ∈ zip(findnz(S)...)
            S[i,j] = 0.0
        end
        S[diagind(S)] .= z
        return nothing
     end

    # primative 1: u -> C*u
    function C_u!(v, u)
        mul!(v, C_const, u)
    end

    #primative 2: (u, z) -> (A*z)u
    function Aadj_zu!(v, u, z)
        v .= u .* z
    end

    #primative 3: u -> A(uu^T)
    function A_uu!(z, u)
        z .= u.*u
    end

    return C, b, A!, A_adj!, n, d, scale_X, scale_C
end

function maxcut_scgal(G, filename, trial_data; R=10)
    trial_data[filename] = Dict()
    C, b, A!, A_adj!, n, d, scale_X, scale_C = setup_maxcut_scgal(G)

    ηt(t) = 2.0/(t + 1.0)
    δt(t) = 1.0
    tt = @timed scgal_full(
        C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
        max_iters=2_000,
        print_iter=1_000,
        R=R,
        logging=true,
        # logging_primal=true,
        # ηt=ηt,
        # δt=δt
    )

    trial_data[filename]["dim"] = (
        n=n,
        nonzeros=nnz(G)
    )
    trial_data[filename]["stats"] = tt

    return nothing
end

function run_trials_size()
    filenames = readdir(joinpath(dirname(@__DIR__), "data/gset"))
    trial_data = Dict()

    for filename in filenames
        G = graph_from_file(joinpath(dirname(@__DIR__), "data/gset/$filename"))
        @printf("On %3s: %8s %12s\n", filename, size(G,1), nnz(G))
        maxcut_scgal(G, filename, trial_data)
    end
    return trial_data
end

trial_data = run_trials_size()
BSON.bson(joinpath(@__DIR__, "output/trials_size.bson"), trial_data)

ext: jl


text:
"""
    collision_checker(veh_a::Entity, veh_b::Entity)
    collision_checker(veh_a, veh_b, veh_a_def::AbstractAgentDefinition, veh_b_def::AbstractAgentDefinition)
return True if `veh_a` and `veh_b` collides.
Relies on the parallel axis theorem.
"""
function collision_checker(veh_a::Entity, veh_b::Entity)
    return collision_checker(veh_a.state, veh_b.state, veh_a.def, veh_b.def)
end

function collision_checker(veh_a, veh_b, veh_a_def::AbstractAgentDefinition, veh_b_def::AbstractAgentDefinition)
    center_a = posg(veh_a)
    center_b = posg(veh_b)
    l_a = length(veh_a_def)
    w_a = width(veh_a_def)
    l_b = length(veh_b_def)
    w_b = width(veh_b_def)
    # first fast check:
    @fastmath begin
        Δ = sqrt((posg(veh_a).x - posg(veh_b).x)^2 + (posg(veh_a).y - posg(veh_b).y)^2)
        r_a = sqrt(l_a*l_a/4 + w_a*w_a/4)
        r_b = sqrt(l_b*l_b/4 + w_b*w_b/4)
    end
    if Δ ≤ r_a + r_b
        # fast check is true, run parallel axis theorem
        Pa = polygon(center_a, veh_a_def)
        Pb = polygon(center_b, veh_b_def)
        return overlap(Pa, Pb)
    end
    return false
end

"""
    collision_checker(scene::Frame{Entity{S,D,I}}, egoid::I) where {S, D<:AbstractAgentDefinition, I}
return true if any entity in the scene collides with the entity of id `egoid`.
"""
function collision_checker(scene::Frame{Entity{S,D,I}}, egoid::I) where {S, D<:AbstractAgentDefinition, I}
    ego = get_by_id(scene, egoid)
    for veh in scene
        if veh.id != egoid
            if collision_checker(veh, ego)
                return true
            end
        end
    end
    return false
end

"""
    polygon(pos::VecSE2{Float64}, veh_def::AbstractAgentDefinition)
    polygon(x::Float64,y::Float64,theta::Float64, length::Float64, width::Float64)
returns a 4x2 static matrix corresponding to a rectangle around a car
centered at `pos` and of dimensions specified by `veh_def`
"""
function polygon(pos::VecSE2{Float64}, veh_def::AbstractAgentDefinition)
    x, y ,θ = pos
    return polygon(x, y, θ, length(veh_def), width(veh_def))
end

function polygon(x::Float64,y::Float64,theta::Float64, length::Float64, width::Float64)
    l_ = length/2
    w_ = width/2

    # compute each vertex:
    @fastmath begin
        #top left
        Ax = x + l_*cos(theta) - w_*sin(theta)
        Ay = y + l_*sin(theta) + w_*cos(theta)

        #top right
        Bx = x + l_*cos(theta) + w_*sin(theta)
        By = y + l_*sin(theta) - w_*cos(theta)

        # bottom right
        Cx = x - l_*cos(theta) + w_*sin(theta)
        Cy = y - l_*sin(theta) - w_*cos(theta)

        # bttom left
        Dx = x - l_*cos(theta) - w_*sin(theta)
        Dy =  y - l_*sin(theta) + w_*cos(theta)
    end
    P = @SMatrix [Ax Ay;
                  Bx By;
                  Cx Cy;
                  Dx Dy]

    return P
end

### POLYGON INTERSECTION USING PARALLEL AXIS THEOREM

"""
    overlap(poly_a::SMatrix{4, 2, Float64}, poly_b::SMatrix{4, 2, Float64})
Check if two convex polygons overlap, using the parallel axis theorem
a polygon is a nx2 matrix where n in the number of verteces
http://gamemath.com/2011/09/detecting-whether-two-convex-polygons-overlap/
  /!\\ edges needs to be ordered
"""
function overlap(poly_a::SMatrix{M, 2, Float64},
                poly_b::SMatrix{N, 2, Float64}) where {M,N}

    if find_separating_axis(poly_a, poly_b)
        return false
    end
    if find_separating_axis(poly_b, poly_a)
        return false
    end

    return true

end

"""
    find_separating_axis(poly_a::SMatrix{4, 2, Float64}, poly_b::SMatrix{4, 2, Float64})
build a list of candidate separating axes from the edges of a
  /!\\ edges needs to be ordered
"""
function find_separating_axis(poly_a::SMatrix{M, 2, Float64},
                              poly_b::SMatrix{N, 2, Float64}) where {M, N}
    n_a = size(poly_a)[1]
    n_b = size(poly_b)[1]
    axis = zeros(2)
    @views previous_vertex = poly_a[end,:]
    for i=1:n_a # loop through all the edges n edges
        @views current_vertex = poly_a[i,:]
        # get edge vector
        edge = current_vertex - previous_vertex

        # rotate 90 degrees to get the separating axis
        axis[1] = edge[2]
        axis[2] = -edge[1]

        #  project polygons onto the axis
        a_min,a_max = polygon_projection(poly_a, axis)
        b_min,b_max = polygon_projection(poly_b, axis)

        # check separation
        if a_max < b_min
            return true#,current_vertex,previous_vertex,edge,a_max,b_min

        end
        if b_max < a_min
            return true#,current_vertex,previous_vertex,edge,a_max,b_min

        end

        @views previous_vertex = poly_a[i,:]
    end

    # no separation was found
    return false
end

"""
    polygon_projection(poly::SMatrix{4, 2, Float64}, axis::Vector{Float64})
return the projection interval for the polygon poly over the axis axis
"""
function polygon_projection(poly::SMatrix{N, 2, Float64},
                            axis::Vector{Float64}) where {N}
        n_a = size(poly)[1]
        @inbounds @fastmath @views d1 = dot(poly[1,:],axis)
        @inbounds @fastmath @views d2 = dot(poly[2,:],axis)
        # initialize min and max
        if d1<d2
            out_min = d1
            out_max = d2
        else
            out_min = d2
            out_max = d1
        end

        for i=1:n_a
            @inbounds @fastmath @views d = dot(poly[i,:],axis)
            if d < out_min
                out_min = d
            elseif d > out_max
                out_max = d
            end
        end
        return out_min,out_max
end


text:
@marginalrule NormalMeanVariance(:out_μ_v) (m_out::NormalDistributionsFamily, m_μ::PointMass, m_v::PointMass) = begin
    return (out = prod(ProdAnalytical(), NormalMeanVariance(mean(m_μ), mean(m_v)), m_out), μ = m_μ, v = m_v)
end

@marginalrule NormalMeanVariance(:out_μ_v) (m_out::PointMass, m_μ::NormalDistributionsFamily, m_v::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, NormalMeanVariance(mean(m_out), mean(m_v))), v = m_v)
end

@marginalrule NormalMeanVariance(:out_μ_v) (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily, m_v::PointMass) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_μ, W_μ     = weightedmean_precision(m_μ)

    W_bar = inv(mean(m_v))

    W  = [ W_out + W_bar -W_bar; -W_bar W_μ + W_bar ]
    xi = [ xi_out; xi_μ ]

    return (out_μ = MvNormalWeightedMeanPrecision(xi, W), v = m_v)
end

@marginalrule NormalMeanVariance(:out_μ) (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily, q_v::Any) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_μ, W_μ     = weightedmean_precision(m_μ)

    W_bar = inv(mean(q_v))

    W  = [ W_out + W_bar -W_bar; -W_bar W_μ + W_bar ]
    xi = [ xi_out; xi_μ ]

    return MvNormalWeightedMeanPrecision(xi, W)
end

@marginalrule NormalMeanVariance(:out_μ) (m_out::PointMass, m_μ::UnivariateNormalDistributionsFamily, q_v::Any) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), NormalMeanVariance(mean(m_out), mean(q_v)), m_μ))
end

@marginalrule NormalMeanVariance(:out_μ) (m_out::UnivariateNormalDistributionsFamily, m_μ::PointMass, q_v::Any) = begin
    return (out = prod(ProdAnalytical(), NormalMeanVariance(mean(m_μ), mean(q_v)), m_out), μ = m_μ)
end


text:
#= ForwardDiff is removed as a dependency from this version
import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, PT, finish
import MechanicalSketch: settext, place_image, PALETTE, color_with_lumin, snapshot
import MechanicalSketch: circle, ∙, arrow_nofill, Point, @layer, sethue, draw_expr
import MechanicalSketch: latexify, arrow, generate_complex_potential_vortex
import MechanicalSketch: generate_complex_potential_source, ∇, color_matrix
import MechanicalSketch: fontsize, draw_legend, draw_streamlines, ∇_rectangle
import MechanicalSketch: @eq_draw, @ev_draw, x_y_iterators_at_pixels
import MechanicalSketch: matrix_to_function, SA
import ColorSchemes:     isoluminant_cgo_80_c38_n256

@import_expand(~m, s, N, °)
include("test_functions_33.jl")
ϕ = ϕ_33
include("test_functions_39.jl")

# Scaling and placement
physwidth = 10.0m
physheight = 4.0m
framedheight = physheight * 1.2
totheight = framedheight * 4
totwidth = totheight * WI / HE
Δx = totwidth / 4
Δy = framedheight / 2
Oadj = O + (-3EM, EM /3)
cent1 = Oadj + ( Δx, 3Δy)
cent2 = Oadj + ( Δx, Δy)
cent3 = Oadj + ( Δx, -Δy)
cent4 = Oadj + ( Δx, -3Δy)

empty_figure( backgroundcolor = color_with_lumin(PALETTE[5], 90),
      hue = "black", height = totheight);
@eq_draw "Potential field, defined previously:" :(ϕ(𝐫))  init= true
@eq_draw "Steady flow velocity field at r:" :( 𝐕(𝐫) = 𝐢∙u(𝐫) + 𝐣∙v(𝐫) = ∇ϕ(𝐫))
Δx, Δy = physwidth / 2, physheight / 2
rng = string([(-Δx, -Δy), (Δx, Δy)])
@ev_draw "We evaluate r ∈ rectangle $rng" :(𝐕 = ∇(ϕ))

@eq_draw("Local acceleration, steady flow is", quote
    𝐚 = d𝐕 / dt = d( 𝐢∙u(𝐫) + 𝐣∙v(𝐫)) / dt
end)
@eq_draw("A particle moves, so r = r(x, y, t)", quote
    𝐚 = δ𝐕 / δt ∙ dt / dt + δ𝐕 / δx ∙ dx / dt + δ𝐕 / δy ∙ dy / dt
end)
@eq_draw("This field is steady over time, so the first term is zero, and dx / dt = u", quote
    𝐚 = δ𝐕 / δx ∙ u + δ𝐕 / δy ∙ v = ∇𝐕∙(u, v)
end)

@eq_draw("The last step is valid because the axes are orthogonal", quote
    𝐚 = δu / δx ∙ u +  δv / δx ∙ v = ∇(u)∙ u + ∇(v)∙ v
end)

#= To speed up the calculations, we'll strip and reapply interpolation information
xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
u_c = matrix_to_function([𝐕(x, y)[1] for y in ys, x in xs])
u(r) = u_c(real(r), imag(r))
v_c = matrix_to_function([𝐕(x, y)[2] for y in ys, x in xs])
v(r) = u_c(real(r), imag(r))
∇_rectangle(v; physwidth, physheight)
∇_rectangle(v; physwidth, physheight)
∇(ϕ; physwidth, physheight)
=#
cpt = O + ( EM, -HE / 2 + 2EM)
@eq_draw("Defining functions", quote
    u(r) = 𝐕(real(r), imag(r))[1]
    v(r) = 𝐕(real(r), imag(r))[2]
    ∇u(r) = ∇(u)(real(r), imag(r))
    ∇v(r) = ∇(v)(real(r), imag(r))
#    𝐚(r) = ∇u(r).∙u(r) + ∇v(r).∙v(r)
#    𝐚(x, y) = 𝐚(complex(x, y))
end)

r₀ = complex(0.0m, 0.0m)
#a₀ , u₀, v₀ = 𝐚(r₀), u(r₀), v(r₀)
#@show(a₀ , u₀, v₀);
#u₀, v₀ = u(r₀), v(r₀)
#@show(u₀, v₀);

xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
import MechanicalSketch.MechanicalUnits.Length
xss = SA[collect(xs)[1:2]]
yss = SA[collect(ys)[1:2]]

#xss = SA{typeof(0.0)}[collect(1.0:1.1:10.0)[1:2]]

um = SA{typeof(u₀)}[𝐕(x, y)[1] for y in ys, x in xs]
vm = [𝐕(x, y)[1] for y in ys, x in xs]
xss = SA[xs[1:2]]
yss = SA[ys[1:2]]
foom() = [∇(u)(x, y)[1] for y in yss, x in xss]
∇um = foom()
#ERROR: MethodError: no method matching
#zero(::Type{Tuple{Quantity{Float64,  ᵀ⁻¹, FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}, Quantity{Float64,  ᵀ⁻¹, FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}}})
# The above can work, but is incredibly slow. Maybe we can improve by defining the missing method?
import MechanicalSketch.MechanicalUnits.Unitfu
import MechanicalSketch.MechanicalUnits.FreeUnits
import MechanicalSketch.MechanicalUnits.dimension
import MechanicalSketch.MechanicalUnits.Time
import MechanicalSketch.MechanicalUnits.Dimensions
import MechanicalSketch.MechanicalUnits.Dimension
import MechanicalSketch.MechanicalUnits.Level
import MechanicalSketch.Quantity
TD = dimension(1.0s)
ITD = dimension(1.0/s)
typeof(InvTime)

typeof(Time)
InvTime isa Time^-1

#It seems we havent defined zero(T), where T =
# Tuple{Unitfu.Quantity{Float64,  ᵀ⁻¹, Unitfu.FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}, Unitfu.Quantity{Float64,  ᵀ⁻¹, Unitfu.FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}}
# Tuple{Quantity{Float64,  ᵀ⁻¹, FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}, Quantity{Float64,  ᵀ⁻¹, FreeUnits{(s⁻¹,),  ᵀ⁻¹, nothing}}}
# Tuple{Quantity{Float64,  ITD, FreeUnits{(s⁻¹,),  ITD, nothing}}, Quantity{Float64,  ITD, FreeUnits{(s⁻¹,),  ITD, nothing}}}
# Tuple{Quantity{Float64,  ITD, U}, Quantity{Float64,  ITD, U}} where U
typeof(1.0s) <: Quantity{T,  TD, U} where {T, U}
typeof(1.0s) <: Quantity{Float64,  TD, U} where {U}
typeof((1.0s, 1.0s)) <: Tuple{Quantity{Float64,  TD, U}, Quantity{Float64,  TD, U}} where {U}
typeof((1.0/s, 1.0/s)) <: Tuple{Quantity{Float64,  ITD, U}, Quantity{Float64,  ITD, U}} where {U}

# But zero is not defined for Tuples anywhere.... We have 44 methods.

zero((1.0s, 1.0s))
zero(typeof((1.0, 1.0)))
∇vm

zero(1.0m)
zero(typeof(1.0m))

fooam() = [𝐚(x, y) for y in ys, x in xs]
am = @time fooam()
# Show velocity field

colmat1, legend1 = color_matrix(𝐕; centered = true, maxlegend = 1.0m/s, name = Symbol("‖V‖"));
ulp, lrp = place_image(cpt, colmat1; centered = false)
@layer begin
    fontsize(25)
    draw_legend(ulp + (0.5EM, 0.5EM) + (physwidth, 0.0m), legend1)
    # To better indicate positive direction for this static image:
    sethue(PALETTE[8])
    draw_streamlines(cent2, 𝐕, probability = 0.001)
end
snapshot()

# Show acceleration field
cpt = O + (ulp[1], lrp[2] - EM)

@time colmat2, legend2 = color_matrix(𝐚; centered = true, noofbins = 7, name = Symbol("‖a‖"));
println("The calculation is done?")
ulp, lrp = place_image(cpt, colmat2; centered = false)
@layer begin
    fontsize(25)
    draw_legend(ulp + (0.5EM, 0.5EM) + (physwidth, 0.0m), legend2)
    # To better indicate positive direction for this static image:
    sethue(PALETTE[8])
    draw_streamlines(cent2, 𝐚, probability = 0.001)
end

snapshot(fname = joinpath(@__DIR__, "test_38.png"))
=#


text:
"""
    FactorisationWrapper

A small wrapper around a a factorisation `fact`, allowing to represent
multiplication by the inverse of `fact`.
"""
struct FactorisationWrapper{Fact}
    fact::Fact
end

Base.size(fact::FactorisationWrapper) = size(fact.fact)
Base.eltype(fact::FactorisationWrapper) = eltype(fact.fact)
mul!(dest, fact::FactorisationWrapper, u) = ldiv!(dest, fact.fact, u)

#TODO
#struct AdaptiveFilter{T<:Real} <: AbstractFilter{T}
#
#end

"""
    ConstantFilter

Represents the action of a modal filter on values in a nodal basis with fixed
strength.
"""
struct ConstantFilter{T<:Real,Nodal2Modal,Modal2Nodal,Tmp,FilterFunction} <: AbstractFilter{T}
    coefficients::Vector{T}
    nodal2modal::Nodal2Modal
    modal2nodal::Modal2Nodal
    tmp::Tmp
    filter::FilterFunction

    function ConstantFilter(coefficients::Vector{T}, nodal2modal::Nodal2Modal, modal2nodal::Modal2Nodal, tmp::Tmp, filter::FilterFunction) where {T<:Real,Nodal2Modal,Modal2Nodal,Tmp,FilterFunction}
        @argcheck length(coefficients) == size(tmp,1)

        new{T,Nodal2Modal,Modal2Nodal,Tmp,FilterFunction}(coefficients, nodal2modal, modal2nodal, tmp, filter)
    end
end

(filter::ConstantFilter)(u) = filter(u, filter.tmp)

function (filter::ConstantFilter)(u::AbstractVector, tmp::AbstractVector)
    @unpack coefficients, nodal2modal, modal2nodal = filter
    @boundscheck begin
        length(coefficients) == length(tmp)
    end

    mul!(tmp, nodal2modal, u)
    @inbounds tmp .*= coefficients
    mul!(u, modal2nodal, tmp)

    nothing
end

function (filter::ConstantFilter)(u::AbstractMatrix, tmp::AbstractVector)
    @boundscheck begin
        length(tmp) == size(u,1)
    end

    for j in Base.OneTo(size(u,2))
        v = view(u, :, j)
        @inbounds filter(v, tmp)
    end

    nothing
end

function (filter::ConstantFilter)(u::AbstractMatrix, tmp::AbstractMatrix)
    @unpack coefficients, nodal2modal, modal2nodal = filter
    @boundscheck begin
        size(tmp) == size(u)
    end

    mul!(tmp, nodal2modal, u)
    @inbounds tmp .*= coefficients
    mul!(u, modal2nodal, tmp)

    nothing
end

"""
    ExponentialFilter

Represents the exponential filter function `σ(η) = exp(-α*η^p)`.
"""
struct ExponentialFilter{T<:Real} <: AbstractFilterFunction
    α::T
    p::Int

    function ExponentialFilter(α::T, p::Int=2) where {T<:Real}
        α < 0 && warn("α should be nonnegative [α = $α].")
        p < 0 && warn("p should be nonnegative [p = $p].")
        new{T}(α, p)
    end
end

function ExponentialFilter(::Type{T}=Float64, p::Int=2) where {T<:Real}
    α = -log(eps(T))
    ExponentialFilter(α, p)
end

function ExponentialFilter(p::Int, ::Type{T}=Float64) where {T<:Real}
    α = -log(eps(T))
    ExponentialFilter(α, p)
end

function set_filter_coefficients!(coefficients::Vector, filter::ExponentialFilter)
    T = eltype(coefficients)
    N = length(coefficients) - 1
    @unpack α, p = filter

    @inbounds for i in Base.OneTo(length(coefficients))
        η = T(i-1) / N
        coefficients[i] = exp(-α*η^p)
    end

    nothing
end


text:
## ----------------------------------------------------------------------------
# EP biomass corr
let
    ps = Plots.Plot[]
    for method in ALL_METHODS
        p = plot(title = string(iJR.PROJ_IDER, " method: ", method),
            xlabel = "exp biom", ylabel = "model biom")
        ep_vals = DAT[method, :ep, :flx, iJR.BIOMASS_IDER, EXPS]
        eperr_vals = DAT[method, :eperr, :flx, iJR.BIOMASS_IDER, EXPS]
        Hd_vals = DAT[method, :Hd, :flx, iJR.BIOMASS_IDER, EXPS]
        m, M = myminmax([Hd_vals; ep_vals])
        margin = abs(M - m) * 0.1
        scatter!(p, Hd_vals, ep_vals;
            yerr = eperr_vals,
            label = "", color = :black,
            alpha = 0.7, ms = 7,
            xlim = [m - margin, M + margin],
            ylim = [m - margin, M + margin],
        )
        push!(ps, p)
    end
    layout = (1, length(ps))
    mysavefig(ps, "obj_val_ep_corr"; layout)
end

## -------------------------------------------------------------------
# correlations
let

    tot_ps = Plots.Plot[]
    for method in ALL_METHODS
        # total corr
        let
            ep_vals = DAT[method, :ep, :flx, FLX_IDERS, EXPS] .|> abs
            ep_errs = DAT[method, :eperr, :flx, FLX_IDERS, EXPS] .|> abs
            Hd_vals = DAT[method, :Hd, :flx, FLX_IDERS, EXPS] .|> abs

            color = [ider_colors[ider] for ider in FLX_IDERS, exp in EXPS]
            scatter_params = (;label = "", color, ms = 7, alpha = 0.7)
            # ep corr
            p = plot(title = "$(iJR.PROJ_IDER) (EP) $method",
                ylabel = "model abs flx",
                xlabel = "exp abs flx",
            )
            scatter!(p, Hd_vals, ep_vals; yerr = ep_errs, scatter_params...)
            all_vals = [ep_vals; Hd_vals] |> sort!
            plot!(p, all_vals, all_vals; ls = :dash, color = :black, label = "")
            push!(tot_ps, deepcopy(p))
        end

        # per ider
        let
            for ider in FLX_IDERS
                ep_vals = DAT[method, :ep, :flx, ider, EXPS] .|> abs
                ep_errs = DAT[method, :eperr, :flx, ider, EXPS] .|> abs
                Hd_vals = DAT[method, :Hd, :flx, ider, EXPS] .|> abs

                color = ider_colors[ider]
                # ep corr
                p = plot(title = "$(iJR.PROJ_IDER) (EP) $method",
                    ylabel = "model abs flx",
                    xlabel = "exp abs flx",
                )
                scatter!(p, Hd_vals, ep_vals; yerr = ep_errs,
                    label = "", color, ms = 7, alpha = 0.7
                )
                bounds = DAT[method, :bounds, :flx, ider, EXPS]
                all_vals = [ep_vals; Hd_vals] |> sort!
                plot!(p, all_vals, all_vals;
                    ls = :dash, color = :black, label = "",
                    # xlim = [lb, ub], ylim = [lb, ub]
                )
                mysavefig(p, "corr"; ider, method)
            end
        end
    end

    layout = (1, length(tot_ps))
    mysavefig(tot_ps, "flx_tot_corr"; layout)

end


text:
"""
    EnhancedConsoleLogger(stream::IO==stderr; kwargs...)

Replacement for the standard library `ConsoleLogger` that adds several usability
improvements.

# Additional keyword arguments
The `EnhancedConsoleLogger` supports all of the standard keyword arguments accepted
by `ConsoleLogger`. In addition, special behavior is defined for the following arguments:
 *  `_pid`: prints in a column on the top right side of a log message. This keyword is
    used internally by `WorkerLogger` to display the originating process for a
    log message.
 *  `_overwrite`: If `true`, then repetitions of the log message will be printed over to
    avoid filling the screen with log messages.
 *  `progress`: For a progress between 0 and 1, draw a progress bar on the right side
    of the log message.
 *  `_showlocation`: If true, print the location of the log message. Defaults to true for
    `Debug`, `Warning`, and `Error` logs, and to false for `Info` logs.
"""
mutable struct EnhancedConsoleLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel

    width::Int
    show_limited::Bool
    message_limits::Dict{Any, Int}

    last_id::Symbol
    last_length::Int
end
function EnhancedConsoleLogger(stream::IO=stderr; show_limited=true, width=80,
                               min_level::LogLevel=ProgressLevel)
    EnhancedConsoleLogger(stream, min_level, width, show_limited, Dict{Any, Int}(),
                          :nothing, 0)
end

"""
    progress_string(progress, width)

Returns a string with a completion percentage and a progress bar. Uses Unicode characters
to render sub-character progress increments.
"""
function progress_string(progress, width)
    progress_chars = [' ', '▏', '▎', '▍', '▌', '▋',  '▊', '▉', '█']

    progress_bar_width = width - 9
    progress_width = progress_bar_width*progress

    if progress_width >= progress_bar_width
        return @sprintf "%5.1f%% ║%s" 100 '█'^progress_bar_width
    end
    full_blocks = floor(Int, progress_width)
    remaining_width = progress_width - full_blocks
    trailing_char = progress_chars[floor(Int, remaining_width*8)+1]

    blank_blocks = progress_bar_width-full_blocks-1
    blank_blocks < 0 && (blank_blocks = 0)

    @sprintf "%5.1f%% ║%s" progress*100 ('█'^full_blocks)*trailing_char*(' '^blank_blocks)
end

function log_label(level)
    level == Logging.Warn && return "Warning"
    level == ProgressLevel && return "Progress"
    return string(level)
end

function log_color(level)
    color = Logging.default_logcolor(level)
    level == ProgressLevel && (color = :green)
    return color
end

function log_location(mod, group, id, file, line)
    location = ""
    mod !== nothing && (location *= "$mod")
    if file !== nothing
        mod !== nothing && (location *= " ")
        location *= Base.contractuser(file)
        if line !== nothing
            location *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    return location
end

function log_message(logger::EnhancedConsoleLogger, msg, kwargs)
    msglines = [(indent=0,msg=l) for l in split(chomp(string(msg)), '\n')]
    dsize = displaysize(logger.stream)
    if !isempty(kwargs)
        valbuf = IOBuffer()
        rows_per_value = max(1, dsize[1]/(length(kwargs)+1))
        valio = IOContext(IOContext(valbuf, logger.stream),
                          :displaysize => (rows_per_value, dsize[2]-5),
                          :limit => logger.show_limited)
        for (key,val) in pairs(kwargs)
            key == :progress      && continue
            key == :_pid          && continue
            key == :_overwrite    && continue
            key == :_showlocation && continue
            Logging.showvalue(valio, val)
            vallines = split(String(take!(valbuf)), '\n')
            if length(vallines) == 1
                push!(msglines, (indent=2,msg=SubString("$key = $(vallines[1])")))
            else
                push!(msglines, (indent=2,msg=SubString("$key =")))
                append!(msglines, ((indent=3,msg=l) for l in vallines))
            end
        end
    end

    return msglines
end

function handle_message(logger::EnhancedConsoleLogger, level, message, mod, group,
                        id, file, line; maxlog=nothing, kwargs...)
    if maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    label = log_label(level)
    color = log_color(level)
    msglines = log_message(logger, message, kwargs)

    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)

    should_overwrite = get(kwargs, :_overwrite, level == ProgressLevel)
    if should_overwrite && id == logger.last_id
        print(iob, "\u1b[A"^logger.last_length)
    end

    pid_string = haskey(kwargs, :_pid) ? "║"*string(kwargs[:_pid]) : ""
    justify_width = min(logger.width, displaysize(logger.stream)[2])

    i = 1
    print_progress = (level == ProgressLevel)
    print_location = get(kwargs, :_showlocation,
                         level < ProgressLevel || level >= Logging.Warn)
    while i <= length(msglines)
        indent = msglines[i][1]
        msg    = msglines[i][2]
        linewidth = 2
        printstyled(iob, "║ ", bold=true, color=color)

        if i == 1 && !isempty(label)
            printstyled(iob, label, " ", bold=true, color=color)
            linewidth += length(label) + 1
        end
        print(iob, ' '^indent, msg)
        linewidth += indent + length(msg)

        if print_progress
            prog_string = haskey(kwargs, :progress) ? progress_string(kwargs[:progress], 30) : "NO PROGRESS PROVIDED "
            prog_color  = haskey(kwargs, :progress) ? :green : :red
            npad = justify_width - linewidth - length(prog_string) -
                   (i == 1 ? length(pid_string) : 0) - 3

            if npad < 0
                insert!(msglines, i+1, (indent=0,msg=SubString("")))
            else
                printstyled(iob, ' ', '·'^npad, ' ', color=:light_black)
                printstyled(iob, prog_string, color=prog_color, bold=true)
                linewidth += npad + 32

                print_progress = false

                if i == length(msglines) && print_location
                    push!(msglines, (indent=0,msg=SubString("")))
                end
            end
        end

        if i == length(msglines) && print_location
            location = log_location(mod, group, id, file, line)

            npad = justify_width - linewidth - length(location) -
                   (i == 1 ? length(pid_string) : 0) - 4

            if npad < 0
                push!(msglines, (indent=0,msg=SubString("")))
            else
                printstyled(iob, ' ', '·'^npad, ' ', location, color=:light_black)
                linewidth += npad + length(location) + 2

                print_location = false
            end
        end

        npad = max(0, justify_width - linewidth - (i == 1 ? length(pid_string) : 0) - 1)
        printstyled(iob, ' '^npad, bold=true, color=color)
        if i == 1 && pid_string != ""
            printstyled(iob, pid_string, bold=true, color=color)
            linewidth += length(pid_string)
        end
        printstyled(iob, "║", bold=true, color=color)
        println(iob)

        i += 1
    end

    logger.last_id = id
    logger.last_length = length(msglines)
    write(logger.stream, take!(buf))
    nothing
end

shouldlog(::EnhancedConsoleLogger, args...) = true
min_enabled_level(logger::EnhancedConsoleLogger) = logger.min_level
catch_exceptions(::EnhancedConsoleLogger) = false


text:
using UnconstrainedMPCCActif
using ActifMPCCmod

UnconstrainedMPCCActif_success=true
println("On teste le module UnconstrainedMPCCActif")

r=1.0;s=1.0;t=0.0; #relaxation KDB

nlp=ADNLPModel(f,ones(2), lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1)
nlp_ma=ADNLPModel(x->mpcc.mp.f(x)+norm(G(x)-x[3])^2+norm(H(x)-x[4])^2, [-1.0;1.0;-r;1.0],lvar=[mpcc.mp.meta.lvar;-r*ones(2*mpcc.nb_comp)], uvar=[mpcc.mp.meta.uvar;Inf*ones(2*mpcc.nb_comp)])
ma=ActifMPCCmod.MPCC_actif(nlp_ma,r,s,t,1)

#Bilan : donner une sortie bilan
if UnconstrainedMPCCActif_success==true
 println("UnconstrainedMPCCActif.jl passes the test !")
else
 println("UnconstrainedMPCCActif.jl contains some error")
end


text:
_scip = JuMP.optimizer_with_attributes(SCIP.Optimizer, "display/verblevel" => 0)
_glpk = GLPK.Optimizer

@testset "Solve semi decomposition with nogoods on y (callback)" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    cbtoposolver = GndStr.CallbackTopo(_scip, _scip)

    @testset "low flow: very easy instance" begin
        inst = Instance(nodes, 1 * demand, bounds, diams)
        result = optimize(inst, topo, cbtoposolver)
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, smallest diameter enough:
        @test zsol[11, 1] ≈ 1.0
        @test zsol[13, 1] ≈ 1.0
        @test zsol[4, 1] ≈ 1.0
        @test sum(zsol) ≈ 3.0 # all others 0

        # uniq flow solution
        @test qsol[11] ≈ 50.0
        @test qsol[13] ≈ 20.0
        @test qsol[4] ≈ 30.0
        @test sum(qsol) ≈ 100.0 # all others 0
    end

    @testset "medium flow: difficult instance" begin
        inst = Instance(nodes, 10 * demand, bounds, diams)
        result = optimize(inst, topo, cbtoposolver)
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, need one large diameter
        @test sum(zsol[2, :]) == 1
        @test sum(zsol[7, :]) == 1
        @test sum(zsol[11, :]) == 1
        @test sum(zsol[13, :]) == 1
        @test sum(zsol[:, 1]) == 1
        @test sum(zsol[:, 2]) == 3
        @test sum(zsol) ≈ 4.0 # all others 0
    end

    @testset "high flow on triangle: infeasible instance" begin
        inst = Instance([Node(0,0), Node(50,0)],
                        [-1000, 1000],
                        [Bounds(60,80), Bounds(60,80)],
                        [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                        [Arc(1,3), Arc(1,2), Arc(2,3),
                         Arc(3,1), Arc(2,1), Arc(3,2)])

        result = optimize(inst, topo, cbtoposolver)
        @test result.status == MOI.INFEASIBLE
    end
end

@testset "Solve GBD (callback)" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    cbgbdsolver = GndStr.CallbackGBD(_scip, _glpk)

    @testset "low flow: very easy instance" begin
        inst = Instance(nodes, 1 * demand, bounds, diams)
        result = optimize(inst, topo, cbgbdsolver)
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, smallest diameter enough:
        @test zsol[11, 1] ≈ 1.0
        @test zsol[13, 1] ≈ 1.0
        @test zsol[4, 1] ≈ 1.0
        @test sum(zsol) ≈ 3.0 # all others 0

        # uniq flow solution
        @test qsol[11] ≈ 50.0
        @test qsol[13] ≈ 20.0
        @test qsol[4] ≈ 30.0
        @test sum(qsol) ≈ 100.0 # all others 0
    end

    @testset "medium flow: difficult instance" begin
        inst = Instance(nodes, 10 * demand, bounds, diams)
        result = optimize(inst, topo, cbgbdsolver)
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, need one large diameter
        @test sum(zsol[2, :]) == 1
        @test sum(zsol[7, :]) == 1
        @test sum(zsol[11, :]) == 1
        @test sum(zsol[13, :]) == 1
        @test sum(zsol[:, 1]) == 1
        @test sum(zsol[:, 2]) == 3
        @test sum(zsol) ≈ 4.0 # all others 0
    end

    @testset "high flow on triangle: infeasible instance" begin
        inst = Instance([Node(0,0), Node(50,0)],
                        [-1000, 1000],
                        [Bounds(60,80), Bounds(60,80)],
                        [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                        [Arc(1,3), Arc(1,2), Arc(2,3),
                         Arc(3,1), Arc(2,1), Arc(3,2)])

        result = optimize(inst, topo, cbgbdsolver)
        @test result.status == MOI.INFEASIBLE
    end
end