text:
using .JuliennedArrays

# import TensorCast: julienne_slice, julienne_glue, newcode, jcollect # because Revise doesn't track these

if isdefined(JuliennedArrays, :julienne) # old notation

    @inline julienne_slice(A::AbstractArray, code::Tuple) where {T,N} =
        collect(julienne(A, code))

    @inline julienne_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N} =
        align(A, code)

elseif isdefined(JuliennedArrays, :Slices) # new notation

    @inline julienne_slice(A::AbstractArray{T,N}, code::Tuple) where {T,N} =
        Slices(A, newcode(code)...)

    @inline julienne_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N} =
        jcollect(Align(A, newcode(code)...))

    @generated newcode(code::Tuple) = Expr(:tuple, [ c == Colon ? True() : False() for c in code.parameters ]...)

    # collect(Align(... and copy(Align(... both give errors
    function jcollect(A::JuliennedArrays.Align{T,N}) where {T,N}
        out = Array{T}(undef, size(A))
        @inbounds for I in CartesianIndices(A)
            out[I] = A[I.I...]
        end
        out
    end

else
    @warn "TensorCast failed to load functions from JuliennedArrays"
end


text:
enumtypes(::Type{ALIGN_PACKED}) = (UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64, UInt128, Int128)

alignof(::Type{ALIGN_PACKED}, ::Type{T}) where {T<:Union{Int8, Int16, Int32, Int64, Int128}} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{T}) where {T<:Union{UInt8, UInt16, UInt32, UInt64, UInt128}} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{T}) where {T<:Union{Cfloat, Cdouble, Clongdouble}} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{T}) where {T<:Union{Ptr, Cstring}} = 1

alignof(::Type{ALIGN_PACKED}, ::Type{CC}) where {CC<:Cconst} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{CE}) where {CE<:Cenum} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{CA}) where {CA<:Carray} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{CA}) where {CA<:Caggregate} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}} = 1
alignof(::Type{ALIGN_PACKED}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Cenum, <:Calignment, <:Tuple}} = 1

enumtypes(::Type{ALIGN_NATIVE}) = (UInt32, Int32, UInt64, Int64, UInt128, Int128)

alignof(::Type{ALIGN_NATIVE}, ::Type{CC}) where {CC<:Cconst} = alignof(ALIGN_NATIVE, nonconst(CC))
alignof(::Type{ALIGN_NATIVE}, ::Type{CE}) where {CE<:Cenum} = alignof(ALIGN_NATIVE, eltype(CE))
alignof(::Type{ALIGN_NATIVE}, ::Type{CA}) where {CA<:Carray} = alignof(ALIGN_NATIVE, eltype(CA))
alignof(::Type{ALIGN_NATIVE}, ::Type{CA}) where {CA<:Caggregate} = Ctypelayout(CA).align
alignof(::Type{ALIGN_NATIVE}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}} = Ctypelayout(spec).align
alignof(::Type{ALIGN_NATIVE}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Cenum, <:Calignment, <:Tuple}} = alignof(ALIGN_NATIVE, Cenumlayout(spec).type)

const (_i8a, _i16a, _i32a, _i64a, _f32a, _f64a) = let
	(i8a, i16a, i32a, i64a, f32a, f64a) = refs = ((Ref{UInt}() for i in 1:6)...,)
	ccall("jl_native_alignment",
		Nothing,
		(Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}),
		i8a, i16a, i32a, i64a, f32a, f64a
	)
	(Int(r[]) for r in refs)
end
alignof(::Type{ALIGN_NATIVE}, ::Type{UInt8})   = _i8a
alignof(::Type{ALIGN_NATIVE}, ::Type{UInt16})  = _i16a
alignof(::Type{ALIGN_NATIVE}, ::Type{UInt32})  = _i32a
alignof(::Type{ALIGN_NATIVE}, ::Type{UInt64})  = _i64a
alignof(::Type{ALIGN_NATIVE}, ::Type{Float32}) = _f32a
alignof(::Type{ALIGN_NATIVE}, ::Type{Float64}) = _f64a
alignof(::Type{ALIGN_NATIVE}, ::Type{<:Ptr})   = alignof(ALIGN_NATIVE, sizeof(Ptr{Cvoid}) == sizeof(UInt32) ? UInt32 : UInt64)
alignof(::Type{ALIGN_NATIVE}, ::Type{Cstring}) = alignof(ALIGN_NATIVE, Ptr)
alignof(::Type{ALIGN_NATIVE}, ::Type{S}) where {S<:Signed} = alignof(ALIGN_NATIVE, unsigned(S))
alignof(::Type{ALIGN_NATIVE}, ::Type{UInt128}) = 2*alignof(ALIGN_NATIVE, UInt64)
alignof(::Type{ALIGN_NATIVE}, ::Type{Clongdouble}) = 2*alignof(ALIGN_NATIVE, Cdouble)

padding(::Type{ALIGN_PACKED}, offset::Int, align::Int) = (align%8) == 0 ? padding(ALIGN_NATIVE, offset, align) : 0
padding(::Type{ALIGN_PACKED}, offset::Int, bits::Int, typ) = bits == 0 ? padding(ALIGN_PACKED, offset, checked_alignof(ALIGN_PACKED, typ)*8) : 0

padding(::Type{ALIGN_NATIVE}, offset::Int, align::Int) = -offset & (align - 1)
function padding(::Type{ALIGN_NATIVE}, offset::Int, bits::Int, typ)
	pad = padding(ALIGN_NATIVE, offset, checked_alignof(ALIGN_NATIVE, typ)*8)
	return 0 < bits <= pad ? 0 : pad
end

function checked_alignof(args...)
	a = alignof(args...)
	a == 0 || a == nextpow(2, a) || error("Alignment must be a power of 2")
	return a
end

mutable struct Cenumlayout
	type::DataType
	min::Integer
	max::Integer
	values::Dict{Symbol, Integer}

	Cenumlayout() = new(Nothing, 0, 0, Dict{Symbol, Integer}())
end
Cenumlayout(::Type{CE}) where {CE<:Cenum} = Cenumlayout(Ctypespec(CE))

function _addvalue(layout::Cenumlayout, ::Type{Pair{sym, val}}, ::Type{spec}) where {spec<:Ctypespec, sym, val}
	haskey(layout.values, sym) && error("Encountered a duplicate value name `$(sym)` in enum specification")
	layout.values[sym] = val

	(min, max) = extrema(values(layout.values))
	for typ in enumtypes(strategy(spec))
		if typemin(typ) <= min && max <= typemax(typ)
			layout.type = typ
			layout.min = typ(min)
			layout.max = typ(max)
			return
		end
	end
	error("Unable to determine suitable enumeration storage type for `$(sym) = $(val)`")
end

@generated function Cenumlayout(::Type{spec}) where {spec<:Ctypespec{<:Any, <:Cenum, <:Calignment, <:Tuple}}
	layout = Cenumlayout()

	values = specification(spec)
	while values !== Tuple{}
		value = Base.tuple_type_head(values)
		values = Base.tuple_type_tail(values)

		_addvalue(layout, value, spec)
	end

	return layout
end

struct Ctypefield
	ind::Int  # field index
	type::DataType  # field type
	size::Int  # in bits:  0 means use sizeof(type)*8, >0 means bit field size
	offset::Int  # in bits
end

mutable struct Ctypelayout
	align::Int  # in bytes
	size::Int  # in bits
	offset::Int  # in bits
	fields::Dict{Symbol, Ctypefield}  # symbol => Ctypefield

	Ctypelayout() = new(1, 0, 0, Dict{Symbol, Ctypefield}())
end
Ctypelayout(::Type{CC}) where {CC<:Cconst} = Ctypelayout(nonconst(CC))
Ctypelayout(::Type{CA}) where {CA<:Caggregate} = Ctypelayout(Ctypespec(CA))

_offset(::Type{spec}, offset::Int) where {spec<:Ctypespec} = _offset(kind(spec), offset)
_offset(::Type{Cstruct}, offset::Int) = offset
_offset(::Type{Cunion}, offset::Int) = 0

_size(::Type{spec}, args...) where {spec<:Ctypespec} = _size(kind(spec), args...)
_size(::Type{Cstruct}, args...) = +(args...)
_size(::Type{Cunion}, args...) = max(args...)

_field(::Type{T}) where {T} = error("Unable to handle field type `$(T)` found in type specification")
_field(::Type{Pair{sym, Tuple{T}}}) where {sym, T} = (sym, T, 0)
_field(::Type{Pair{sym, Tuple{T, B}}}) where {sym, T, B} = (sym, T, B)
_field(::Type{Pair{sym, TS}}) where {sym, TS<:Ctypespec} = (sym, TS, 0)
_field(::Type{TS}) where {TS<:Ctypespec} = (nothing, TS, 0)

function _addfield(layout::Ctypelayout, sym::Symbol, field::Ctypefield)
	haskey(layout.fields, sym) && error("Encountered a duplicate field name `$(sym)` in type specification")
	layout.fields[sym] = field
end

_addfield(layout::Ctypelayout, ::Nothing, typ) = error("Encountered an unnamed field of type `$(typ)` in type specification")

function _addfield(layout::Ctypelayout, ::Nothing, ::Type{spec}, bits) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}}
	for (sym, field) in Ctypelayout(spec).fields
		_addfield(layout, sym, Ctypefield(length(layout.fields), (type(spec) <: Cconst ? Cconst(field.type) : field.type), field.size, layout.offset+field.offset))
	end
	return sizeof(type(spec))*8
end

function _addfield(layout::Ctypelayout, sym::Symbol, ::Type{spec}, bits) where {spec<:Ctypespec{<:Any, <:Copaques, <:Calignment, <:Tuple}}
	_addfield(layout, sym, Ctypefield(length(layout.fields), type(spec), 0, layout.offset))
	return sizeof(type(spec))*8
end

function _addfield(layout::Ctypelayout, sym::Symbol, typ, bits)
	_addfield(layout, sym, Ctypefield(length(layout.fields), typ, bits, layout.offset))
	return iszero(bits) ? sizeof(typ)*8 : bits
end

function _addfield(layout::Ctypelayout, ::Type{T}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}, T}
	(sym, typ, bits) = _field(T)
	pad = padding(strategy(spec), layout.offset, bits, typ)
	align = checked_alignof(strategy(spec), typ)
	layout.offset = _size(spec, layout.offset, pad)

	bits = _addfield(layout, sym, typ, bits)

	layout.align = max(layout.align, align)
	layout.size = _size(spec, layout.size, pad + bits)
end

function _addfield(layout::Ctypelayout, ::Type{Calignment{align}}, ::Type{spec}) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}, align}
	pad = padding(strategy(spec), layout.offset, align*8)
	layout.align = max(layout.align, align)
	layout.size = _size(spec, layout.size, pad)
end

@generated function Ctypelayout(::Type{spec}) where {spec<:Ctypespec{<:Any, <:Caggregate, <:Calignment, <:Tuple}}
	layout = Ctypelayout()

	fields = specification(spec)
	while fields !== Tuple{}
		field = Base.tuple_type_head(fields)
		fields = Base.tuple_type_tail(fields)

		layout.offset = _offset(spec, layout.size)
		_addfield(layout, field, spec)
	end

	layout.size += padding(strategy(spec), layout.size, layout.align*8)
	layout.size += -layout.size & (8-1)  # ensure size is divisible by 8

	return layout
end

ext: jl


text:
function rectangle_cell()
    img = zeros(Bool, 200, 200)
    orig_poligon = DSB2018.Model.Area{Bool}()
    for k = 35:45, j = 35 : 45
        img[k, j] = true
        orig_poligon[CartesianIndex(k, j)] = true
    end

    poligon = convert(DSB2018.Model.Area{Bool}, img)

    mdl = DSB2018.Model.XModel(8, 0.2, true, 0.2, 0.9);
    mdl.base_color = 0.8
    mdl.σ = NaN

    x = CartesianIndex(37, 37)
    @assert img[x]

    c = Cell(x, mdl, DSB2018.Model.Area{Bool}());
    DSB2018.Model.Molecular.grow!(c, poligon)

    imshow(c, size(img))
    imshow(img)

    x = collect(keys(orig_poligon))
    y = collect(keys(c.area))

    @show length(x), length(y)
    return (DSB2018.Errors.intersection_over_union(x, y), c)
end

rectangle_cell()


text:
import Stochy: pass_store

@test pass_store(:(function f(); end), :x) == :(function f(x); end)
@test pass_store(:(function f(y); end), :x) == :(function f(x,y); end)
@test pass_store(:(function f(y,z); end), :x) == :(function f(x,y,z); end)

@test pass_store(:(f()), :x) == :(f(x))
@test pass_store(:(f(y)), :x) == :(f(x,y))
@test pass_store(:(f(y,z)), :x) == :(f(x,y,z))

@test expr_eq(pass_store(:(()->false), :x), :((x,)->false))
@test expr_eq(pass_store(:((y)->false), :x), :((x,y)->false))
@test expr_eq(pass_store(:((y,z)->false), :x), :((x,y,z)->false))

@test expr_eq(pass_store(:(f(g())), :x), :(f(x,g(x))))

@test expr_eq(pass_store(:(f()()), :x), :(f(x)(x)))

# Don't pass store to primitives, thunk constructor or trampoline.
@assert :+ in Stochy.primitives
@test expr_eq(pass_store(:(1+2), :x), :(1+2))

expr = :(Stochy.trampoline(Stochy.Thunk(()->false)))
@test expr_eq(pass_store(expr, :x), expr)

# Other syntax.
expr = quote
    if x > 0 || y == :test && a.b < 0
        [true,false][1]
    else
        "none"
    end
end

@test expr_eq(pass_store(expr, :x), expr)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   The NRLMSIS-00 empirical atmosphere model was developed by Mike
#   Picone, Alan Hedin, and Doug Drob based on the MSISE90 model.
#
#   The MSISE90 model describes the neutral temperature and densities in Earth's
#   atmosphere from ground to thermospheric heights. Below 72.5 km the model is
#   primarily based on the MAP Handbook (Labitzke et al., 1985) tabulation of
#   zonal average temperature and pressure by Barnett and Corney, which was also
#   used for the CIRA-86. Below 20 km these data were supplemented with averages
#   from the National Meteorological Center (NMC). In addition, pitot tube,
#   falling sphere, and grenade sounder rocket measurements from 1947 to 1972
#   were taken into consideration. Above 72.5 km MSISE-90 is essentially a
#   revised MSIS-86 model taking into account data derived from space shuttle
#   flights and newer incoherent scatter results. For someone interested only in
#   the thermosphere (above 120 km), the author recommends the MSIS-86 model.
#   MSISE is also not the model of preference for specialized tropospheric work.
#   It is rather for studies that reach across several atmospheric boundaries.
#   (quoted from http://nssdc.gsfc.nasa.gov/space/model/atmos/nrlmsise00.html)
#
#   This Julia version of NRLMSISE-00 was converted from the C version
#   implemented and maintained by Dominik Brodowski <devel@brodo.de> and
#   available at http://www.brodo.de/english/pub/nrlmsise/index.html .
#
#   The source code is available at the following git:
#
#       https://git.linta.de/?p=~brodo/nrlmsise-00.git;a=tree
#
#   The conversion also used information available at the FORTRAN source code
#   available at
#
#       https://ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/nrlmsise00/
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://www.brodo.de/space/nrlmsise/index.html
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export conf_nrlmsise00, gtd7, gtd7d, nrlmsise00

################################################################################
#                                  Constants
################################################################################

include("./nrlmsise00_coefs.jl")

################################################################################
#                                     API
################################################################################

"""
    nrlmsise00(JD::Number, alt::Number, g_lat::Number, g_long::Number [, f107A::Number, f107::Number, ap::Union{Number,AbstractVector}]; output_si::Bool = true, dversion::Bool = true)

**NRLMSISE-00**

Neutral Atmosphere Empirical Model from the surface to lower exosphere.

This routine computes the NRLMSISE-00 outputs (see `NRLMSISE00_Output`) using
the configurations in the structure `nrlmsise00` (see `NRLMSISE00_Structure`).

Notice that the NRLMSISE-00 will be run using the default flags (see
`NRLMSISE00_DEFAULT_FLAGS`). The user can only change the value of
`flags[:output_m_kg]` using the keyword `output_si` to select whether the output
must be converted to SI units. If more control is needed, then the user must
manually call the function `conf_nrlmsise00` and then call `gtd7` or `gtd7d`
with the desired flags.

If the space indices `f107A`, `f107`, and `ap` are missing, then they will be
obtained from the online databases (see `init_space_indices()`).

# Args

* `JD`: Julian Day [UTC].
* `alt`: Altitude [m].
* `g_lat`: Geodetic latitude [rad].
* `g_long`: Geodetic longitude [rad].
* `f107A`: 81 day average of F10.7 flux (centered on day of year `JD`).
* `f107`: Daily F10.7 flux for previous day.
* `ap`: Magnetic index (daily) if it is a number. If it is an array, then see
        **Remarks**.

# Keywords

* `output_si`: (OPTIONAL) If `true`, then the output units will be [m⁻³] for
               species number density and [kg/m⁻³] for the total density.
               Otherwise, the units will be [cm⁻³] and [g/cm⁻³], respectively.
* `dversion`: (OPTIONAL) If `true`, run `gtd7d`. Otherwise, run `gtd7` (see
              **Remarks**).

# Returns

An instance of the structure `NRLMSISE00_Output`. The result in variable
`den_Total` depends on the value of `dversion` (see **Remarks**, **Notes on
input variables**).

# Remarks

1. The densities of `O`, `H`, and `N` are set to `0` below `72.5 km`.
2. The exospheric temperature `T_exo` is set to global average for altitudes
   below `120 km`. The `120 km` gradient is left at global average value for
   altitudes below `72.5 km`.
3. Anomalous oxygen is defined as hot atomic oxygen or ionized oxygen that can
   become appreciable at high altitudes (`> 500 km`) for some ranges of inputs,
   thereby affection drag on satellites and debris. We group these species under
   the term **Anomalous Oxygen**, since their individual variations are not
   presently separable with the drag data used to define this model component.

## AP

If `ap` is a `Vector`, then it must be a vector with 7 dimensions as described
below:

| Index | Description                                                                   |
|-------|:------------------------------------------------------------------------------|
|     1 | Daily AP.                                                                     |
|     2 | 3 hour AP index for current time.                                             |
|     3 | 3 hour AP index for 3 hours before current time.                              |
|     4 | 3 hour AP index for 6 hours before current time.                              |
|     5 | 3 hour AP index for 9 hours before current time.                              |
|     6 | Average of eight 3 hour AP indices from 12 to 33 hours prior to current time. |
|     7 | Average of eight 3 hour AP indices from 36 to 57 hours prior to current time. |

## Notes on input variables

`f107` and `f107A` values used to generate the model correspond to the 10.7 cm
radio flux at the actual distance of the Earth from the Sun rather than the
radio flux at 1 AU. The following site provides both classes of values:

    ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/

`f107`, `f107A`, and `ap` effects are neither large nor well established below
80 km and these parameters should be set to 150, 150, and 4 respectively.

If `dversion` is `true`, then the total mass `den_Total` (see
`NRLMSISE00_Output`) is the sum of the mass densities of the species `He`, `O`,
`N₂`, `O₂`, `Ar`, `H`, and `N`, but **does not** include anomalous oxygen.

If `dversion` is `false`, then total mass `den_Total` (see `NRLMSISE00_Output`)
is the effective total mass density for drag and is the sum of the mass
densities of all species in this model **including** the anomalous oxygen.

"""
function nrlmsise00(JD::Number, alt::Number, g_lat::Number, g_long::Number;
                    output_si::Bool = true, dversion::Bool = true)

    # Get the space indices. If the altitude is lower than 80km, set to default
    # according to the instructions in NRLMSISE-00 source code.
    if alt < 80e3
        f107A = 150.0
        f107  = 150.0
        ap    = 4.0
    else
        # TODO: The online version of NRLMSISE-00 seems to use 90 days, whereas
        # the NRLMSISE-00 source code mentions 81 days.
        f107A = get_space_index(Val(:F10Madj), JD; window = 90)
        f107  = get_space_index(Val(:F10adj),  JD-1)
        ap    = get_space_index(Val(:Ap),      JD)
    end

    # Call the NRLMSISE-00 model.
    nrlmsise00(JD, alt, g_lat, g_long, f107A, f107, ap; output_si = output_si,
               dversion = dversion)
end

function nrlmsise00(JD::Number, alt::Number, g_lat::Number, g_long::Number,
                    f107A::Number, f107::Number,
                    ap::Union{Number,AbstractVector}; output_si::Bool = true,
                    dversion::Bool = true)

    # Constant
    dgtr = 1.74533e-2 # Convert degrees to radians.

    # Convert the Julian Day to Date.
    (Y, M, D, h, m, s) = jd_to_date(JD)

    # Get the number of days since the beginning of the year.
    doy = round(Int,date_to_jd(Y, M, D, 0, 0, 0) - date_to_jd(Y, 1, 1, 0, 0, 0)) + 1

    # Get the number of seconds since the beginning of the day.
    Δds = 3600h + 60m + s

    # Get the local apparent solar time [hours].
    #
    # TODO: To be very precise, I think this should also take into consideration
    # the Equation of Time. However, the online version of NRLMSISE-00 does not
    # use this.
    lst = Δds/3600 + g_long*12/π

    # Create the input structure for NRLMSISE-00 converting the arguments.
    #
    # Notice that, for this algorithm, the year **is not** used.
    nrlmsise00d = conf_nrlmsise00(Y, doy, Δds, alt/1000, g_lat/dgtr,
                                  g_long/dgtr, lst, f107A, f107, ap,
                                  NRLMSISE00_Flags(output_m_kg = output_si))

    # Call the NRLMSISE-00 model.
    nrlmsise00_out = (dversion) ? gtd7d(nrlmsise00d) : gtd7(nrlmsise00d)

    nrlmsise00_out
end

################################################################################
#                             Auxiliary Functions
################################################################################

"""
    conf_nrlmsise00(year::Int, doy::Int, sec::Number, alt::Number, g_lat::Number, g_long::Number, lst::Number, f107A::Number, f107::Number, ap::[Number, AbstractVector], flags::NRLMSISE00_Flags = NRLMSISE00_Flags())

Create the structure with the proper configuration to call the NRLMSISE-00
model.

Notice that the input variables have the same units of the original model.

# Args

* `year`: Year (currently ignored).
* `doy`: Day of year.
* `sec`: Seconds in day [UT].
* `alt`: Altitude [km].
* `g_lat`: Geodetic latitude [deg].
* `g_long`: Geodetic longitude [deg].
* `lst`: Local apparent solar time (hours).
* `f107A`: 81 day average of F10.7 flux (centered on day of year `doy`).
* `f107`: Daily F10.7 flux for previous day.
* `ap`: Magnetic index (daily) if it is a number. If it is an array, then see
        **Remarks**.
* `flags`: (OPTIONAL) An instance of the structure `NRLMSISE00_Flags` with the
            configuration flags for NRLMSISE00. If omitted, then the default
            configurations will be used.

# Returns

An instance of the structure `NRLMSISE00_Structure`.

# Remarks

If `ap` is a `Vector`, then it must be a vector with 7 dimensions as described
below:

| Index | Description                                                                   |
|-------|:------------------------------------------------------------------------------|
|     1 | Daily AP.                                                                     |
|     2 | 3 hour AP index for current time.                                             |
|     3 | 3 hour AP index for 3 hours before current time.                              |
|     4 | 3 hour AP index for 6 hours before current time.                              |
|     5 | 3 hour AP index for 9 hours before current time.                              |
|     6 | Average of eight 3 hour AP indices from 12 to 33 hours prior to current time. |
|     7 | Average of eight 3 hour AP indices from 36 to 57 hours prior to current time. |

## Notes on input variables

UT, Local Time, and Longitude are used independently in the model and are not of
equal importance for every situation. For the most physically realistic
calculation these three variables should be consistent (`lst=sec/3600 +
g_long/15`). The Equation of Time departures from the above formula for
apparent local time can be included if available but are of minor importance.

`f107` and `f107A` values used to generate the model correspond to the 10.7 cm
radio flux at the actual distance of the Earth from the Sun rather than the
radio flux at 1 AU. The following site provides both classes of values:

    ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/

`f107`, `f107A`, and `ap` effects are neither large nor well established below
80 km and these parameters should be set to 150, 150, and 4 respectively.

"""
function conf_nrlmsise00(year::Int,
                         doy::Int,
                         sec::Number,
                         alt::Number,
                         g_lat::Number,
                         g_long::Number,
                         lst::Number,
                         f107A::Number,
                         f107::Number,
                         ap::Number,
                         flags::NRLMSISE00_Flags = NRLMSISE00_Flags())

    # Constants
    # =========

    dgtr = 1.74533e-2 # Convert degrees to radians.
    hr   = 0.2618     # Convert hour angle to radians.

    # Compute auxiliary variables
    # ===========================
    df  = f107  - f107A
    dfa = f107A - 150

    stloc, ctloc  = sincos(1hr*lst)
    s2tloc,c2tloc = sincos(2hr*lst)
    s3tloc,c3tloc = sincos(3hr*lst)

    # Compute Legendre polynomials.
    #
    # TODO: Not all coefficients are used. Hence, we could gain performance
    # here.
    plg = legendre(Val(:conv), pi/2-g_lat*dgtr, 8, 8, false)'

    # Latitude variation of gravity
    # =============================
    #
    # None for flags.time_independent = false.
    xlat      = (!flags.time_independent) ? 45.0 : Float64(g_lat)
    gsurf, re = _glatf(xlat)

    # Create and return the structure
    # ===============================

    # This is necessary to avoid type instability as reported here:
    #
    #   https://github.com/mauro3/Parameters.jl/issues/58
    #
    nrlmsise00d::NRLMSISE00_Structure{Float64} =
        NRLMSISE00_Structure{Float64}(year     = year,
                                      doy      = doy,
                                      sec      = sec,
                                      alt      = alt,
                                      g_lat    = g_lat,
                                      g_long   = g_long,
                                      lst      = lst,
                                      f107A    = f107A,
                                      f107     = f107,
                                      ap       = ap,
                                      flags    = flags,
                                      re       = re,
                                      gsurf    = gsurf,
                                      df       = df,
                                      dfa      = dfa,
                                      plg      = plg,
                                      stloc    = stloc,
                                      ctloc    = ctloc,
                                      s2tloc   = s2tloc,
                                      c2tloc   = c2tloc,
                                      s3tloc   = s3tloc,
                                      c3tloc   = c3tloc)

    nrlmsise00d
end

function conf_nrlmsise00(year::Int, doy::Int, sec::Number, alt::Number,
                         g_lat::Number, g_long::Number, lst::Number,
                         f107A::Number, f107::Number, ap::AbstractVector,
                         flags::NRLMSISE00_Flags = NRLMSISE00_Flags())

    # Constants
    # =========

    dgtr = 1.74533e-2 # Convert degrees to radians.
    hr   = 0.2618     # Convert hour angle to radians.

    # Compute the new flags
    # =====================
    flags.use_ap_array = true

    # Compute auxiliary variables
    # ===========================
    df     = f107  - f107A
    dfa    = f107A - 150

    stloc, ctloc  = sincos(1hr*lst)
    s2tloc,c2tloc = sincos(2hr*lst)
    s3tloc,c3tloc = sincos(3hr*lst)

    # Compute Legendre polynomials.
    #
    # TODO: Not all coefficients are used. Hence, we could gain performance
    # here.
    plg = legendre(Val(:conv), pi/2-g_lat*dgtr, 8, 8, false)'

    # Latitude variation of gravity
    # =============================
    #
    # None for flags.time_independent = false.
    xlat      = (!flags.time_independent) ? 45.0 : Float64(g_lat)
    gsurf, re = _glatf(xlat)

    # Create and return the structure
    # ===============================

    # This is necessary to avoid type instability as reported here:
    #
    #   https://github.com/mauro3/Parameters.jl/issues/58
    #
    nrlmsise00d::NRLMSISE00_Structure{Float64} =
        NRLMSISE00_Structure{Float64}(year     = year,
                                      doy      = doy,
                                      sec      = sec,
                                      alt      = alt,
                                      g_lat    = g_lat,
                                      g_long   = g_long,
                                      lst      = lst,
                                      f107A    = f107A,
                                      f107     = f107,
                                      ap_array = ap,
                                      flags    = flags,
                                      re       = re,
                                      gsurf    = gsurf,
                                      df       = df,
                                      dfa      = dfa,
                                      plg      = plg,
                                      stloc    = stloc,
                                      ctloc    = ctloc,
                                      s2tloc   = s2tloc,
                                      c2tloc   = c2tloc,
                                      s3tloc   = s3tloc,
                                      c3tloc   = c3tloc)

    nrlmsise00d
end

################################################################################
#                               Public Functions
################################################################################

"""
    gtd7(nrlmsise00d::NRLMSISE00_Structure{T}) where T<:Number

**NRLMSISE-00**

Neutral Atmosphere Empirical Model from the surface to lower exosphere.

This routine computes the NRLMSISE-00 outputs (see `NRLMSISE00_Output`) using
the configurations in the structure `nrlmsise00` (see `NRLMSISE00_Structure`).

# Args

* `nrlmsise00d`: An instance of `NRLMSISE00_Structure`.

# Returns

An instance of structure `NRLMSISE00_Output` with the outputs.

In this case, the total mass `den_Total` (see `NRLMSISE00_Output`) is the sum of
the mass densities of the species `He`, `O`, `N₂`, `O₂`, `Ar`, `H`, and `N`, but
**does not** include anomalous oxygen.

# Remarks

1. The densities of `O`, `H`, and `N` are set to `0` below `72.5 km`.
2. The exospheric temperature `T_exo` is set to global average for altitudes
   below `120 km`. The `120 km` gradient is left at global average value for
   altitudes below `72.5 km`.
3. Anomalous oxygen is defined as hot atomic oxygen or ionized oxygen that can
   become appreciable at high altitudes (`> 500 km`) for some ranges of inputs,
   thereby affection drag on satellites and debris. We group these species under
   the term **Anomalous Oxygen**, since their individual variations are not
   presently separable with the drag data used to define this model component.

"""
function gtd7(nrlmsise00d::NRLMSISE00_Structure{T}) where T<:Number

    @unpack_NRLMSISE00_Structure nrlmsise00d

    # Constants
    # =========

    mn2  = 4
    zn2  = SVector{4,T}([72.5,55.0,45.0,32.5])
    mn3  = 5
    zn3  = SVector{5,T}([32.5,20.0,15.0,10.0,0.0])
    zmix = T(62.5)

    # Initialization of variables
    # ===========================
    meso_tn2  = MVector{4,T}(undef)
    meso_tn3  = MVector{5,T}(undef)
    meso_tgn2 = MVector{2,T}(undef)
    meso_tgn3 = MVector{2,T}(undef)

    @inbounds begin
        # Latitude variation of gravity
        # =============================

        xmm = pdm[3,5]

        # Thermosphere / Mesosphere (above zn2[1])
        # ========================================

        (alt < zn2[1]) && (nrlmsise00d.alt = zn2[1])

        out_thermo = gts7(nrlmsise00d)
        nrlmsise00d.alt = alt

        # Unpack the values again because `gts7` may have modified `nrlmsise00d`.
        @unpack_NRLMSISE00_Structure nrlmsise00d

        # If we are above `zn2[1]`, then we do not need to compute anything else.
        (alt >= zn2[1]) && return out_thermo

        # Unpack the output values from the thermospheric portion.
        @unpack_NRLMSISE00_Output out_thermo

        # Check if we must convert to SI.
        dm28m = (flags.output_m_kg) ?  dm28*1e6 : dm28

        # Lower Mesosphere / Upper Stratosphere (between `zn3[1]` and `zn2[1]`)
        # =====================================================================

        meso_tgn2[1] = meso_tgn1_2

        meso_tn2[1]  = meso_tn1_5

        meso_tn2[2]  =  pma[1,1]*pavgm[1]/
                        ( 1 - flags.all_tn2_var*_glob7s(pma[1,:], nrlmsise00d) )
        meso_tn2[3]  =  pma[2,1]*pavgm[2]/
                        ( 1 - flags.all_tn2_var*_glob7s(pma[2,:], nrlmsise00d) )
        meso_tn2[4]  =  pma[3,1]*pavgm[3]/
                        ( 1 - flags.all_tn2_var*flags.all_tn3_var*_glob7s(pma[3,:], nrlmsise00d) )
        meso_tn3[1]  = meso_tn2[4]

        meso_tgn2[2] = pavgm[9]*pma[10,1]*
                        ( 1 + flags.all_tn2_var*flags.all_tn3_var*_glob7s(pma[10,:], nrlmsise00d))*
                        meso_tn2[4]^2/(pma[3,1]*pavgm[3])^2

        # Lower Stratosphere and Troposphere (below `zn3[1]`)
        # ===================================================

        if alt < zn3[1]
            meso_tgn3[1] = meso_tgn2[2]
	    	meso_tn3[2]  = pma[4,1]*pavgm[4]/
                           ( 1 - flags.all_tn3_var*_glob7s(pma[4,:], nrlmsise00d) )
	    	meso_tn3[3]  = pma[5,1]*pavgm[5]/
                           ( 1 - flags.all_tn3_var*_glob7s(pma[5,:], nrlmsise00d) )
	    	meso_tn3[4]  = pma[6,1]*pavgm[6]/
                           ( 1 - flags.all_tn3_var*_glob7s(pma[6,:], nrlmsise00d) )
	    	meso_tn3[5]  = pma[7,1]*pavgm[7]/
                           ( 1 - flags.all_tn3_var*_glob7s(pma[7,:], nrlmsise00d) )
	    	meso_tgn3[2] = pma[8,1]*pavgm[8]*
                           ( 1 + flags.all_tn3_var*_glob7s(pma[8,:], nrlmsise00d))*
                           meso_tn3[5]*meso_tn3[5]/(pma[7,1]*pavgm[7])^2
        end

        # Linear Transition to Full Mixing Below `zn2[1]`
        # ===============================================

        dmc  = (alt > zmix) ? 1 - (zn2[1]-alt)/(zn2[1] - zmix) : T(0)
        dz28 = den_N2

        # N2 Density
        # ==========

        dmr = den_N2 / dm28m - 1

        den_N2, tz = _densm(re, gsurf, alt, dm28m, xmm, T(0), zn3, meso_tn3,
                            meso_tgn3, zn2, meso_tn2, meso_tgn2)

        den_N2 *= 1 + dmr*dmc

        # He Density
        # ==========

        dmr    = den_He / (dz28 * pdm[1,2]) - 1
        den_He = den_N2 * pdm[1,2] * ( 1 + dmr*dmc )

        # O Density
        # =========

        den_O  = T(0)
        den_aO = T(0)

        # O2 Density
        # ==========

        dmr    = den_O2 / (dz28 * pdm[4,2]) - 1
        den_O2 = den_N2 * pdm[4,2] * ( 1 + dmr*dmc )

        # Ar Density
        # ==========

        dmr    = den_Ar / (dz28 * pdm[5,2]) - 1
        den_Ar = den_N2 * pdm[5,2] * ( 1 + dmr*dmc )

        # H Density
        # =========

        den_H = T(0)

        # N Density
        # =========

        den_N = T(0)

        # Total Mass Density
        # ==================

        den_Total = 1.66e-24( 4den_He +
                             16den_O  +
                             28den_N2 +
                             32den_O2 +
                             40den_Ar +
                               den_H  +
                             14den_N)

        # Check if we must convert the units to SI.
        (flags.output_m_kg) && (den_Total /= 1000)

        # Temperature at Selected Altitude
        # ================================

        dd, T_alt = _densm(re, gsurf, alt, T(1), T(0), tz, zn3, meso_tn3,
                           meso_tgn3, zn2, meso_tn2, meso_tgn2)

        # Create output structure and return.
        #
        # This is necessary to avoid type instability as reported here:
        #
        #   https://github.com/mauro3/Parameters.jl/issues/58
        #

        nrlmsise00_out::NRLMSISE00_Output{T} =
            NRLMSISE00_Output{T}(den_N     = den_N,
                                 den_N2    = den_N2,
                                 den_O     = den_O,
                                 den_aO    = den_aO,
                                 den_O2    = den_O2,
                                 den_H     = den_H,
                                 den_He    = den_He,
                                 den_Ar    = den_Ar,
                                 den_Total = den_Total,
                                 T_exo     = T_exo,
                                 T_alt     = T_alt,
                                 flags     = flags)
    end

    nrlmsise00_out
end

"""
    gtd7d(nrlmsise00d::NRLMSISE00_Structure{T}) where T<:Number

**NRLMSISE-00**

Neutral Atmosphere Empirical Model from the surface to lower exosphere.

This routine computes the NRLMSISE-00 outputs (see `NRLMSISE00_Output`) using
the configurations in the structure `nrlmsise00` (see `NRLMSISE00_Structure`).

# Args

* `nrlmsise00d`: An instance of `NRLMSISE00_Structure`.

# Returns

An instance of structure `NRLMSISE00_Output` with the outputs.

In this case, the total mass `den_Total` (see `NRLMSISE00_Output`) is the
effective total mass density for drag and is the sum of the mass densities of
all species in this model **including** the anomalous oxygen.

# Remarks

1. The densities of `O`, `H`, and `N` are set to `0` below `72.5 km`.
2. The exospheric temperature `T_exo` is set to global average for altitudes
   below `120 km`. The `120 km` gradient is left at global average value for
   altitudes below `72.5 km`.
3. Anomalous oxygen is defined as hot atomic oxygen or ionized oxygen that can
   become appreciable at high altitudes (`> 500 km`) for some ranges of inputs,
   thereby affection drag on satellites and debris. We group these species under
   the term **Anomalous Oxygen**, since their individual variations are not
   presently separable with the drag data used to define this model component.

"""
function gtd7d(nrlmsise00d::NRLMSISE00_Structure{T}) where T<:Number
    # Call `gt7d` to compute the NRLMSISE-00 outputs.
    out = gtd7(nrlmsise00d)

    # Update the computation of the total mass density.
    den_Total = 1.66e-24( 4out.den_He +
                         16out.den_O  +
                         28out.den_N2 +
                         32out.den_O2 +
                         40out.den_Ar +
                           out.den_H  +
                         14out.den_N  +
                         16out.den_aO)

    # Check if we must convert the units to SI.
    (nrlmsise00d.flags.output_m_kg) && (den_Total /= 1000)

    # Create the new output and return.
    #
    # This is necessary to avoid type instability as reported here:
    #
    #   https://github.com/mauro3/Parameters.jl/issues/58
    #

    # Densities.
    #
    # Originally, we were using `reconstruct` to copy all values except one:
    #
    #   NRLMSISE00_Output(out; den_Total = den_Total)
    #
    # However, it turns out to be a very slow method. Hence, it was change to
    # this version.
    nrlmsise00_out::NRLMSISE00_Output{T} =
        NRLMSISE00_Output(den_N     = out.den_N,
                          den_N2    = out.den_N2,
                          den_O     = out.den_O,
                          den_aO    = out.den_aO,
                          den_O2    = out.den_O2,
                          den_H     = out.den_H,
                          den_He    = out.den_He,
                          den_Ar    = out.den_Ar,
                          den_Total = den_Total,
                          T_exo     = out.T_exo,
                          T_alt     = out.T_alt,
                          flags     = out.flags)

    nrlmsise00_out
end

################################################################################
#                              Private Functions
################################################################################

include("./nrlmsise00_priv.jl")


text:
"""
    LRFinderResult(lrs, losses)

Result of the learning rate finder [`lrfind`](#). Use `plot`
to visualize.
"""
struct LRFinderResult
    losses
    lrs
    estimators
    estimates
end

"""
    lrfind(learner[, dataiter; kwargs...]) -> LRFinderResult

Run the learning rate finder. Exponentially increases the learning rate from a very
low value to a very high value and uses the losses to estimate an optimal learning
rate. Return a [`LRFinderResult`](#).

## Keyword arguments

- `nsteps = 100`: maximum number of steps to run the learning rate finder for
- `startlr = 1e-7`: minimum learning rate
- `endlr = 10`: maximum learning rate
- `divergefactor`: stop finder early if loss goes higher than lowest loss times
    this factor
- `estimators = [Steepest(), MinDivByTen()]`: list of [`LREstimator`](#)s
"""
function lrfind(
        learner,
        dataiter = learner.data.training;
        nsteps = 100,
        startlr = 1e-7,
        endlr = 10,
        divergefactor = 4,
        estimators = [Steepest(), MinDivByTen()])
    losses = Float64[]
    lrs = Float64[]
    bestloss = Inf
    scheduler = FluxTraining.removecallback!(learner, Scheduler)  # remove current `Scheduler` so it does not interfere
    modelcheckpoint = deepcopy(cpu(learner.model))

    withfields(
        learner,
        model = modelcheckpoint,
        params = params(modelcheckpoint),
        optimizer = deepcopy(learner.optimizer)
        ) do

        FluxTraining.runepoch(learner, TrainingPhase()) do _
            for (i, batch) in zip(1:nsteps, dataiter)
                lr = startlr * (endlr / startlr) ^ (i / nsteps)
                learner.optimizer.eta = lr

                state = step!(learner, TrainingPhase(), batch)

                push!(losses, state.loss)
                push!(lrs, learner.optimizer.eta)
                bestloss = min(state.loss, bestloss)
                if state.loss > bestloss * divergefactor
                    throw(FluxTraining.CancelEpochException("Learning rate finder diverged."))
                end
            end
        end
    end

    # Restore previous `Learner` state
    isnothing(scheduler) || FluxTraining.addcallback!(learner, scheduler)
    FluxTraining.model!(learner, modelcheckpoint)
    return LRFinderResult(losses, lrs, estimators)
end

# `LREstimator`s allow giving back different suggestions for a learning rate.

"""
    abstract type LREstimator

Estimator for an optimal learning rate. Needs to implement [`estimatelr`](#).

See [`Steepest`](#) and [`MinDivByTen`](#).
"""
abstract type LREstimator end

"""
    estimatelr(::LREstimator, losses, lrs)

Estimate the optimal learning rate using `losses` and `lrs`.
"""
function estimatelr end

"""
    Steepest <: LREstimator

Estimate the optimal learning rate to be where the gradient of the loss
is the steepest, i.e. the decrease is largest.
"""
struct Steepest <: LREstimator
    beta
end
Steepest() = Steepest(0.98)

function estimatelr(est::Steepest, losses, lrs)
    slosses = smoothvalues(losses, est.beta)
    grads = (slosses[2:end] .- slosses[1:end-1]) ./ log.(lrs[2:end] .- lrs[1:end-1])
    i = length(lrs) ÷ 3
    lr = lrs[i:end][argmax(grads[i:end])]
    return lr
end

"""
    MinDivByTen <: LREstimator

Estimate the optimal learning rate to be value at the minimum loss divided by 10.
"""
struct MinDivByTen <: LREstimator
    beta
end
MinDivByTen() = MinDivByTen(0.98)

function estimatelr(est::MinDivByTen, losses, lrs)
    i = length(lrs) ÷ 3
    lr = lrs[i:end][argmin(smoothvalues(losses, est.beta)[i:end])] / 10
    return lr
end

LRFinderResult(losses, lrs) = LRFinderResult(losses, lrs, [Steepest(), MinDivByTen()])
LRFinderResult(losses, lrs, estimators) = LRFinderResult(
    losses,
    lrs,
    estimators,
    [estimatelr(est, losses, lrs) for est in estimators]
)

# Printing and plotting

function Base.show(io::IO, lrfindresult::LRFinderResult)
    println(io, "LRFindResult(")
    for (est, v) in zip(lrfindresult.estimators, lrfindresult.estimates)
        println(io, "    ", est, " => ", v)
    end
    print(io, ")")
end

function Makie.plot(result::LRFinderResult)
    ticks = [round((10.)^i, digits=abs(i)) for i in -10:2]
    fig = Figure()
    ax = Axis(
        title = "Learning rate finder",
        titlesize = 20,
        fig[1, 1],
        xscale = log,
        xticks = (ticks, string.(ticks)),
        xminorticks = IntervalsBetween(5),
        xminorgridvisible=true,
        ygridcolor = :white,
        xlabelsize = 14,
        ylabelsize = 14,
        ylabelcolor = :gray,
        xtickcolor= :gray,
        xticklabelcolor= :black,
        xticklabelsize= 12,
        ytickcolor= :gray,
        yticklabelcolor= :gray,
        yticklabelsize= 12,
        ylabel = "Loss",
        xlabel = "Learning rate (log)")

    lines!(
        result.lrs,
        smoothvalues(result.losses, 0.98),
        color = :black,
    )

    hidespines!(ax)

    # plot suggestions
    ls = []
    for (estim, val) in zip(result.estimators, result.estimates)
        push!(ls, vlines!(ax, [val]))

    end

    leg = Legend(
        fig[2, 1], ls,
        ["$(estim): $(round(val, sigdigits=3))" for (estim, val) in zip(result.estimators, result.estimates)],
        framevisible=false,
        labelsize=14,
        orientation = :horizontal,)
    leg.tellheight[] = true
    leg.tellwidth[] = false
    fig
end

# Utilities

"""
    smoothvalues(xs, β)

Apply exponential smoothing with parameter `b` to vector `xs`.
"""
function smoothvalues(xs, β)
    res = similar(xs)
    val = 0
    for (i, x) in enumerate(xs)
        val = (val * β) + (x * (1 - β))
        res[i] = val/((1)-β^i)
    end
    res
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "broyden" begin
    dep=nep_gallery("dep1");
    S,V=broyden(dep)
    # Broyden returns a Schur factorization so check with MM
    @test opnorm(compute_MM(dep,S,V))<sqrt(eps())
    # test addconj
    S,V=broyden(dep,addconj=true,pmax=5)
    # Test by computing the eigenpairs
    λv,X = eigen(S)
    V_eigvecs=V*X;
    verify_lambdas(4, dep, λv, V_eigvecs, sqrt(eps()))
end


text:
# Julia wrapper for header: /usr/local/include/libavcodec/avdct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avcodec_dct_alloc,
    avcodec_dct_init,
    avcodec_dct_get_class

function avcodec_dct_alloc()
    ccall((:avcodec_dct_alloc,libavcodec),Ptr{AVDCT},())
end

function avcodec_dct_init(arg1)
    ccall((:avcodec_dct_init,libavcodec),Cint,(Ptr{AVDCT},),arg1)
end

function avcodec_dct_get_class()
    ccall((:avcodec_dct_get_class,libavcodec),Ptr{AVClass},())
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
function α(s::Union{AbstractVector{S},S},K::Symmetric{T}) where {S<:Integer,T<:Real}
    N = size(K,2)
    s = broadcast(min,s,N) # any value s>N will get an α score of 1
    eig = reverse(eigvals(K, N-maximum(s)+1:N)) # get the largest eigenvalues (as many as we need) and sort them in decreasing order
    num = cumsum(eig)[s] # ∑ᵢ₌₁ˢ λᵢ² ∀s
    denom = tr(K) # ∑_ᵢ₌₁ⁿ λᵢ²
    broadcast(sqrt,num/denom)
end
α(s::Union{AbstractVector{S},S}, X::AbstractMatrix) where {S<:Integer} = α(s,Symmetric(X'*X))

_resamplerow(x) = rand(x,size(x))
_resamplematrix(X) = mapslices(_resamplerow, X, dims=2)

function _αbootstrap(X::AbstractMatrix, s::Union{AbstractVector{S},S}, nbrIter::Integer) where {S<:Integer}
    sum(i->α(s,_resamplematrix(X)), 1:nbrIter) # bootstrap and sum
end

function projectionscore(X::AbstractMatrix, s::Union{AbstractVector{S},S}; nbrIter::Integer=10) where {S<:Integer}
    αB = NaN
    if nprocs()==1 || nbrIter==1
        αB=_αbootstrap(X,s,nbrIter)/nbrIter # no threading needed
    else
        # try to minimize data movement by spawning once per process
        W = workers()
        nbrWorkers = min( length(W), nbrIter ) # never use more than nbrIter workers
        nbrIterPerWorker = [ ((d,r)=divrem(nbrIter,nbrWorkers); d+(r>=i)) for i=1:nbrWorkers ] # divide iterations evenly among workers

        # TODO: do with pmap instead?
        refs = Array{Any}(nbrWorkers)
        for i=1:nbrWorkers
            refs[i] = @spawnat W[i] _αbootstrap(X,s,nbrIterPerWorker[i])
        end
        αB = sum(fetch, refs)/nbrIter # take mean over all bootstrapped α's
    end
    α(s,X) - αB
end

# Assumes X is sorted by decreasing σ!
function _αfiltered(X, s, σ, σThresholds)
    αs = zeros(length(s), length(σThresholds))

    # Build Kernel Matrix gradually, starting with the variables with highest σ
    K = zeros(size(X,2),size(X,2))
    prevInd = 1
    for i = length(σThresholds):-1:1
        # get range of variables to add
        currInd = prevInd
        while currInd<length(σ) && σ[currInd] >= σThresholds[i]
            currInd += 1
        end

        varRange = prevInd:currInd-1

        if !isempty(varRange)
            X2 = X[varRange,:]
            K += X2'X2
            αs[:,i] = α(s, Symmetric(K))
        elseif i<length(σThresholds)
            αs[:,i] = αs[:,i+1]
        end
        prevInd = currInd
    end

    αs
end

function _αfilteredsum(X,s,σ,σThresholds,nbrIter)
    sum(i->_αfiltered(_resamplematrix(X),s,σ,σThresholds), 1:nbrIter)
end

_stdnormalized(X) = (σ=dropdims(std(X,dims=2); dims=2); σ/maximum(σ))

function projectionscorefiltered(X::AbstractMatrix, s::Union{AbstractVector{S},S}, σThresholds::AbstractVector;
                                             nbrIter::Integer=10, σ=_stdnormalized(X)) where {S<:Integer}
    @assert all(diff(σThresholds).>=0) "σThresholds must be increasing"
    @assert length(σ) == size(X,1)

    # get rid of those variables that will not be used at any threshold
    baseFilter = σ.>=σThresholds[1]
    X = X[baseFilter,:]
    σ = σ[baseFilter]

    if length(σ)==0
        @warn "No variables remain after filtering at lowest given threshold."
        return zeros(length(s), length(σThresholds))
    end

    # sort X so that the variables are ordered by decreasing σ
    σPerm = sortperm(σ, rev=true)
    σ = σ[σPerm]
    X = X[σPerm,:]

    α = _αfiltered(X,s,σ,σThresholds)

    αB = NaN

    if nprocs()==1 || nbrIter==1  # no threading needed
        αB = _αfilteredsum(X,s,σ,σThresholds,nbrIter)
    else
        # try to minimize data movement by spawning once per process
        W = workers()
        nbrWorkers = min( length(W), nbrIter ) # never use more than nbrIter workers
        nbrIterPerWorker = [ ((d,r)=divrem(nbrIter,nbrWorkers); d+(r>=i)) for i=1:nbrWorkers ] # divide iterations evenly among workers

        refs = Array{Any}(undef,nbrWorkers)
        for i=1:nbrWorkers
            refs[i] = @spawnat W[i] _αfilteredsum(X,s,σ,σThresholds,nbrIterPerWorker[i])
        end
        αB = sum(fetch, refs) # take mean over all bootstrapped α's
    end

    α - αB/nbrIter
end


text:
# This script provides an example of simulating changepoint data,
# using PELT to find these changepoints, and plotting the results
#
# Requires Plots package

using Distributions, Plots
using Changepoints

#########################
# Simulate changepoints #
#########################

n = 1000        # Length of sample
λ = 100         # Frequencey of changes

# Mu parameter sampled from a N(0,10) distribution for each changepoint
μ = Normal(0,5)
σ = 1.0 # Sigma parameter fixed

# Construct sample with changepoints where
# at each changepoint a new value is sample for μ
sample, cps = @changepoint_sampler n λ Normal(μ, σ)

########################
# Finding changepoints #
########################

# Find changepoints via PELT
# assuming Normal likelihood with changing mean
# parameter and fixed standard deviation of σ
pelt_cps, pelt_cost = @PELT sample Normal(:?, σ)

# Find changepoints via Binary Segmentation
bs_cps, bs_cost = @BS sample Normal(:?, σ)

# If plotting multiple sets of changepoints which need to labelled
# separately, it is best not to use `changepoint_plot` function
# but to instead first plot data
# and then add changepoints with `vline!` function
p = plot(sample, label="Simulated data")
vline!(p, cps, label="True changepoints")
vline!(p, pelt_cps, label="PELT changepoints")
vline!(p, bs_cps, label="BS changepoints")
display(p)


text:
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk
using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")
x = y = z = [-2, -1, 0, 1, 2];
rect = VTKRectilinearData((x, y, z));
rect.cell_data["Cell scalar"] = reshape([rand() for i in 1:num_of_cells(rect)], cell_extents(rect));
unstruct = VTKUnstructuredData(rect)
pyunstruct = PyVTK(unstruct)
dash_vtk_utils = pyimport("dash_vtk.utils")

dash_vtk_utils.to_mesh_state
state = dash_vtk_utils.to_mesh_state(pyunstruct)

content = vtk_view([
        vtk_geometryrepresentation([
        vtk_mesh(state=state)
    ]),
])

app = dash()

app.layout = html_div() do
    html_div(
    style = Dict(
        "width" => "100%",
        "height" => "400px"
    ),
    children=[content],
    )
end

run_server(app, "0.0.0.0", debug = true)


text:
mutable struct RobustMaximinObjectiveCore
    points::Vector{Float64}
    standard_deviation::Array{Float64,2}
    IIF::Array{Float64,2}
    RobustMaximinObjectiveCore() =
        new(Vector{Float64}(), Array{Float64,2}(undef, 0, 0), Array{Float64,2}(undef, 0, 0))
end


text:
function controller(t, z, walker)
    tt = t - walker.t₀
    θ₁ = z[1]
    θ₂ = z[3]
    ω₁ = z[2]
    ω₂ = z[4]
    M = walker.M; m = walker.m; I=walker.I; l = walker.l; c = walker.c;
    g = walker.g; γ = walker.γ

    t₀ = 0
    tf = walker.tf
    q₀ = walker.Θ₂
    qf = walker.Θf
    q₀̇ = walker.Θ₂̇
    qḟ = 0.0

    AA = [  1 t₀ t₀^2   t₀^3    t₀^4    t₀^5;
            1 tf tf^2   tf^3    tf^4    tf^5;
            0  1 2*t₀ 3*t₀^2  4*t₀^3  5*t₀^4;
            0  1 2*tf 3*tf^2  4*tf^3  5*tf^4;
            0  0    2   6*t₀ 12*t₀^2 20*t₀^3;
            0  0    2   6*tf 12*tf^2 20*tf^3]
    bb = [q₀; qf; q₀̇ ; qḟ; 0.; 0.]
    xx = AA\bb
    a0 = xx[1]; a1 = xx[2]; a2 = xx[3];
    a3 = xx[4]; a4 = xx[5]; a5 = xx[6];
    if (tt>tf)
        tt = tf
    end

    θ₂ᵣ = a5*tt^5 + a4*tt^4 + a3*tt^3 + a2*tt^2 + a1*tt + a0
    θ₂̇ᵣ = 5*a5*tt^4 + 4*a4*tt^3 + 3*a3*tt^2 + 2*a2*tt + a1
    θ₂̈ᵣ =   20*a5*tt^3 + 12*a4*tt^2 + 6*a3*tt + 2*a2

    if t >= 4.5
        # println("t = ",t)
        θ₂ᵣ = 2*θ₂ᵣ
    end

    #partial feedback linearization
    M11 = 2*I + M*l^2 + 2*c^2*m + 2*l^2*m - 2*c*l*m - 2*c*l*m*cos(θ₂)
    M12 = I + c^2*m - c*l*m*cos(θ₂)
    M21 = I + c^2*m - c*l*m*cos(θ₂)
    M22 = I + c^2*m
    Ms = [M11 M12; M21 M22]

    N1 = c*g*m*sin(θ₁ - γ + θ₂) + M*g*l*sin(γ - θ₁) - c*g*m*sin(γ - θ₁) + 2*g*l*m*sin(γ - θ₁) + c*l*m*ω₂^2*sin(θ₂) + 2*c*l*m*ω₁*ω₂*sin(θ₂)
    N2 = c*m*(g*sin(θ₁ - γ + θ₂) - l*ω₁^2*sin(θ₂))
    Ns = [N1; N2]

    Kp = walker.Kp
    Kd = 2*sqrt(Kp)
    B = [0; 1]
    Sc = [0 1]
    e = θ₂-θ₂ᵣ
    ė = ω₂ - θ₂̇ᵣ;
    v = [θ₂̈ᵣ - Kp*e - Kd*ė]
    Minv = Ms\[1. 0.; 0. 1.]
    u = (Sc*Minv*B)\(v+Sc*Minv*Ns)
    println("torque = ",u)

    return (u, θ₂ᵣ, θ₂̇ᵣ, θ₂̈ᵣ)

end


text:
import Base.mean

export mean, mean_rb

# Mean of MCMCChain
mean(c::MCMCChain, par::Ranges=1:ncol(c.samples)) = Float64[mean(c.samples[:, i]) for i = par]

mean(c::MCMCChain, par::Real) = mean(c, par:par)

# Mean of MCMCChain using Rao-Blackwell samples
function mean_rb_hmc(c::MCMCChain, par::Ranges=1:ncol(c.samples))
  nsamples, npars, nleaps = nrow(c.samples), length(par), length(c.diagnostics["leaps"][1])-1

  w = Array(Float64, nsamples, nleaps)
  sums = Array(Float64, nsamples, npars)

  for i = 1:nsamples
    for j = 1:nleaps
      w[i, j] = exp(c.diagnostics["leaps"][i][1].H-c.diagnostics["leaps"][i][j+1].H)
    end
  end

  for i = 1:nsamples
    for j = 1:npars
      s = c.samples[i, j]
      for k = 1:nleaps
        s += w[i, k]*c.diagnostics["leaps"][i][k+1].pars[j]
      end

      sums[i, j] = s/(nleaps+1)
    end
  end

  mean(sums, 1)
end

mean_rb_hmc(c::MCMCChain, par::Real) = mean_rb_hmc(c, par:par)

function mean_rb(c::MCMCChain, p::Union(Real, Ranges)=1:ncol(c.samples), s::Symbol=:hmc)
  if s == :hmc
    mean_rb_hmc(c, p)
  end
end


text:
function test_set_rxn()
    model = MetNets.toy_model()

    rxn = MetNets.Rxn("new_rxn";
        lb = -rand(), ub = rand(), c = rand())
    for met in model.mets
        push!(rxn.mets, met)
        push!(rxn.S, rand([0,-1, 1]))
    end

    rxni = 1
    MetNets.set_rxn!(model, rxni, rxn)

    # Test
    @test model.rxns[rxni] == rxn.id
    @test model.lb[rxni] == rxn.lb
    @test model.ub[rxni] == rxn.ub
    @test model.c[rxni] == rxn.c

    for (met, s) in zip(rxn.mets, rxn.S)
        meti = MetNets.metindex(model, met)
        @test model.S[meti, rxni] == s
    end
end
test_set_rxn()


text:
__precompile__()

## REFERENCE: https://discourse.julialang.org/t/using-as-a-wildcard-in-backtick-commands/6094/5

module ShellCommands

export @cmd_cmd, @ps_cmd, @busybox_cmd, @sh_cmd

macro cmd_cmd(s_str)
    s_expr = Meta.parse(string('"', escape_string(s_str), '"'))
    return esc(:(Cmd(Base.cmd_gen(("cmd", "/s", "/c", string('"', $s_expr, '"'))), windows_verbatim=true)))
end

macro ps_cmd(s_str)
    s_expr = Meta.parse(string('"', escape_string(s_str), '"'))
    return esc(:(Base.cmd_gen(("powershell", "-Command", $s_expr))))
end

macro busybox_cmd(s_str)
    s_expr = Meta.parse(string('"', escape_string(s_str), '"'))
    return esc(:(Base.cmd_gen(("busybox", "sh", "-c", $s_expr))))
end

macro sh_cmd(s_str)
    s_expr = Meta.parse(string('"', escape_string(s_str), '"'))
    return esc(:(Base.cmd_gen(("sh", "-c", $s_expr))))
end

end

max_stars_repo_name: jacob-roth/ShellCommands

max_issues_repo_name: jacob-roth/ShellCommands

max_forks_repo_name: jacob-roth/ShellCommands


text:
@testset "Built-in predicates" begin

# Test handling of and/N and or/N
clauses = @julog [
    shocked(P) <<= and(surprised(P), upset(P)),
    upset(P) <<= or(unhappy(P), angry(P)),
]

# Avery is surprised. Are they shocked? (No.)
facts = @julog [surprised(avery) <<= true]
@test resolve(@julog(shocked(avery)), [facts; clauses])[1] == false
# Bailey is surprised and angry. Are they shocked? (Yes.)
facts = @julog [surprised(bailey) <<= true, angry(bailey) <<= true]
@test resolve(@julog(shocked(bailey)), [facts; clauses])[1] == true
# Casey is unhappy. Is anyone upset? (Yes, Casey.)
facts = @julog [unhappy(casey) <<= true]
@test @varsub({P => casey}) in resolve(@julog(upset(P)), [facts; clauses])[2]

# Test handling of unifies/2 and not/1
clauses = @julog [
    child(zeus, kronos) <<= true,
    child(hera, kronos) <<= true,
    sibling(A, B) <<= child(A, C) & child(B, C) & not(unifies(A, B)),
]

# Is Hera her own sibling? (No)
@test resolve(@julog(sibling(hera, hera)), clauses)[1] == false
# Is Hera the sibling of Zeus? (Yes)
@test resolve(@julog(sibling(hera, zeus)), clauses)[1] == true

# Test the is/2 predicate (see also math_ops.jl for more tests)
clauses = @julog [
    square(X, Y) <<= is(Y, X * X),
    cube(X, Y) <<= is(Y, X * X * X)
]

# Is the cube of -1 the negation of its square? (Yes)
@test resolve(@julog([square(-1, S), cube(-1, C), S == -C]), clauses)[1] == true
# Is the cube of 2 twice its square? (Yes)
@test resolve(@julog([square(2, S), cube(2, C), C == 2*S]), clauses)[1] == true

# Test negation and double negation
clauses = @julog [red(roses) <<= true, blue(violets) <<= true]
@test resolve(@julog(not(red(roses))), clauses)[1] == false
@test resolve(@julog(not(blue(violets))), clauses)[1] == false
@test resolve(@julog(not(not(red(roses)))), clauses)[1] == true
@test resolve(@julog(not(not(blue(violets)))), clauses)[1] == true

# Test exists/2, forall/2, implies/2
clauses = @julog [
    human(pythagoras) <<= true,
    human(pyrrho) <<= true,
    human(zeno) <<= true,
    human(epicurus) <<= true,
    human(aristotle) <<= true,
    human(plato) <<= true,
    god(hera) <<= true,
    god(zeus) <<= true,
    god(aphrodite) <<= true,
    god(ares) <<= true,
    person(X) <<= human(X),
    person(X) <<= god(X),
    mortal(X) <<= human(X),
    immortal(X) <<= god(X)
]

# Is there a human who is mortal? (Yes.)
@test resolve(@julog(exists(human(X), mortal(X))), clauses)[1] == true
# Is there a god who is immortal? (Yes.)
@test resolve(@julog(exists(god(X), immortal(X))), clauses)[1] == true
# Is there a person who is mortal? (Yes.)
@test resolve(@julog(exists(person(X), mortal(X))), clauses)[1] == true

# Are all persons mortal? (No.)
@test resolve(@julog(forall(person(X), mortal(X))), clauses)[1] == false
# Are all humans mortal? (Yes.)
@test resolve(@julog(forall(human(X), mortal(X))), clauses)[1] == true
# Are all gods immortal? (Yes.)
@test resolve(@julog(forall(god(X), immortal(X))), clauses)[1] == true

# Is it true that if Hera is mortal, she's human? (Yes, since she's not mortal.)
@test resolve(@julog(imply(mortal(hera), human(hera))), clauses)[1] == true
# Is it true that if Hera is immortal, she's human? (No.)
@test resolve(@julog(imply(immortal(hera), human(hera))), clauses)[1] == false
# Is it true that for all people, being a god implies immortality? (Yes.)
@test resolve(@julog(forall(person(X), god(X) => immortal(X))), clauses)[1] == true
# Is it true that for all people, being a person implies immortality? (No.)
@test resolve(@julog(forall(person(X), person(X) => mortal(X))), clauses)[1] == false

# Of those who are persons, which are immortal?
sat, subst = resolve(@julog(person(X) => immortal(X)), clauses)
ans = Set([@varsub({X => hera}), @varsub({X => zeus}),
           @varsub({X => aphrodite}), @varsub({X => ares})])
@test ans == Set(subst)

# Test findall/3 and countall/2
humans = @julog list[pythagoras, pyrrho, zeno, epicurus, aristotle, plato]
@test resolve(@julog(findall(X, human(X), :humans)), clauses)[1] == true
@test resolve(@julog(countall(human(X), 6)), clauses)[1] == true
gods = @julog list[hera, zeus, aphrodite, ares]
@test resolve(@julog(findall(X, god(X), :gods)), clauses)[1] == true
@test resolve(@julog(countall(god(X), 4)), clauses)[1] == true
persons = @julog list[pythagoras, pyrrho, zeno, epicurus, aristotle, plato,
                      hera, zeus, aphrodite, ares]
@test resolve(@julog(findall(X, person(X), :persons)), clauses)[1] == true
@test resolve(@julog(countall(person(X), 10)), clauses)[1] == true

# Test cut and fail by preventing infinite loops
clauses = @julog [
    fakeloop1(A) <<= fakeloop1(B),
    fakeloop1(A) <<= cut,
    fakeloop2(A) <<= fail & fakeloop2(B),
]

@test resolve(@julog(fakeloop1(0)), clauses)[1] == true
@test resolve(@julog(fakeloop2(0)), clauses)[1] == false

# Test the meta-call predicate call/N
clauses = @julog [
    test(x, y) <<= true,
    pred(test) <<= true,
    metatest1(A, B, C) <<= call(A, B, C),
    metatest2(A, B, C) <<= pred(A) & call(A, B, C)
]

@test_throws Exception resolve(@julog(call(P, A, B)), clauses)
@test @varsub({A => x, B => y}) in resolve(@julog(call(test, A, B)), clauses)[2]
@test @varsub({A => x, B => y}) in resolve(@julog(call(test(A), B)), clauses)[2]
@test @varsub({B => y}) in resolve(@julog(call(test(x), B)), clauses)[2]
@test @varsub({A => x}) in resolve(@julog(call(test(A), y)), clauses)[2]

@test_throws Exception resolve(@julog(metatest1(P, x, y)), clauses)
@test @varsub({B => y}) in resolve(@julog(metatest1(test, x, B)), clauses)[2]
@test @varsub({P => test}) in resolve(@julog(metatest2(P, x, y)), clauses)[2]
@test @varsub({B => y}) in resolve(@julog(metatest2(test, x, B)), clauses)[2]

end


text:
function move_order_to_invest(orderId, client, user_id, symbol)
    order = client.get_order(symbol=symbol, orderId=orderId)
    trades = client.get_my_trades(symbol=order["symbol"], orderId=order["orderId"])

    if order["price"] == "0.00000000"
        order["price"] = string(parse(Float64, order["cummulativeQuoteQty"]) / parse(Float64, order["executedQty"]))
    end
    farm_coin = split(order["symbol"], "EUR")[1]
    base_coin = "EUR"
    invest = Dict(
        "amount" => order["origQty"],
        "farm_coin" => farm_coin,
        "base_coin" => base_coin,
        "mode" => "now",
        "t" => string(unix2datetime(order["time"]/1000)),
        "symbol" => order["symbol"],
        "value" => order["cummulativeQuoteQty"],
        "order" => order,
        "trades" => trades,
        "status" => "Done"
    )

    y, m = split(split(string(invest["t"]), "T")[1], "-")[1:2]
    file = string(joinpath(Pinguu.base_dir[2:end], "invest/users/$user_id/filled/", y, "/", m, "/", invest["symbol"], "_", invest["t"], ".json"))
    !isdir(dirname(file)) ? mkpath(dirname(file)) : ""
    open(file, "w") do f
        JSON.print(f, invest, 4)
    end
    return invest
end
export move_order_to_invest


text:
using Plots
plotly()

using simulator.fresnelltools
using simulator.materials.spesifics
using simulator.matrixcore

θs = 0:1e-3:π/2
θ = 0.806
λs = 350e-9:1e-9:1800e-9
λ = 1434e-9

d_ags = 10e-9:1e-9:200e-9
d_ag = 50e-9
d_sis = 10e-9:1e-9:200e-9
d_si = 20e-9
d_pds = 10e-9:1e-9:200e-9
d_pd = 5e-9

n_fiber = SiO2_core_Sellmeier
n_ag = Ag
n_silicafilm = SiO2_thinfilm_Ciprian
n_palladium(λ) = Pd_Werner(λ) - 0im
n_air = 1

function structure(λ, θ, d_ag, d_si, d_pd)
    interface_fiber_ag_te, interface_fiber_ag_tm, θag = FresnellBoundrary(n_fiber(λ), n_ag(λ), θ)
    bulk_ag = FresnellSlab(n_ag(λ), 2*π/λ, d_ag, θag)
    interface_ag_si_te, interface_ag_si_tm, θsi = FresnellBoundrary(n_ag(λ), n_silicafilm(λ), θag)
    bulk_si = FresnellSlab(n_silicafilm(λ), 2*π/λ, d_si, θsi)
    interface_si_pd_te, interface_si_pd_tm, θpd = FresnellBoundrary(n_silicafilm(λ), n_palladium(λ), θsi)
    bulk_pd = FresnellSlab(n_palladium(λ), 2*π/λ, d_pd, θpd)
    interface_pd_air_te, interface_pd_air_tm, _ = FresnellBoundrary(n_palladium(λ), n_air, θpd)

    te_system = CascadeScattering([interface_fiber_ag_te, bulk_ag, interface_ag_si_te, bulk_si, interface_si_pd_te, bulk_pd, interface_pd_air_te])
    tm_system = CascadeScattering([interface_fiber_ag_tm, bulk_ag, interface_ag_si_tm, bulk_si, interface_si_pd_tm, bulk_pd, interface_pd_air_tm])

    te_system, tm_system
end

function structure_inside(λ, θ, d_ag, d_si, d_pd)
    bulk_fiber(d) = FresnellSlab(n_fiber(λ), 2*π/λ, d, θ)
    interface_fiber_ag_te, interface_fiber_ag_tm, θ_ag = FresnellBoundrary(n_fiber(λ), n_ag(λ), θ)
    interface_fiber_ag_tm_f(d) = interface_fiber_ag_tm
    interface_fiber_ag_te_f(d) = interface_fiber_ag_te
    bulk_ag(d) = FresnellSlab(n_ag(λ), 2*π/λ, d, θ_ag)
    interface_ag_si_te, interface_ag_si_tm, θ_si = FresnellBoundrary(n_ag(λ), n_silicafilm(λ), θ_ag)
    interface_ag_si_tm_f(d) = interface_ag_si_tm
    interface_ag_si_te_f(d) = interface_ag_si_te
    bulk_si(d) = FresnellSlab(n_silicafilm(λ), 2*π/λ, d, θ_si)

    # interface_si_pd_te, interface_si_pd_tm, θ_pd = FresnellBoundrary(n_silicafilm(λ), n_palladium(λ), θ_si)
    # interface_si_pd_te_f(d) = interface_si_pd_te
    # interface_si_pd_tm_f(d) = interface_si_pd_tm
    # bulk_pd(d) = FresnellSlab(n_palladium(λ), 2*π/λ, d, θ_pd)
    # interface_pd_air_te, interface_pd_air_tm, θ_air = FresnellBoundrary(n_palladium(λ), n_air, θ_air)
    # interface_pd_air_te_f(d) = interface_pd_air_te
    # interface_pd_air_tm_f(d) = interface_pd_air_tm
    # bulk_air(d) = FresnellSlab(n_air, 2*π/λ, d, θ_air)

    interface_si_air_te, interface_si_air_tm, θ_air = FresnellBoundrary(n_silicafilm(λ), n_air, θ_si)
    interface_si_air_te_f(d) = interface_si_air_te
    interface_si_air_tm_f(d) = interface_si_air_tm
    bulk_air(d) = FresnellSlab(n_air, 2*π/λ, d, θ_air)

    interrogated = (Up, Un) -> abs(Up + Un)^2

    u_values_tm, d_values_tm = Interrogator(
        [bulk_fiber, interface_fiber_ag_tm_f, bulk_ag, interface_ag_si_tm_f, bulk_si, interface_si_air_tm_f, bulk_air],
        [10e-9, 0, d_ag, 0, d_si, 0, 10e-9], 1e-10,
        interrogated)
    u_values_te, d_values_te = Interrogator(
        [bulk_fiber, interface_fiber_ag_te_f, bulk_ag, interface_ag_si_te_f, bulk_si, interface_si_air_te_f, bulk_air],
        [10e-9, 0, d_ag, 0, d_si, 0, 10e-9], 1e-10,
        interrogated)

    u_values_te, d_values_te, u_values_tm, d_values_tm
end

R_te(λ, θ, d_ag, d_si, d_pd) = abs((structure(λ, θ, d_ag, d_si, d_pd)[1] * [1, 0])[2])^2
R_tm(λ, θ, d_ag, d_si, d_pd) = abs((structure(λ, θ, d_ag, d_si, d_pd)[2] * [1, 0])[2])^2
R_ratio(λ, θ, d_ag, d_si, d_pd) = R_te(λ, θ, d_ag, d_si, d_pd) / R_tm(λ, θ, d_ag, d_si, d_pd)

interrogated_u_te, interrogated_d_te, interrogated_u_tm, interrogated_d_tm = structure_inside(λ, θ, d_ag, d_si, d_pd)

# plot(λs, R_ratio.(λs, θ, d_ag, d_si, d_pd), label="θ=0.856", title="TE Reflection / TM Reflection at d=53nm")
# plot(θs, (θ) -> (R_ratio.(λ, θ + 1e-9im, d_ag, d_si, d_pd)))
# plot(d_ags, R_ratio.(λ, θ, d_ags, d_si, d_pd))
# plot(d_sis, R_ratio.(λ, θ, d_ag, d_sis, d_pd))
# plot(d_pds, R_ratio.(λ, θ, d_ag, d_si, d_pds))

# plot(λs, R_te.(λs, θ, d_ag, d_si, d_pd), label="te, Pd-structure, Different Pd", title="Reflection at silica-palladium-air", legend=:outertopright)
# plot!(λs, R_tm.(λs, θ, d_ag, d_si, d_pd), label="tm, Pd-structure, Different Pd", xlabel="wavelength", ylabel="reflected power")

# plot(interrogated_d_te .* 1e9, interrogated_u_te, label="te, λ=$λ", title="Sum of fields at plasmone conditions", legend=:topright, xlabel="distance in structure", ylabel="total power")
# plot!(interrogated_d_tm .* 1e9, interrogated_u_tm, label="tm λ=$λ")


text:
mutable struct xNES_state{T}
    ημ::T
    ησ::T
    samples::Int
    termination::Termination{T}
    function xNES_state{T}(d::Integer;P...) where T
        ημ=T(1)
        ησ=T( (9+3*log(d))/(5*d*√d) )
        samples=4 + ceil(Int, log(3*d))
        Term=Termination{T}()
        S=new{T}(ημ,ησ,samples,Term)
        for k in keys(P)
            hasfield(sNES_state{T},k) && setfield!(S,k,P[k])
            hasfield(Termination{T},k) && setfield!(Term,k,P[k])
        end
        S
    end
end

struct xNES <: OptMethod end

function δ(T,i,j)
    ifelse(i==j,one(T),zero(T))
end

function exponential_nes(f,μ0::AbstractVector{T},A::AbstractMatrix{T},params::xNES_state{T}) where T
    n=params.samples
    ημ=params.ημ
    ηB=ησ=params.ησ
    μ=copy(μ0)
    d=length(μ)
    Z=Array{T}(undef,d,n)
    σ=abs(det(A))^(1/d)
    F=fill(f(μ),n)
    B=A./σ
    idx=collect(1:n)
    u=utility_function{T}(n)
    Gδ=zeros(T,d)
    GM=zeros(T,d,d)
    tmp_μ=copy(μ)
    iterations=0
    @inbounds while !termination_criteria(params.termination,σ,iterations)
        randn!(Z)
        for i in 1:n
            #F[i] = f(μ .+ σ .* B * Z[:,i])
            mul!(tmp_μ,B,view(Z,:,i))
            tmp_μ .= μ .+ tmp_μ .* σ
            F[i] = f(tmp_μ)
        end
        sortperm!(idx,F)
        #Gδ=sum(u[i] .* Z[:,idx[i]] for i in 1:n)
        #GM=sum(u[i] .* (Z[:,idx[i]] * (Z[:,idx[i]]') - I) for i in 1:n)
        let i=1
            j = idx[i]
            for k2 in 1:d
                Gδ[k2] = u[i] * Z[k2,j]
                for k1 in 1:d
                    GM[k1,k2] = u[i] * (Z[k1,j] * Z[k2,j] - δ(T,k1,k2))
                end
            end
        end
        for i in 2:n
            j = idx[i]
            for k2 in 1:d
                Gδ[k2] += u[i] * Z[k2,j]
                for k1 in 1:d
                    GM[k1,k2] += u[i] * (Z[k1,j] * Z[k2,j] - δ(T,k1,k2))
                end
            end
        end
        #=
        Gσ = tr(GM) / d
        GB=GM - Gσ * I
        μ=μ .+ ημ * σ *B *Gδ
        σ=σ * exp(ησ/2 * Gσ)
        B=B * exp(ηB/2 .* GB)
        =#
        Gσ = tr(GM) / d
        for i in 1:d
            GM[i,i]-=Gσ
        end
        GM .*= ηB/2
        Gσ *= ησ/2
        mul!(tmp_μ,B,Gδ)
        μ .+= ημ .* σ .* tmp_μ
        σ *= exp(Gσ)
        mul!(GM,B,exp(GM))
        B.=GM
        iterations+=1
    end
    return (sol=μ, cost=f(μ))
end

function exponential_nes(f,μ::AbstractVector{T},σ::AbstractVector{T},params::xNES_state{T}) where T
    A=diagm(σ)
    exponential_nes(f,μ,A,params)
end

function exponential_nes(f,μ::AbstractVector{T},σ::T,params::xNES_state{T}) where T
    exponential_nes(f,μ,fill(σ,length(μ)),params)
end

function exponential_nes(f,μ::AbstractVector{T},σ; P...) where T
    exponential_nes(f,μ,σ,xNES_state{T}(length(μ);P...))
end

pickmethod(::Type{xNES}) = exponential_nes


text:
macro CSI_str(str)
    return :(string("\x1b[", $(esc(str)), "m"))
end
const TYPE_COLOR = CSI"36"
# const RED_COLOR = CSI"35"
# const BLUE_COLOR = CSI"34;1"
# const YELLOW_COLOR = CSI"33"
const ORANGE_COLOR = CSI"38;5;208"
const NO_COLOR = CSI"0"

import Base: summary
Base.summary(DS::AbstractDataSet) = string(TYPE_COLOR, nameof(typeof(DS)), NO_COLOR, " with N=$(Npoints(DS)), xdim=$(xdim(DS)) and ydim=$(ydim(DS))")

GeneratedFromAutoDiff(F::Function) = occursin("Autodmodel", string(nameof(typeof(F))))
GeneratedFromAutoDiff(F::ModelMap) = GeneratedFromAutoDiff(F.Map)
GeneratedFromSymbolic(F::Function) = occursin("SymbolicModel", string(nameof(typeof(F))))
GeneratedFromSymbolic(F::ModelMap) = GeneratedFromSymbolic(F.Map)

###### Useful info: Autodmodel? Symbolic? StaticArray output? In-place?
function Base.summary(DM::AbstractDataModel)
    # Also use "RuntimeGeneratedFunction" string from build_function in ModelingToolkit.jl
    string(TYPE_COLOR, nameof(typeof(DM)),
    NO_COLOR, " containing ",
    TYPE_COLOR, nameof(typeof(Data(DM))),
    NO_COLOR, ". Model jacobian: ",
    ORANGE_COLOR, (GeneratedFromAutoDiff(dPredictor(DM)) ? "AutoDiff" : (GeneratedFromSymbolic(dPredictor(DM)) ? "Symbolic" : "manually provided")),
    NO_COLOR)
end

# http://docs.junolab.org/stable/man/info_developer/#
# hastreeview, numberofnodes, treelabel, treenode
TreeViews.hastreeview(x::Union{AbstractDataSet,AbstractDataModel,ModelMap}) = true
TreeViews.numberofnodes(x::AbstractDataSet) = 4
TreeViews.numberofnodes(x::AbstractDataModel) = 4
TreeViews.numberofnodes(x::ModelMap) = 4
TreeViews.numberofnodes(x::CompositeDataSet) = 1
TreeViews.numberofnodes(x::GeneralizedDataSet) = 1

function TreeViews.treelabel(io::IO, DS::Union{AbstractDataSet,AbstractDataModel}, mime::MIME"text/plain" = MIME"text/plain"())
    show(io, mime, Text(Base.summary(DS)))
end
# To hide the treenode display, simply return missing:
# treenode(x::Foo, i::Int) = missing

import Base: show
#### Need proper show() methods for DataSet, DataModel, ModelMap
#### Show Distribution types for DataSetExact
function Base.show(io::IO, mime::MIME"text/plain", DS::AbstractDataSet)
    println(io, "$(nameof(typeof(DS))) with N=$(Npoints(DS)), xdim=$(xdim(DS)) and ydim=$(ydim(DS)):")
    print(io, "x-data: ");    show(io, mime, xdata(DS));    print(io, "\n")
    if DS isa DataSetExact
        if typeof(xsigma(DS)) <: AbstractVector
            println(io, "Standard deviation associated with x-data:")
            show(io, mime, xsigma(DS))
        else
            println(io, "Covariance matrix associated with x-data:")
            show(io, mime, xsigma(DS))
        end
        print(io, "\n")
    end
    print(io, "y-data: ");    show(io, mime, ydata(DS));    print(io, "\n")
    if typeof(ysigma(DS)) <: AbstractVector
        println(io, "Standard deviation associated with y-data:")
        show(io, mime, ysigma(DS))
    else
        println(io, "Covariance matrix associated with y-data:")
        show(io, mime, ysigma(DS))
    end
end

function Base.show(io::IO, DS::AbstractDataSet)
    println(io, "$(nameof(typeof(DS))) with N=$(Npoints(DS)), xdim=$(xdim(DS)) and ydim=$(ydim(DS)):")
    print(io, "x-data: ");    show(io, xdata(DS));    print(io, "\n")
    if DS isa DataSetExact
        if typeof(xsigma(DS)) <: AbstractVector
            println(io, "Standard deviation associated with x-data:")
            show(io, xsigma(DS))
        else
            println(io, "Covariance matrix associated with x-data:")
            show(io, xsigma(DS))
        end
        print(io, "\n")
    end
    print(io, "y-data: ");    show(io, ydata(DS));    print(io, "\n")
    if typeof(ysigma(DS)) <: AbstractVector
        println(io, "Standard deviation associated with y-data:")
        show(io, ysigma(DS))
    else
        println(io, "Covariance matrix associated with y-data:")
        show(io, ysigma(DS))
    end
end

function Base.show(io::IO, mime::MIME"text/plain", GDS::GeneralizedDataSet)
    println(io, "$(nameof(typeof(GDS))) with N=$(Npoints(GDS)), xdim=$(xdim(GDS)) and ydim=$(ydim(GDS)):")
    print(io, "Combined x-y data: ");    show(io, mime, GetMean(dist(GDS)));    print(io, "\n")
    print(io, "Combined x-y covariance: ");    show(io, mime, Sigma(dist(GDS)));    print(io, "\n")
end
function Base.show(io::IO, GDS::GeneralizedDataSet)
    println(io, "$(nameof(typeof(GDS))) with N=$(Npoints(GDS)), xdim=$(xdim(GDS)) and ydim=$(ydim(GDS)):")
    print(io, "Combined x-y data: ");    show(io, GetMean(dist(GDS)));    print(io, "\n")
    print(io, "Combined x-y covariance: ");    show(io, Sigma(dist(GDS)));    print(io, "\n")
end

function Base.show(io::IO, DM::AbstractDataModel)
    Expr = SymbolicModel(DM)
    IsLin = try IsLinearParameter(DM) catch; nothing end
    println(io, "$(nameof(typeof(DM))) containing a $(nameof(typeof(Data(DM))))")
    Jac = if GeneratedFromAutoDiff(dPredictor(DM))
        "generated via automatic differentiation"
    elseif GeneratedFromSymbolic(dPredictor(DM))
        "symbolically provided"
    else
        "manually provided by user"
    end
    println(io, "Model jacobian " * Jac)
    if DM isa DataModel
        println(io, "Maximum Likelihood Estimate: $(MLE(DM))")
        println(io, "Maximal value of log-likelihood: $(LogLikeMLE(DM))")
    end
    Expr[1] == 'y' && println(io, "Model Expr:  $Expr")
    !isnothing(IsLin) && println(io, "Model parametrization linear in n-th parameter: $(IsLin)")
end


text:
#generate amino acid alignments
function simAlignments(numFiles,numSeqs,parentdir)
	infiles = Vector{String}(undef, numFiles)
	for f in 1:numFiles
		file, io = mktemp(parentdir)
    lenSeq = rand(100:1000)
		infiles[f] = file
    w = FASTX.FASTA.Writer(io)
		for n in 1:numSeqs
			seq = BioSequences.randaaseq(lenSeq) # all seqs need to have same length
			rec = FASTX.FASTA.Record(string(n),seq)
			write(w, rec)
		end
    close(io)
	end
	return(infiles)
end

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
"""
    ChannelData

Stores the results from a single channel.
"""
struct ChannelData
    channelname::AbstractString
    molecules::Vector{<:DataEntity}
    neighbormolecules::Vector{<:DataEntity}
end

"""
    ChannelDataOptimizing

Stores the results from many molecule lists. Intended for parameter optimizing only.
"""
struct ChannelDataOptimizing
    channelname::AbstractString
    molecules::Vector{Vector{<:DataEntity}}
    neighbormolecules::Vector{<:DataEntity}
end


text:
################################################################################
# Displacements
################################################################################
function displacement(joint::Rotational,
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        vmat=true)

    q = inv(joint.orientation_offset) * inv(qa) * qb
    vmat ? (return Vmat(q)) : (return q)
end

function displacement_jacobian_configuration_unstable(relative::Symbol, joint::Rotational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        attjac::Bool=true, vmat::Bool=true) where T
    X = szeros(T, 3, 3)
    if relative == :parent
		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
		attjac && (Q *= LVᵀmat(qa))
    elseif relative == :child
		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
		attjac && (Q *= LVᵀmat(qb))
	end
	vmat && (Q = Vmat() * Q)
	return X, Q
end

function displacement_jacobian_configuration(relative::Symbol, joint::Rotational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
		) where T
    X = szeros(T, 3, 3)
    if relative == :parent
		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
    elseif relative == :child
		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
	end
	return X, Vmat() * Q
end

# function displacement_jacobian_configuration(relative::Symbol, joint::Rotational{T},
#         xa::AbstractVector, qa::Quaternion,
#         xb::AbstractVector, qb::Quaternion;
#         attjac::Bool=true, vmat=true) where T
#     X = szeros(T, 3, 3)
#     if relative == :parent
# 		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
# 		attjac && (Q *= LVᵀmat(qa))
#     elseif relative == :child
# 		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
# 		attjac && (Q *= LVᵀmat(qb))
# 	end
# 	vmat && (Q = Vmat() * Q)
# 	return X, Q
# end

################################################################################
# Coordinates
################################################################################
function minimal_coordinates(joint::Rotational,
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion)

    return nullspace_mask(joint) * rotation_vector(displacement(joint, xa, qa, xb, qb, vmat=false))
end

function minimal_coordinates_jacobian_configuration(relative::Symbol, joint::Rotational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        attjac::Bool=true) where T

    A = nullspace_mask(joint)
    q = displacement(joint, xa, qa, xb, qb, vmat=false)
    X, Q = displacement_jacobian_configuration_unstable(relative, joint, xa, qa, xb, qb, attjac=attjac, vmat=false)
    ∂rv∂q = drotation_vectordq(q)

	return A * [X ∂rv∂q * Q]
end

################################################################################
# Set Coordinates
################################################################################
function set_minimal_coordinates!(joint::Rotational,
	pnode::Node, cnode::Node,
	timestep;
	Δθ::AbstractVector=szeros(input_dimension(joint)))

	orientation_offset = joint.orientation_offset
	qa = pnode.state.q2
	Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
	Δq = axis_angle_to_quaternion(Aᵀ * Δθ)
	qb = qa * orientation_offset * Δq
	set_maximal_configurations!(cnode; x=cnode.state.x2, q = qb)
	set_previous_configuration!(cnode, timestep)
	return nothing
end

################################################################################
# Velocities
################################################################################
function minimal_velocities(joint::Rotational,
		xa::AbstractVector, va::AbstractVector, qa::Quaternion, ωa::AbstractVector,
		xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ωb::AbstractVector,
		timestep)

	orientation_offset = joint.orientation_offset
	A = nullspace_mask(joint)

	# 1 step backward in time
	qa1 = next_orientation(qa, -ωa, timestep)
	qb1 = next_orientation(qb, -ωb, timestep)

	q = inv(orientation_offset) * inv(qa) * qb
	q1 = inv(orientation_offset) * inv(qa1) * qb1
	return A * rotation_vector(inv(q1) * q) ./ timestep
end

function minimal_velocities_jacobian_configuration(relative::Symbol, joint::Rotational{T},
	xa::AbstractVector, va::AbstractVector, qa::Quaternion, ωa::AbstractVector,
	xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ωb::AbstractVector,
	timestep) where T

	orientation_offset = joint.orientation_offset
	A = nullspace_mask(joint)
	nu = input_dimension(joint)

	# 1 step backward in time
	qa1 = next_orientation(qa, -ωa, timestep)
	qb1 = next_orientation(qb, -ωb, timestep)

	q = inv(orientation_offset) * inv(qa) * qb
	q1 = inv(orientation_offset) * inv(qa1) * qb1

	X = szeros(T, nu, 3)

	if relative == :parent
		Q = 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Rmat(q) * Tmat() * Rmat(qb1) * Lmat(inv(orientation_offset)) * Tmat() * rotational_integrator_jacobian_orientation(qa, -ωa, timestep, attjac=false)
		Q += 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Lmat(inv(q1)) * Rmat(qb) * Lmat(inv(orientation_offset)) * Tmat()
		Q *= LVᵀmat(qa)
	elseif relative == :child
		Q = 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Rmat(q) * Tmat() * Lmat(inv(orientation_offset) * inv(qa1)) * rotational_integrator_jacobian_orientation(qb, -ωb, timestep, attjac=false)
		Q += 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Lmat(inv(q1) * inv(orientation_offset) * inv(qa))
		Q *= LVᵀmat(qb)
	end

	return [X Q]
end

function minimal_velocities_jacobian_velocity(relative::Symbol, joint::Rotational{T},
	xa::AbstractVector, va::AbstractVector, qa::Quaternion, ωa::AbstractVector,
	xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ωb::AbstractVector,
	timestep) where T

	orientation_offset = joint.orientation_offset
	A = nullspace_mask(joint)
	nu = input_dimension(joint)

	# 1 step backward in time
	qa1 = next_orientation(qa, -ωa, timestep)
	qb1 = next_orientation(qb, -ωb, timestep)

	q = inv(orientation_offset) * inv(qa) * qb
	q1 = inv(orientation_offset) * inv(qa1) * qb1

	V = szeros(T, nu, 3)
	if relative == :parent
		Ω = 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Rmat(q) * Tmat() * Lmat(inv(orientation_offset)) * Rmat(qb1) * Tmat() * -rotational_integrator_jacobian_velocity(qa, -ωa, timestep)
	elseif relative == :child
		Ω = 1.0 / timestep * A * drotation_vectordq(inv(q1) * q) * Rmat(q) * Tmat() * Lmat(inv(orientation_offset) * inv(qa1)) * -rotational_integrator_jacobian_velocity(qb, -ωb, timestep)
	end
	return [V Ω]
end


text:
"""
    get_pMag(T::Real, MASS::Real)

Description of get_pMag
--------------------------
Computes magnitude of momentum as:
``` |p| = √((MASS + T)² - MASS²) ```.

Input arguments are:
+ T    : kinetic energy in [MeV]
+ MASS : particle mass in [MeV/c²]
"""
function get_pMag(T::Real, MASS::Real)
	return sqrt( (MASS+T)^2 - MASS^2)
end

"""
    get_px(pMag, ϕ, θ)

Description of get_px
--------------------------
Computes x-component of the momentum as:
``` pₓ = |p|*cos(ϕ)*sin(θ) ```

Input arguments are:
+ pMag : momentum magnitude in [MeV]
+ ϕ    : azimuth angle in [rad]
+ θ    : polar angle in [rad]
"""
function get_px(pMag, ϕ, θ)
    return pMag*cos(ϕ)*sin(θ)
end

"""
    get_py(pMag, ϕ, θ)

Description of get_py
--------------------------
Computes y-component of the momentum as:
``` p𝔂 = |p|*sin(ϕ)*sin(θ) ```

Input arguments are:
+ pMag : momentum magnitude in [MeV]
+ ϕ    : azimuth angle in [rad]
+ θ    : polar angle in [rad]
"""
function get_py(pMag, ϕ, θ)
    return pMag*sin(ϕ)*sin(θ)
end

"""
    ```get_pz(pMag, θ)```

Description of get_pz
--------------------------
Computes z-component of the momentum as:
`` p𝐳 = |p|*cos(θ) ``

Input arguments are:
+ pMag : momentum magnitude in [MeV]
+ θ    : polar angle in [rad]
"""
function get_pz(pMag, θ)
    return pMag*cos(θ)
end

"""
```get_first_vector(T::Real, MASS::Real)```

Description of ```get_first_vector```
-------------------------------
Returns the momentum vector of the first electron. Momentum magnitude is calculated from the electron kinetic energy (T).
The direction is generated uniformly in a sphere.

The momentum components are (in order):
    + pMag - Magnitude of the momentum vector (given by the electron Energy)
    + px   - x-component
    + py   - y-component
    + pz   - z-component
    + θ    - polar angle  (used in calculation of the second momentum vector)
    + ϕ    - azimuthal angle (used in calculation of the second momentum vector)
"""
function get_first_vector(T::Real, MASS::Real)
    pMag = get_pMag(T, MASS)
    ϕ = sample_phi()
    θ = sample_theta()
    px = get_px(pMag, ϕ, θ)
    py = get_py(pMag, ϕ, θ)
    pz = get_pz(pMag, θ)
    return [pMag, px, py, pz, θ, ϕ]
end

"""
```get_pPrime(T::Real, MASS::Real)```

Description of get_pPrime
-------------------------------
Returns the momentum (prime) vector of the second electron. The vector is generated as a random point from a circle projected onto the x-y plane.
The circle radius is given by ``r = \\sin(θdif)|p1|``, where ``\\theta_{dif}`` is the generated angle between the two electrons and ``|p1|``
is the magnitude of the momentum vector of the first electron. The process assumes p1 to be directed along the z-axis (hence projection onto x-y plane).

The vector components are (in order):
    + ``\\sin(\\theta_{dif})\\cos(\\phi_2)pMag`` - x-component
    + ``\\sin(\\theta_{dif})\\sin(\\phi_2)pMag`` - y-component
    + ``\\cos(\\theta_{dif})pMag``               - z-component

```get_pPrime(T::Real, MASS::Real, θdif)```

Generates momentum (prime) vector of the second electron, where θdif is predefined.

"""
function get_pPrime(T::Real, MASS::Real)
    pMag = get_pMag(T, MASS)
    θdif = sample_theta_dif()
    ϕ2   = sample_phi()

    return [ sin(θdif)*cos(ϕ2)*pMag ; sin(θdif)*sin(ϕ2)*pMag ; cos(θdif)*pMag ] # returns column vector of p2'
end

function get_pPrime(T::Real, MASS::Real, θdif)
    pMag = get_pMag(T, MASS)
    ϕ2   = sample_phi()

    return [ sin(θdif)*cos(ϕ2)*pMag ; sin(θdif)*sin(ϕ2)*pMag ; cos(θdif)*pMag ] # returns column vector of p2'
end

"""
```get_second_vector(T::Real, MASS::Real, p1::DataFrameRow) ```   p1 given as dataframe row
```get_second_vector(T::Real, MASS::Real, p1::Vector{<:Real})``` p1 given as vector

Description of get_second_vector
--------------------------------
Returns the momentum vector of the second electron. Momentum magnitude is calculated from the electron kinetic energy (T).
The direction is generated from sampling θdif angle (see ```sample_theta_dif()```). In order to calculate the direction: first a p2Prime vector is generated
(see ```get_pPrime(T::Real, MASS::Real)```) by sampling uniformly from a circle projected onto x-y plane,
where the circle is determined by the p1 vector directed along z-axis. Afterwards, p2Prime is rotated along 1. y-axis, 2. z-axis to obtain p2 in the direction
of a cone around p1, where the cone is given by θdif. The equation for p2 is:

`` p2 = Rz*Ry*p2Prime ``

The momentum components are (in order):
    + pMag - Magnitude of the momentum vector (given by the electron Energy)
    + px   - x-component
    + py   - y-component
    + pz   - z-component
    + θ    - polar angle  (used in calculation of the second momentum vector)
    + ϕ    - azimuthal angle (used in calculation of the second momentum vector)
"""
function get_second_vector(T::Real, MASS::Real, p1::DataFrameRow)
    p2Prime = get_pPrime(T, MASS) # vector generated by randomly chosing point on circle projected onto xy plane

    p2 = [
        cos(p1.phi)*( p2Prime[1]*cos(p1.theta ) + p2Prime[3]*sin(p1.theta) ) - p2Prime[2]*sin(p1.phi)
        sin(p1.phi)*( p2Prime[1]*cos(p1.theta ) + p2Prime[3]*sin(p1.theta) ) + p2Prime[2]*cos(p1.phi)
        p2Prime[3]*cos(p1.theta ) - p2Prime[1]*sin(p1.theta )
    ]

    p2Mag = get_pMag(T, MASS)
    θ2 = acos(p2[3]/p2Mag)
    ϕ2 = atan(p2[2],p2[1])   # computes atan(y/x) while taking care of undefined

    while (p2[1] == NaN)
        get_second_vector(T, MASS, p1)
    end

    return [ p2Mag, p2[1], p2[2], p2[3], θ2, ϕ2 ]
end

function get_second_vector(T::Real, MASS::Real, p1::Vector{<:Real})
    p2Prime = get_pPrime(T, MASS) # vector generated by randomly chosing point on circle projected onto xy plane

    p2 = [
        cos(p1[6] )*( p2Prime[1]*cos(p1[5] ) + p2Prime[3]*sin(p1[5]) ) - p2Prime[2]*sin(p1[6])
        sin(p1[6] )*( p2Prime[1]*cos(p1[5] ) + p2Prime[3]*sin(p1[5]) ) + p2Prime[2]*cos(p1[6])
        p2Prime[3]cos(p1[5] ) - p2Prime[1]*sin(p1[5] )
    ]

    p2Mag = get_pMag(T, MASS)
    θ2 = acos(p2[3]/p2Mag)
    ϕ2 = atan(p2[2],p2[1])   # computes atan(y/x) while taking care of undefined

    while (p2[1] == NaN)
        get_second_vector(T, MASS, p1)
    end

    return [ p2Mag, p2[1], p2[2], p2[3], θ2, ϕ2 ]
end

function get_second_vector(T::Real, MASS::Real, p1::DataFrameRow, θdif::Real)
    p2Prime = get_pPrime(T, MASS, θdif) # vector generated by randomly chosing point on circle projected onto xy plane

    p2 = [
        cospi(p1.phi /π)*( p2Prime[1]*cospi(p1.theta /π) + p2Prime[3]*sinpi(p1.theta/π) ) - p2Prime[2]*sinpi(p1.phi/π)
        sinpi(p1.phi /π)*( p2Prime[1]*cospi(p1.theta /π) + p2Prime[3]*sinpi(p1.theta/π) ) + p2Prime[2]*cospi(p1.phi/π)
        p2Prime[3]cospi(p1.theta /π) - p2Prime[1]*sinpi(p1.theta /π)
    ]

    p2Mag = get_pMag(T, MASS)
    θ2 = acos(p2[3]/p2Mag)
    ϕ2 = atan(p2[2],p2[1])   # computes atan(y/x) while taking care of undefined

    while (p2[1] == NaN)
        get_second_vector(T, MASS, p1, θdif)
    end

    return [ p2Mag, p2[1], p2[2], p2[3], θ2, ϕ2 ]
end

function get_second_vector(T::Real, MASS::Real, p1::Vector{<:Real},  θdif::Real)
    p2Prime = get_pPrime(T, MASS, θdif) # vector generated by randomly chosing point on circle projected onto xy plane

    p2 = [
        cospi(p1[6] /π)*( p2Prime[1]*cospi(p1[5] /π) + p2Prime[3]*sinpi(p1[5]/π) ) - p2Prime[2]*sinpi(p1[6]/π)
        sinpi(p1[6] /π)*( p2Prime[1]*cospi(p1[5] /π) + p2Prime[3]*sinpi(p1[5]/π) ) + p2Prime[2]*cospi(p1[6]/π)
        p2Prime[3]cospi(p1[5] /π) - p2Prime[1]*sinpi(p1[5] /π)
    ]

    p2Mag = get_pMag(T, MASS)
    θ2 = acos(p2[3]/p2Mag)
    ϕ2 = atan(p2[2],p2[1])   # computes atan(y/x) while taking care of undefined

    while (p2[1] == NaN)
        get_second_vector(T, MASS, p1, θdif)
    end

    return [ p2Mag, p2[1], p2[2], p2[3], θ2, ϕ2 ]
end


text:
# (31+42)::Float64 #> TypeError(:typeassert,"",Float64,73)
# convert(Int64, 7.01) #> ERROR: InexactError()
convert(Int64, 7.0) #> 7
Int64(7.0) #> 7
# Int64(7.01) #> ERROR: InexactError: Int64(Int64, 7.01)
# convert(Int64, "CV")
# ERROR: MethodError: Cannot `convert` an object of type String to an object of type Int64
# convert(Int64, "123") # same ERROR as above

promote(1, 2.5, 3//4) #> (1.0,2.5,0.75)
promote(1.5, im) #> (1.5 + 0.0im,0.0 + 1.0im)
promote(true, 1.0) #> (1.0, 1.0)
promote_type(Int8, UInt16) #> UInt16


text:
using Test

import sbp
using sbp.Operator: read_operator, to_matrix
using sbp.Grid:grid_xp, grid_xm
using sbp.Test:test_first_derivative, test_quadrature, test_interpolation,
        test_first_derivative_sbp, test_interpolation_sbp
include("operators.jl")

using . OP2019: build_operators, load_description

m = 10

desc = load_description()
xp, xm, h, Dp, Dm, Hp, Hm, Pp, Pm = build_operators(desc, m)

@time begin

@testset "Dp_42" begin
        ml = size(desc.Dp.left, 1)
        mr = size(desc.Dp.right, 1)
        test_first_derivative(Dp, xp, xm, 4, 2, ml, mr)
end

@testset "Dm_42" begin
        ml = size(desc.Dm.left, 1)
        mr = size(desc.Dm.right, 1)
        test_first_derivative(Dm, xm, xp, 4, 2, ml, mr)
end

@testset "Pp_42" begin
        ml = size(desc.Pp.left, 1)
        mr = size(desc.Pp.right, 1)
        test_interpolation(Pp, xp, xm, 3, 1, ml, mr)
end

@testset "Pm_42" begin
        ml = size(desc.Pm.left, 1)
        mr = size(desc.Pm.right, 1)
        test_interpolation(Pm, xm, xp, 3, 1, ml, mr)
end

@testset "Hp_42" begin
        test_quadrature(Hp, xp, 3)

end

@testset "Hm_42" begin
        test_quadrature(Hm, xm, 3)

end

@testset "D_SBP" begin
        B = 0*Dp
        B[1,1] = -1
        B[end,end] = 1
        test_first_derivative_sbp(Hp, Dp, Hm, Dm, B)
end

@testset "P_SBP" begin
        test_interpolation_sbp(Hp, Pp, Hm, Pm)
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
"""
    LCS(a, b)

# Brief:
    A function that finds the longest common subsequence of two strings a and b.

# For more informations:
    https://www.programiz.com/dsa/longest-common-subsequence

# Contributed by:- [Nikola Mircic](https://github.com/Nikola-Mircic)
"""

function LCS(a::String, b::String)
    # Create a matrix to store previous results
    # sol[i+1, j+1] is length of lcs for a[1:i] and b[1:j]
    sol = zeros(length(a)+1, length(b)+1)

    # Calculate sol
    for i in 1:length(a)
        for j in 1:length(b)
            if a[i] == b[j]
                # If a character a[i] is equal to a character b[j],
                # new length of lcs will be equal to length without it plus one
                sol[i+1, j+1] = sol[i, j]+1
            else
                # If characters are different, length will be maximum of two previously calculated lcs
                # One is without a[i], and second one is lcs when b[i] is not inluded
                sol[i+1, j+1] = max(sol[i+1, j], sol[i, j+1])
            end
        end
    end

    result = ""

    # Start going trough the matrix from the last value and check which character
    # increases the length of lcs, and witch is copied value
    i = length(a)+1
    j = length(b)+1
    while i>1 && j>1
        if sol[i, j] == sol[i, j-1]
            # Copied left value
            j-=1
        elseif sol[i, j] == sol[i-1, j]
            # Copied value above
            i-=1
        else
            # Value is increased here
            result = b[j-1] * result # or a[i-1] * result
            i-=1
            j-=1
        end
    end

    return result
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
#
# Test BayesNet generation
#

let
N = 10
mxpa = 4
mxst = 5

bn = rand_discrete_bn(N, mxpa, mxst)

for i = 1:N
    s = Symbol("N", i)
    @test s in names(bn)
end

for nn in names(bn)
    @test ncategories(bn, nn) <= mxst
    @test length(parents(bn, nn)) <= mxpa
end

@test_throws ArgumentError rand_bn_inference(bn, 20, 16)

nq = 3
ne = 5

inf = rand_bn_inference(bn, nq, ne)
@test isempty(intersect(inf.query, keys(inf.evidence)))
end


text:
using Test
using MixedModelsSim
using Random
using DSP
include("test_utilities.jl")
τ = 1.5
fs = 18
evt, epoch_dat = simulate_lmm(
    MersenneTwister(1),
    τ,
    fs,
    β = [0.0, -1.0],
    σs = [1, 1, 1],
    σ = 1,
    n_sub = 20,
    n_item = 30,
    noise_type = "AR-exponential",
);

#mres,res = fit(UnfoldLinearMixedModel,f,evt,dat[chIx:chIx,:,:].*1e6 ,times,contrasts=Dict(:category => EffectsCoding(), :condition => EffectsCoding()));

times = range(0, stop = τ - 1 / fs, step = 1 / fs)

f2 = @formula 0 ~ 1 + stimType + (1 + stimType | subj) + (1 | item)
mres, res = fit(UnfoldLinearMixedModel, f2, evt, epoch_dat, times)

p_df = cluster_permutation_test(mres, epoch_dat, times, 2)

@test p_df[1, :pval] == 0.036

#--- Testing cluster_permutation call
tRange = 1:length(times)
nPerm = 10
permDat = cluster_permutation(mres, epoch_dat, tRange, 2, nPerm)
@test size(permDat, 2) == nPerm
@test size(permDat, 1) == length(tRange)

tRange = 5:length(times)-5
permDat = cluster_permutation(mres, epoch_dat, tRange, 2, nPerm)
@test size(permDat, 1) == length(tRange)

#------ testing pymne_cluster
z_obs = [
    m.z for m in coefpvalues(mres.modelinfo) if
    String(m.coefname) == coefnames(mres.X.formulas)[2][coeffOfInterest]
]

#-- clusterpermutation
obs_cluster = pymne_cluster(z_obs, 4)
@test all(abs.(obs_cluster[2]) .> 4) # clusters have to be larger 4
@test obs_cluster[1][1][1].start == 13 # start of cluster
@test obs_cluster[1][1][1].stop == 14 # end of cluster

obs_cluster = pymne_cluster(z_obs, 0.1)
@test all(abs.(obs_cluster[2]) .> 0.1) # all clusters have to be larger 0.1
@test length(obs_cluster[2]) == 12 # if nothing changed, 12 clusters should be found

#-- tfce
clusterFormingThreshold = Dict(:start => 0, :step => 0.2)
obs_cluster = pymne_cluster(z_obs, clusterFormingThreshold)
@test all(obs_cluster[2] .>= 0) # TFCE is always positive
@test findmax(obs_cluster[2])[2] == findmax(abs.(z_obs))[2] # tfce should not move maximum

#gen_han(τ,fs,3)
if 1 == 0
    # plot some design bits
    map = mapping(:subj, :dv, color = :stimType)#,linestyle=:channel)
    AlgebraOfGraphics.data(evt) * visual(Scatter) * map |> draw
end

max_forks_repo_name: unfoldtoolbox/jlunfold