text:
@doc raw"""
    Direction2

An object `d` of the class [`Direction_2`](@ref) is a vector in the
two-dimensional vector space ``ℝ^2`` where we forget about its length.

They can be viewed as unit vectors, although there is no normalization
internally, since this is error prone. Directions are used whenever the
length of a vector does not matter. They also characterize a set of parallel
oriented lines that have the same orientations. For example, you can ask for
the direction orthogonal to an oriented plane, or the direction of an
oriented line. Further, they can be used to indicate angles. The slope of a
direction is [`dy()`](@ref)/[`dx()`](@ref).

There is a total order on directions.

We compare the angles between the positive ``x``-axis and the directions in
counterclockwise order.
"""
Direction2

"""
    Direction2(v::Vector2)

Introduces the direction `d` of vector `v`.
"""
Direction2(v::Vector2)

"""
    Direction2(l::Line2)

Introduces the direction `d` of line `l`.
"""
Direction2(l::Line2)

"""
    Direction2(r::Ray2)

Introduces the direction `d` of ray `r`.
"""
Direction2(r::Ray2)

"""
    Direction2(s::Segment2)

Introduces the direction `d` of segment `s`.
"""
Direction2(s::Segment2)

"""
    Direction2(x::Real, y::Real)

Introduces a direction `d` passing through the origin and the point with
Cartesian coordinates ``(x, y)``.
"""
Direction2(x::Real, y::Real) = Direction2(convert(FT, x), convert(FT, y))

@doc raw"""
    delta(d::Direction2, i::Integer)

Returns values, such that `d == [Direction2](@ref)(delta(d, 0), delta(d, 1))`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``
"""
delta(d::Direction2, i::Integer)

"""
    dx(d::Direction2)

Returns `delta(d, 0)`.
"""
dx(d::Direction2)

"""
    dy(d::Direction2)

Returns `delta(d, 1)`.
"""
dy(d::Direction2)

"""
    counterclockwise_in_between(d::Direction2, d1::Direction2, d2::Direction2)

Returns `true`, iff `d` is not equal to `d1`, and while rotating
counterclockwise starting at `d1`, `d` is reached strictly before `d2` is
reached.

Note that `true` is returned if `d1 == d2`, unless also `d == d1`.
"""
counterclockwise_in_between(d::Direction2, d1::Direction2, d2::Direction2)

"""
    -(d::Direction2)

The direction opposite to `d`.
"""
-(d::Direction2)

"""
    vector(d::Direction2)

Returns a vector that has the same direction as `d`.
"""
vector(d::Direction2)

"""
    transform(d::Direction2, t::AffTransformation2)

Returns the direction obtained by applying `t` on `d`.
"""
transform(d::Direction2, t::AffTransformation2)


text:
function layout_text(
        string::AbstractString, startpos::VecTypes{N, T}, textsize::Number,
        font, align, rotation, model
    ) where {N, T}

    offset_vec = to_align(align)
    ft_font = to_font(font)
    rscale = to_textsize(textsize)
    rot = to_rotation(rotation)

    atlas = get_texture_atlas()
    mpos = model * Vec4f0(to_ndim(Vec3f0, startpos, 0f0)..., 1f0)
    pos = to_ndim(Point{N, Float32}, mpos, 0)
    scales = Vec2f0[glyph_scale!(atlas, c, ft_font, rscale) for c in string]

    positions2d = calc_position(string, Point2f0(0), rscale, ft_font, atlas)
    aoffset = align_offset(Point2f0(0), positions2d[end], atlas, rscale, ft_font, offset_vec)
    aoffsetn = to_ndim(Point{N, Float32}, aoffset, 0f0)
    toffset = calc_offset(string, scales, ft_font, atlas)
    positions = map(positions2d, toffset) do p, o
        pn = rot * (to_ndim(Point{N, Float32}, p .+ o, 0f0) .+ aoffsetn)
        pn .+ (pos)
    end
    positions, scales
end
function text_bb(str, font, size)
    positions, scale = layout_text(
        str, Point2f0(0), size,
        font, Vec2f0(0), Quaternionf0(0,0,0,1), Mat4f0(I)
    )
    union(AABB(positions),  AABB(positions .+ scale))
end

"""
calculates how much `child` rectangle needs to move to not touch the `parent`
"""
function move_from_touch(
        parent::GeometryPrimitive{N, T}, child::GeometryPrimitive{N},
        pad::Vec{N}
    ) where {N, T}
    pmini, cmini = minimum(parent), minimum(child) .- pad
    pmaxi, cmaxi = maximum(parent), maximum(child) .+ pad

    move = ntuple(Val(N)) do i
        posdir = ifelse(cmini[i] < pmini[i], (pmini[i] - cmini[i]), zero(T)) #always positive
        negdir = ifelse(cmaxi[i] > pmaxi[i], (pmaxi[i] - cmaxi[i]), zero(T)) #always minus
        ifelse(posdir > abs(negdir), posdir, negdir) # move in the bigger direction
    end
    Vec{N, T}(move)
end

"""
Moves `child` so that it doesn't touch parent. Leaves a gap to parent defined by `pad`.
"""
function dont_touch(
        parent::GeometryPrimitive{N}, child::GeometryPrimitive{N},
        pad::Vec{N}
    ) where N
    child + move_from_touch(parent, child, pad)
end

"""
    fit_factor_stretch(rect, lims::NTuple{N}) where N

Calculates the stretch factor to fill `rect` in all dimension.
Returns a stretch `N` dimensional fit factor.
"""
function fit_factor_stretch(rect, lims::NTuple{N, Any}) where N
    w = widths(rect)
    stretches = ntuple(Val(N)) do i
        from, to = lims[i]
        w[i] / abs(to - from)
    end
    stretches
end

"""
    fit_factor(rect, lims::NTuple{N}) where N

Calculates the scaling one needs to apply to lims to fit `rect` without changing aspect ratio.
Returns float scaling and the full strech as given by [`fit_factor_stretch`](@ref)
"""
function fit_factor(rect, lims::NTuple{N, Any}) where N
    stretches = fit_factor_stretch(rect, lims)
    minimum(stretches), stretches
end

"""
    fit_ratio(rect, lims)

Calculates the ratio one needs to stretch `lims` in order to get the same aspect ratio
"""
function fit_ratio(rect, lims)
    s, stretches = fit_factor(rect, lims)
    stretches ./ s
end

function align_offset(startpos, lastpos, atlas, rscale, font, align)
    xscale, yscale = glyph_scale!('X', rscale)
    xmove = (lastpos-startpos)[1] + xscale
    if isa(align, GeometryTypes.Vec)
        return -Vec2f0(xmove, yscale) .* align
    elseif align == :top
        return -Vec2f0(xmove/2f0, yscale)
    elseif align == :right
        return -Vec2f0(xmove, yscale/2f0)
    else
        error("Align $align not known")
    end
end

function alignment2num(x::Symbol)
    (x == :center) && return 0.5f0
    (x in (:left, :bottom)) && return 0.0f0
    (x in (:right, :top)) && return 1.0f0
    0.0f0 # 0 default, or better to error?
end

grid(x::Transformable...; kw_args...) = grid([x...]; kw_args...)
function grid(plots::Vector{<: Transformable}; kw_args...)
    N = length(plots)
    grid = close2square(N)
    grid(reshape(plots, grid))
end
function grid(plots::Matrix{<: Transformable}; kw_args...)
    N = length(plots)
    grid = size(plots)
    w, h = (0.0, 0.0)
    pscene = Scene()
    cam2d!(pscene)
    for idx in 1:N
        i, j = Tuple(CartesianIndices(grid)[idx])
        child = plots[idx]
        push!(pscene.children, child)
        translate!(child, w, h, 0.0)
        swidth = widths(boundingbox(child))
        if i == grid[1]
            h += (swidth[2] * 1.1)
            w = 0.0
        else
            w += (swidth[1] * 1.1)
        end
    end
    center!(pscene)
    pscene
end

vbox(plots::Transformable...; kw_args...) = vbox([plots...]; kw_args...)

estimated_space(x, N, w) = 1/N

function Base.resize!(scene::Scene, rect::Rect2D)
    pixelarea(scene)[] = rect
    force_update!()
    yield()
end

function vbox(plots::Vector{T}; kw_args...) where T <: Scene
    N = length(plots)
    w = 0.0
    pscene = Scene()
    area = pixelarea(pscene)

    for idx in 1:N
        p = plots[idx]
        foreach(area) do a
            # TODO this is terrible!
            w2 = widths(a) .* Vec((1/N), 1)
            resize!(p, IRect(minimum(a) .+ Vec((idx - 1) * w2[1], 0), w2))
            center!(p)
        end
        push!(pscene.children, p)
        nodes = map(fieldnames(Events)) do field
            if field != :window_area
                foreach(getfield(pscene.events, field)) do val
                    push!(getfield(p.events, field), val)
                end
            end
        end
    end
    pscene
end

function vbox(plots::Vector{T}; kw_args...) where T <: AbstractPlot
    N = length(plots)
    w = 0.0
    for idx in 1:N
        p = plots[idx]
        translate!(p, w, 0.0, 0.0)
        swidth = widths(boundingbox(p))
        w += (swidth[1] * 1.1)
    end
end

function hbox!(plots::Vector{T}; kw_args...) where T <: AbstractPlot
    N = length(plots)
    h = 0.0
    for idx in 1:N
        p = plots[idx]
        translate!(p, 0.0, h, 0.0)
        swidth = widths(boundingbox(p))
        h += (swidth[2] * 1.2)
    end
end

ext: jl


text:
module Constant
const NONO_SIZE = 25
const ATTEMPT_TRY = 100000
const NEW_GEN = 10
const PAIRS = 10
const PROBLEMS = 1
end


text:
using Documenter, NLIDatasets

makedocs(sitename = "NLIDatasets.jl",
         format = Documenter.HTML(),
         modules = [NLIDatasets],
         pages = ["Home" => "index.md"],
         doctest = true)

deploydocs(repo = "github.com/dellison/NLIDatasets.jl.git")


text:
using HoJBot: replace_newlines, is_real_user, is_bot, flipper, previous, next
using Discord: User

@testset "utilities" begin
    let
        @test replace_newlines("a\nb") == "a b"
        @test replace_newlines("a\nb\n") == "a b "
        @test replace_newlines("a\nb \n\n") == "a b   "
    end

    let user = User(; id=1, bot=false), bot = User(; id=2, bot=true)
        @test is_real_user(user) === true
        @test is_real_user(bot) === false
        @test is_bot(user) === false
        @test is_bot(bot) === true
    end

    # Tainted module here but it's OK since it's just a test
    HoJBot.previous(x::Int) = x - 1
    HoJBot.next(x::Int) = x + 1
    @test flipper(previous, Int)(1) == 0
    @test flipper(next, Int)(1) == 2
end


text:
using Sesam
using Documenter

DocMeta.setdocmeta!(Sesam, :DocTestSetup, :(using Sesam); recursive=true)

makedocs(;
    modules=[Sesam],
    authors="Thomas Wutzler <twutz@bgc-jena.mpg.de> and contributors",
    repo="https://github.com/bgctw/Sesam.jl/blob/{commit}{path}#{line}",
    sitename="Sesam.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/Sesam.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bgctw/Sesam.jl",
    devbranch="main",
)


text:
"""
    OpticalInterface{T<:Real}

Any subclass of OpticalInterface **must** implement the following:

```julia
processintersection(opticalinterface::OpticalInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, ::Bool, firstray::Bool = false) -> Tuple{SVector{N,T}, T, T}
```

See documentation for [`processintersection`](@ref) for details.

These methods are also commonly implemented, but not essential:
```julia
insidematerialid(i::OpticalInterface{T}) -> OpticSim.GlassCat.AbstractGlass
outsidematerialid(i::OpticalInterface{T}) -> OpticSim.GlassCat.AbstractGlass
reflectance(i::OpticalInterface{T}) -> T
transmission(i::OpticalInterface{T}) -> T
```
"""
abstract type OpticalInterface{T<:Real} end
export OpticalInterface

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

"""
    NullInterface{T} <: OpticalInterface{T}

Interface which will be ignored totally by any rays, used only in construction of CSG objects.

```julia
NullInterface(T = Float64)
NullInterface{T}()
```
"""
struct NullInterface{T} <: OpticalInterface{T}
    NullInterface(::Type{T} = Float64) where {T<:Real} = new{T}()
    NullInterface{T}() where {T<:Real} = new{T}()
end

insidematerialid(::NullInterface{T}) where {T<:Real} = glassid(OpticSim.GlassCat.Air)
outsidematerialid(::NullInterface{T}) where {T<:Real} = glassid(OpticSim.GlassCat.Air)
reflectance(::NullInterface{T}) where {T<:Real} = zero(T)
transmission(::NullInterface{T}) where {T<:Real} = one(T)

export NullInterface

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

"""
    ParaxialInterface{T} <: OpticalInterface{T}

Interface describing an idealized planar lens, i.e. one that is thin and with no aberrations.

**In general this interface should not be constructed directly, the `ParaxialLensEllipse` and `ParaxialLensRect` functions should be used to create a [`ParaxialLens`](@ref) object directly.**

```julia
ParaxialInterface(focallength::T, centroid::SVector{3,T}, outsidematerial::Y)
```
"""
struct ParaxialInterface{T} <: OpticalInterface{T}
    focallength::T
    outsidematerial::GlassID
    centroid::SVector{3,T}
    function ParaxialInterface(focallength::T, centroid::SVector{3,T}, outsidematerial::Y) where {Y<:OpticSim.GlassCat.AbstractGlass,T<:Real}
        return new{T}(focallength, glassid(outsidematerial), centroid)
    end
end
export ParaxialInterface

function Base.show(io::IO, a::ParaxialInterface{R}) where {R<:Real}
    print(io, "ParaxialInterface($(a.focallength), $(glassname(a.outsidematerial)))")
end

insidematerialid(a::ParaxialInterface{T}) where {T<:Real} = a.outsidematerial
outsidematerialid(a::ParaxialInterface{T}) where {T<:Real} = a.outsidematerial
reflectance(::ParaxialInterface{T}) where {T<:Real} = zero(T)
transmission(::ParaxialInterface{T}) where {T<:Real} = one(T)

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

"""
    FresnelInterface{T} <: OpticalInterface{T}

Interface between two materials with behavior defined according to the [Fresnel equations](https://en.wikipedia.org/wiki/Fresnel_equations), with a specified reflectance and transmission.
Assumes unpolarized light.

```julia
FresnelInterface{T}(insidematerial, outsidematerial; reflectance = 0, transmission = 1)
```
"""
struct FresnelInterface{T} <: OpticalInterface{T}
    # storing glasses as IDs (integer) rather than the whole thing seems to improve performance significantly, even when the Glass type is a fixed size (i.e. the interface is pointer-free)
    insidematerial::GlassID
    outsidematerial::GlassID
    reflectance::T
    transmission::T

    function FresnelInterface{T}(insidematerial::Z, outsidematerial::Y; reflectance::T = zero(T), transmission::T = one(T)) where {Z<:OpticSim.GlassCat.AbstractGlass,Y<:OpticSim.GlassCat.AbstractGlass,T<:Real}
        return FresnelInterface{T}(glassid(insidematerial), glassid(outsidematerial), reflectance = reflectance, transmission = transmission)
    end

    function FresnelInterface{T}(insidematerialid::GlassID, outsidematerialid::GlassID; reflectance::T = zero(T), transmission::T = one(T)) where {T<:Real}
        @assert zero(T) <= reflectance <= one(T)
        @assert zero(T) <= transmission <= one(T)
        @assert reflectance + transmission <= one(T)
        return new{T}(insidematerialid, outsidematerialid, reflectance, transmission)
    end
end
export FresnelInterface

function Base.show(io::IO, a::FresnelInterface{R}) where {R<:Real}
    print(io, "FresnelInterface($(glassname(a.insidematerial)), $(glassname(a.outsidematerial)), $(a.reflectance), $(a.transmission))")
end

insidematerialid(a::FresnelInterface{T}) where {T<:Real} = a.insidematerial
outsidematerialid(a::FresnelInterface{T}) where {T<:Real} = a.outsidematerial
reflectance(a::FresnelInterface{T}) where {T<:Real} = a.reflectance
transmission(a::FresnelInterface{T}) where {T<:Real} = a.transmission

transmissiveinterface(::Type{T}, insidematerial::X, outsidematerial::Y) where {T<:Real,X<:OpticSim.GlassCat.AbstractGlass,Y<:OpticSim.GlassCat.AbstractGlass} = FresnelInterface{T}(insidematerial, outsidematerial, reflectance = zero(T), transmission = one(T))
reflectiveinterface(::Type{T}, insidematerial::X, outsidematerial::Y) where {T<:Real,X<:OpticSim.GlassCat.AbstractGlass,Y<:OpticSim.GlassCat.AbstractGlass} = FresnelInterface{T}(insidematerial, outsidematerial, reflectance = one(T), transmission = zero(T))
opaqueinterface(::Type{T} = Float64) where {T<:Real} = FresnelInterface{T}(OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, reflectance = zero(T), transmission = zero(T))
export opaqueinterface

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

"""
    ThinGratingInterface{T} <: OpticalInterface{T}

Interface representing an idealized thin grating. `period` is in microns, `vector` should lie in the plane of the surface.
Transmission and reflectance can be specified for an arbitrary number of orders up to 10, selected using the `maxorder` and `minorder` parameters.
If `nothing` then `reflectance` is assumed to be **0** and `transmission` is assumed to be **1**.

```julia
ThinGratingInterface(vector, period, insidematerial, outsidematerial; maxorder = 1, minorder = -1, reflectance = nothing, transmission = nothing)
```
"""
struct ThinGratingInterface{T} <: OpticalInterface{T}
    insidematerial::GlassID
    outsidematerial::GlassID
    vector::SVector{3,T}
    period::T
    maxorder::Int
    minorder::Int
    transmission::SVector{GRATING_MAX_ORDERS,T}
    reflectance::SVector{GRATING_MAX_ORDERS,T}

    function ThinGratingInterface(vector::SVector{3,T}, period::T, insidematerial::Z, outsidematerial::Y; maxorder::Int = 1, minorder::Int = -1, reflectance::Union{Nothing,AbstractVector{T}} = nothing, transmission::Union{Nothing,AbstractVector{T}} = nothing) where {T<:Real,Z<:OpticSim.GlassCat.AbstractGlass,Y<:OpticSim.GlassCat.AbstractGlass}
        @assert maxorder >= minorder
        norders = maxorder - minorder + 1
        @assert norders <= GRATING_MAX_ORDERS "Thin grating is limited to $GRATING_MAX_ORDERS orders"
        @assert zero(T) <= sum(reflectance) <= one(T)
        @assert zero(T) <= sum(transmission) <= one(T)
        @assert zero(T) <= sum(reflectance) + sum(transmission) <= one(T)
        @assert ((reflectance === nothing) || length(reflectance) == norders) && ((transmission === nothing) || length(transmission) == norders)
        if reflectance !== nothing
            sreflectance = vcat(SVector{length(reflectance),T}(reflectance), ones(SVector{GRATING_MAX_ORDERS - length(reflectance),T}))
        else
            sreflectance = zeros(SVector{GRATING_MAX_ORDERS,T})
        end
        if transmission !== nothing
            stransmission = vcat(SVector{length(transmission),T}(transmission), ones(SVector{GRATING_MAX_ORDERS - length(transmission),T}))
        else
            stransmission = ones(SVector{GRATING_MAX_ORDERS,T})
        end
        new{T}(glassid(insidematerial), glassid(outsidematerial), normalize(vector), period, maxorder, minorder, sreflectance, stransmission)
    end
end
export ThinGratingInterface

insidematerialid(a::ThinGratingInterface{T}) where {T<:Real} = a.insidematerial
outsidematerialid(a::ThinGratingInterface{T}) where {T<:Real} = a.outsidematerial
reflectance(a::ThinGratingInterface{T}, order::Int) where {T<:Real} = a.reflectance[order - a.minorder + 1]
transmission(a::ThinGratingInterface{T}, order::Int) where {T<:Real} = a.transmission[order - a.minorder + 1]

function Base.show(io::IO, a::ThinGratingInterface{R}) where {R<:Real}
    print(io, "ThinGratingInterface($(glassname(a.insidematerial)), $(glassname(a.outsidematerial)), $(a.vector), $(a.period), $(a.transmission), $(a.reflectance))")
end

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

"""
`ConvergingBeam`, `DivergingBeam` or `CollimatedBeam`, defines the behavior of a beam in a [`HologramInterface`](@ref).
"""
@enum BeamState ConvergingBeam DivergingBeam CollimatedBeam
export BeamState, ConvergingBeam, DivergingBeam, CollimatedBeam

"""
    HologramInterface{T} <: OpticalInterface{T}

Interface representing a _thick_ hologram (though geometrically thin).
The efficiency, `η`, is calculated using Kogelnik's coupled wave theory so is only valid for the first order.
If the zero order is included then it has efficiency `1 - η`.
Also assumes that the HOE was recorded under similar conditions to the playback conditions, `thickness` is in microns.

`BeatState` arguments can be one of `ConvergingBeam`, `DivergingBeam` and `CollimatedBeam`. In the first two cases `signalpointordir` and `referencepointordir` are 3D point in global coordinate space. For `CollimatedBeam` they are normalized direction vectors.

For reference, see:
- _Coupled Wave Theory for Thick Hologram Gratings_ - H Kogelnik, 1995
- _Sequential and non-sequential simulation of volume holographic gratings_ - M Kick et al, 2018

```julia
HologramInterface(signalpointordir::SVector{3,T}, signalbeamstate::BeamState, referencepointordir::SVector{3,T}, referencebeamstate::BeamState, recordingλ::T, thickness::T, beforematerial, substratematerial, aftermaterial, signalrecordingmaterial, referencerecordingmaterial, RImodulation::T, include0order  = false)
```
"""
struct HologramInterface{T} <: OpticalInterface{T}
    beforematerial::GlassID
    substratematerial::GlassID
    aftermaterial::GlassID
    signalpointordir::SVector{3,T}
    signalbeamstate::BeamState
    referencepointordir::SVector{3,T}
    referencebeamstate::BeamState
    recordingλ::T
    signalrecordingmaterial::GlassID
    referencerecordingmaterial::GlassID
    thickness::T
    RImodulation::T
    include0order::Bool

    function HologramInterface(signalpointordir::SVector{3,T}, signalbeamstate::BeamState, referencepointordir::SVector{3,T}, referencebeamstate::BeamState, recordingλ::T, thickness::T, beforematerial::Z, substratematerial::X, aftermaterial::Y, signalrecordingmaterial::P, referencerecordingmaterial::Q, RImodulation::T, include0order::Bool = false) where {T<:Real,X<:OpticSim.GlassCat.AbstractGlass,Z<:OpticSim.GlassCat.AbstractGlass,Y<:OpticSim.GlassCat.AbstractGlass,P<:OpticSim.GlassCat.AbstractGlass,Q<:OpticSim.GlassCat.AbstractGlass}
        @assert substratematerial !== OpticSim.GlassCat.Air "Substrate can't be air"
        if signalbeamstate === CollimatedBeam
            signalpointordir = normalize(signalpointordir)
        end
        if referencebeamstate === CollimatedBeam
            referencepointordir = normalize(referencepointordir)
        end
        new{T}(glassid(beforematerial), glassid(substratematerial), glassid(aftermaterial), signalpointordir, signalbeamstate, referencepointordir, referencebeamstate, recordingλ, glassid(signalrecordingmaterial), glassid(referencerecordingmaterial), thickness, RImodulation, include0order)
    end

    # 'zero' constructor to use for blank elements in SVector - not ideal...
    HologramInterface{T}() where {T<:Real} = new{T}(NullInterface(T), 0, NullInterface(T), zeros(SVector{3,T}), CollimatedBeam, zeros(SVector{3,T}), CollimatedBeam, zero(T), 0, 0, zero(T), zero(T), false)
end
export HologramInterface

insidematerialid(a::HologramInterface{T}) where {T<:Real} = a.beforematerial
outsidematerialid(a::HologramInterface{T}) where {T<:Real} = a.aftermaterial
substratematerialid(a::HologramInterface{T}) where {T<:Real} = a.substratematerial

function Base.show(io::IO, a::HologramInterface{R}) where {R<:Real}
    print(io, "HologramInterface($(glassname(a.beforematerial)), $(glassname(a.substratematerial)), $(glassname(a.aftermaterial)), $(a.signalpointordir), $(a.signalbeamstate), $(a.referencepointordir), $(a.referencebeamstate), $(a.recordingλ), $(a.thickness), $(a.RImodulation))")
end

"""
    MultiHologramInterface{T} <: OpticalInterface{T}

Interface to represent multiple overlapped [`HologramInterface`](@ref)s on a single surface. Each ray randomly selects an interface to use.

```julia
MultiHologramInterface(interfaces::Vararg{HologramInterface{T}})
MultiHologramInterface(interfaces::Vector{HologramInterface{T}})
```
"""
struct MultiHologramInterface{T} <: OpticalInterface{T}
    interfaces::Ptr{HologramInterface{T}} # pointer to the array of hologram interfaces TODO!! fix this to not be hacky...
    numinterfaces::Int

    MultiHologramInterface(interfaces::Vararg{HologramInterface{T},N}) where {T<:Real,N} = MultiHologramInterface(collect(interfaces))

    function MultiHologramInterface(interfaces::Vector{HologramInterface{T}}) where {T<:Real}
        N = length(interfaces)
        @assert N > 1 "Don't need to use MultiHologramInterface if only 1 hologram"
        p = pointer(interfaces)
        push!(MultiHologramInterfaceRefCache, interfaces)
        new{T}(p, N)
    end
end
export MultiHologramInterface
interface(m::MultiHologramInterface{T}, i::Int) where {T<:Real} = (@assert i <= m.numinterfaces; return unsafe_load(m.interfaces, i))

# need this to keep julia references to the arrays used in MultiHologramInterface so they don't get garbage collected
const MultiHologramInterfaceRefCache = []

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

# a few things need a concrete type to prevent allocations resulting from the ambiguities introduce by an abstract type (i.e. OpticalInterface{T})
# if adding new subtypes of OpticalInterface they must be added to this definition as well (and only paramaterised by T so they are fully specified)
AllOpticalInterfaces{T} = Union{NullInterface{T},ParaxialInterface{T},FresnelInterface{T},HologramInterface{T},MultiHologramInterface{T},ThinGratingInterface{T}}
NullOrFresnel{T} = Union{NullInterface{T},FresnelInterface{T}}

# can't have more than 4 types here or we get allocations because inference on the union stops: https://github.com/JuliaLang/julia/blob/1e6e65691254a7fe81f5da8706bb30aa6cb3f8d2/base/compiler/types.jl#L113
# use the macro to get around this: @unionsplit OpticalInterface T x f(x)
import InteractiveUtils
# TODO REMOVE
macro unionsplit(type, paramtype, var, call)
    st = InteractiveUtils.subtypes(@eval $type)
    orig_expr = Expr(:if, :($var isa $(st[1]){$paramtype}), call)
    expr = orig_expr
    for t in st[2:end]
        push!(expr.args, Expr(:elseif, :($var isa $t{$paramtype}), call))
        expr = expr.args[end]
    end
    push!(expr.args, :(error("unhandled type")))
    return :($(esc(orig_expr)))
end


text:
sizeof("Hello, world!") # gives 13
sizeof("Hellö, wørld!") # gives 15

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function parareal_distributed!(cache::PararealCache, solution::PararealSolution, problem, solver::Parareal; savechunks::Bool = false)
    @↓ iterates, ψ, U, T = solution
    @↓ ℱ, 𝒢, P, K = solver
    @↓ 𝜑, ϵ, Λ, updateΛ = solver.error_check
    # coarse guess
    G = similar(U)
    G[1] = U[1]
    for n = 1:P
        chunk = 𝒢(problem, U[n], T[n], T[n+1])
        G[n+1] = chunk.u[end]
    end
    # main loop
    F = similar(U)
    F[1] = U[1]
    getF(args...) = ℱ(args...).u[end]
    for k = 1:K
        # @↑ solution[k] = U .← U
        solution[k].U .= U
        # for n = 1:k-1
        #     solution[k][n] = solution[k-1][n]
        # end
        # fine run (with Julia's Distributed.jl)
        @sync for n = k:P
            @async F[n+1] = remotecall_fetch(getF, n, problem, U[n], T[n], T[n+1])
        end
        solution[k].F .= F
        # update Lipschitz constant
        Λ = updateΛ ? update_Lipschitz(Λ, U, F) : Λ
        # check convergence
        ψ[k] = 𝜑(solution, k, Λ)
        if ψ[k] ≤ ϵ
            resize!(iterates, k)
            resize!(ψ, k)
            break
        end
        # update (serial)
        for n = k:P
            chunk = 𝒢(problem, U[n], T[n], T[n+1])
            U[n+1] = chunk.u[end] + F[n+1] - G[n+1]
            G[n+1] = chunk.u[end]
        end
        @↑ solution = U, T
    end
    return solution
end


text:
const pyjlnumbertype = pynew()
const pyjlcomplextype = pynew()
const pyjlrealtype = pynew()
const pyjlrationaltype = pynew()
const pyjlintegertype = pynew()

struct pyjlnumber_op{OP}
    op :: OP
end
(op::pyjlnumber_op)(self) = Py(op.op(self))
function (op::pyjlnumber_op)(self, other_::Py)
    if pyisjl(other_)
        other = pyjlvalue(other_)
    else
        other = @pyconvert(Number, other_, return Py(pybuiltins.NotImplemented))
    end
    Py(op.op(self, other))
end
function (op::pyjlnumber_op)(self, other_::Py, other2_::Py)
    if pyisjl(other_)
        other = pyjlvalue(other_)
    else
        other = @pyconvert(Number, other_, return Py(pybuiltins.NotImplemented))
    end
    if pyisjl(other2_)
        other2 = pyjlvalue(other2_)
    else
        other2 = @pyconvert(Number, other2_, return Py(pybuiltins.NotImplemented))
    end
    Py(op.op(self, other, other2))
end
pyjl_handle_error_type(op::pyjlnumber_op, self, exc) = exc isa MethodError && exc.f === op.op ? pybuiltins.TypeError : PyNULL

struct pyjlnumber_rev_op{OP}
    op :: OP
end
function (op::pyjlnumber_rev_op)(self, other_::Py)
    if pyisjl(other_)
        other = pyjlvalue(other_)
    else
        other = @pyconvert(Number, other_, return Py(pybuiltins.NotImplemented))
    end
    Py(op.op(other, self))
end
function (op::pyjlnumber_rev_op)(self, other_::Py, other2_::Py)
    if pyisjl(other_)
        other = pyjlvalue(other_)
    else
        other = @pyconvert(Number, other_, return Py(pybuiltins.NotImplemented))
    end
    if pyisjl(other2_)
        other2 = pyjlvalue(other2_)
    else
        other2 = @pyconvert(Number, other2_, return Py(pybuiltins.NotImplemented))
    end
    Py(op.op(other, self, other2))
end
pyjl_handle_error_type(op::pyjlnumber_rev_op, self, exc) = exc isa MethodError && exc.f === op.op ? pybuiltins.TypeError : PyNULL

pyjlreal_trunc(self::Real) = Py(trunc(Integer, self))
pyjl_handle_error_type(::typeof(pyjlreal_trunc), self, exc::MethodError) = exc.f === trunc ? pybuiltins.TypeError : PyNULL

pyjlreal_floor(self::Real) = Py(floor(Integer, self))
pyjl_handle_error_type(::typeof(pyjlreal_floor), self, exc::MethodError) = exc.f === floor ? pybuiltins.TypeError : PyNULL

pyjlreal_ceil(self::Real) = Py(ceil(Integer, self))
pyjl_handle_error_type(::typeof(pyjlreal_ceil), self, exc::MethodError) = exc.f === ceil ? pybuiltins.TypeError : PyNULL

function pyjlreal_round(self::Real, ndigits_::Py)
    ndigits = pyconvertarg(Union{Int,Nothing}, ndigits_, "ndigits")
    if ndigits === nothing
        Py(round(Integer, self))
    else
        Py(round(self; digits = ndigits))
    end
end
pyjl_handle_error_type(::typeof(pyjlreal_round), self, exc::MethodError) = exc.f === round ? pybuiltins.TypeError : PyNULL

function init_jlwrap_number()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class NumberValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __bool__(self):
            return not self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(iszero))))
        def __add__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(+))), other)
        def __sub__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(-))), other)
        def __mul__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(*))), other)
        def __truediv__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(/))), other)
        def __floordiv__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(÷))), other)
        def __mod__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(%))), other)
        def __pow__(self, other, modulo=None):
            if modulo is None:
                return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(^))), other)
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(powermod))), other, modulo)
        def __lshift__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(<<))), other)
        def __rshift__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(>>))), other)
        def __and__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(&))), other)
        def __xor__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(⊻))), other)
        def __or__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(|))), other)
        def __radd__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(+))), other)
        def __rsub__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(-))), other)
        def __rmul__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(*))), other)
        def __rtruediv__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(/))), other)
        def __rfloordiv__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(÷))), other)
        def __rmod__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(%))), other)
        def __rpow__(self, other, modulo=None):
            if modulo is None:
                return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(^))), other)
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(powermod))), other, modulo)
        def __rlshift__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(<<))), other)
        def __rrshift__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(>>))), other)
        def __rand__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(&))), other)
        def __rxor__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(⊻))), other)
        def __ror__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_rev_op(|))), other)
        def __eq__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(==))), other)
        def __ne__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(!=))), other)
        def __le__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(≤))), other)
        def __lt__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(<))), other)
        def __ge__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(≥))), other)
        def __gt__(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(>))), other)
    class ComplexValue(NumberValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __complex__(self):
            return self._jl_callmethod($(pyjl_methodnum(pycomplex)))
        @property
        def real(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(real))))
        @property
        def imag(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(imag))))
        def conjugate(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(conj))))
    class RealValue(ComplexValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __float__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyfloat)))
        @property
        def real(self):
            return self
        @property
        def imag(self):
            return 0
        def conjugate(self):
            return self
        def __complex__(self):
            return complex(float(self))
        def __trunc__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlreal_trunc)))
        def __floor__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlreal_floor)))
        def __ceil__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlreal_ceil)))
        def __round__(self, ndigits=None):
            return self._jl_callmethod($(pyjl_methodnum(pyjlreal_round)), ndigits)
    class RationalValue(RealValue):
        __slots__ = ()
        __module__ = "juliacall"
        @property
        def numerator(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(numerator))))
        @property
        def denominator(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlnumber_op(denominator))))
    class IntegerValue(RationalValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __int__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyint)))
        def __index__(self):
            return self.__int__()
        @property
        def numerator(self):
            return self
        @property
        def denominator(self):
            return 1
    import numbers
    numbers.Number.register(NumberValue)
    numbers.Complex.register(ComplexValue)
    numbers.Real.register(RealValue)
    numbers.Rational.register(RationalValue)
    numbers.Integral.register(IntegerValue)
    del numbers
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlnumbertype, jl.NumberValue)
    pycopy!(pyjlcomplextype, jl.ComplexValue)
    pycopy!(pyjlrealtype, jl.RealValue)
    pycopy!(pyjlrationaltype, jl.RationalValue)
    pycopy!(pyjlintegertype, jl.IntegerValue)
end

pyjl(v::Number) = pyjl(pyjlnumbertype, v)
pyjl(v::Complex) = pyjl(pyjlcomplextype, v)
pyjl(v::Real) = pyjl(pyjlrealtype, v)
pyjl(v::Rational) = pyjl(pyjlrationaltype, v)
pyjl(v::Integer) = pyjl(pyjlintegertype, v)


text:
__precompile__()

module Xtensor

using CxxWrap

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, package xtensor-julia did not build properly")
end
include(depsfile)

end


text:
using StaticThreadSchedules
using Test

@show Threads.nthreads()

@testset "StaticThreadSchedules.jl" begin

    # divisible by common thread counts
    v = zeros(3 * 5 * 128)
    @threads :static (i->log(i+1)) for i in 1:length(v)
        v[i] = Threads.threadid()
    end
    @test all(v .> 0)

    # a prime number
    v = zeros(71)
    @threads :static (i->log(i+1)) for i in 1:length(v)
        v[i] = Threads.threadid()
    end
    @test all(v .> 0)

    # size 1 vector edge case
    v = zeros(1)
    @threads :static (i->log(i+1)) for i in 1:length(v)
        v[i] = Threads.threadid()
    end
    @test all(v .> 0)

    # call Base.@threads if no third arg
    v = zeros(7)
    @threads :static for i in 1:length(v)
        v[i] = Threads.threadid()
    end
    @test all(v .> 0)
    v = zeros(7)
    @threads for i in 1:length(v)
        v[i] = Threads.threadid()
    end
    @test all(v .> 0)

end


text:
noext_options = (; use_extended_form = false)
noconic_options =
    (; solve_subproblems = false, solve_relaxation = false, use_init_fixed_oa = true)
noext_noconic_options = merge(noext_options, noconic_options)

insts = OrderedDict()

insts["test"] = [
    ((3, VecNegRtdet(), false),),
    ((3, VecNegRtdet(), true),),
    ((3, VecNegRtdet(), true), noconic_options),
    ((3, VecNegRtdet(), false), noext_options),
    ((3, VecNegRtdet(), false), noext_noconic_options),
    ((3, VecNegRtdetEFExp(), false),),
    ((5, VecNegSqrtConj(), false),),
    ((5, VecNegSqrtConj(), true),),
    ((5, VecNegSqrtConjEF(), false),),
]

# continuous, separation only:
function cont_knapsack_insts(f::VecSpecExt, max_n::Int, options::NamedTuple)
    return [((n, f, true), options) for n in vcat(5, 5:5:max_n)]
end

insts["cont_geo"] = cont_knapsack_insts(VecNegRtdet(), 200, noconic_options)
insts["cont_noext_geo"] = cont_knapsack_insts(VecNegRtdet(), 70, noext_noconic_options)
insts["cont_log"] = cont_knapsack_insts(VecNegLog(), 200, noconic_options)
insts["cont_noext_log"] = cont_knapsack_insts(VecNegLog(), 70, noext_noconic_options)
insts["cont_inv"] = cont_knapsack_insts(VecNegSqrtConj(), 200, noconic_options)
insts["cont_noext_inv"] = cont_knapsack_insts(VecNegSqrtConj(), 70, noext_noconic_options)

# integer:
function int_knapsack_insts(f::VecSpecExt, max_n::Int, options::NamedTuple)
    return [((n, f, false), options) for n in vcat(3, 3:3:15, 20:5:35, 40:10:max_n)]
end

insts["nat_geo"] = int_knapsack_insts(VecNegRtdet(), 120, (;))
insts["noext_geo"] = int_knapsack_insts(VecNegRtdet(), 20, noext_options)
insts["ext_geo"] = int_knapsack_insts(VecNegRtdetEFExp(), 120, (;))

return (Knapsack, insts)


text:
function createrandoms(n)
  i = 1
  a = zeros(n)
  while (i <= n)
    r1 = rand()
    r2 = rand()
    if (r2 < asin(sqrt(r1)) / (pi / 2))
      a[i] = r1
      i += 1
    end
  end
  return a
end

function tofile(a)
  outfile = open("generated_rand_numbers", "w")
  for x in a
    write(outfile, join(x), "\n")
  end
  close(outfile)
end

function fromfile()
  infile = open("generated_rand_numbers", "r")
  s = 0
  while !eof(infile)
    s += float(readline(infile))
  end
  close(infile)
  return s
end

n = 2000000

print("io.jl\n")

print("generation... ")
@time (a = createrandoms(n))

print("o time... ")
@time (tofile(a))

print("i time... ")
@time (s = fromfile())
print("result = ", s / n)
print('\n')

max_stars_repo_name: EvgenyNerush/performance


text:
# Unit tests for judiVector
# Philipp Witte (pwitte.slim@gmail.com)
# May 2018
#

using JUDI.TimeModeling, SegyIO, Test, LinearAlgebra
import LinearAlgebra.BLAS.axpy!

function example_rec_geometry(; nsrc=2, nrec=120)
    xrec = range(50f0, stop=1150f0, length=nrec)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nrec)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

# number of sources/receivers
nsrc = 2
nrec = 120
ns = 251

################################################# test constructors ####################################################

@testset "judiVector Unit Tests" begin

    # set up judiVector fr,om array
    dsize = (nsrc*nrec*ns, 1)
    rec_geometry = example_rec_geometry(nsrc=nsrc, nrec=nrec)
    data = randn(Float32, ns, nrec)
    d_obs = judiVector(rec_geometry, data)

    @test isequal(d_obs.nsrc, nsrc)
    @test isequal(typeof(d_obs.data), Array{Array, 1})
    @test isequal(typeof(d_obs.geometry), GeometryIC)
    @test iszero(norm(d_obs.data[1] - d_obs.data[2]))
    @test isequal(size(d_obs), dsize)

    # set up judiVector from cell array
    data = Array{Array}(undef, nsrc)
    for j=1:nsrc
        data[j] = randn(Float32, 251, nrec)
    end
    d_obs =  judiVector(rec_geometry, data)

    @test isequal(d_obs.nsrc, nsrc)
    @test isequal(typeof(d_obs.data), Array{Array, 1})
    @test isequal(typeof(d_obs.geometry), GeometryIC)
    @test iszero(norm(d_obs.data - d_obs.data))
    @test isequal(size(d_obs), dsize)

    # contructor for in-core data container
    block = segy_read("../data/unit_test_shot_records.segy")
    d_block = judiVector(block; segy_depth_key="RecGroupElevation")
    dsize = (prod(size(block.data)), 1)

    @test isequal(d_block.nsrc, nsrc)
    @test isequal(typeof(d_block.data), Array{Array, 1})
    @test isequal(typeof(d_block.geometry), GeometryIC)
    @test isequal(size(d_block), dsize)

    # contructor for in-core data container and given geometry
    rec_geometry = Geometry(block; key="receiver", segy_depth_key="RecGroupElevation")
    d_block = judiVector(rec_geometry, block)

    @test isequal(d_block.nsrc, nsrc)
    @test isequal(typeof(d_block.data), Array{Array, 1})
    @test isequal(typeof(d_block.geometry), GeometryIC)
    @test isequal(rec_geometry, d_block.geometry)
    @test isequal(size(d_block), dsize)

    # contructor for out-of-core data container from single container
    container = segy_scan("../data/", "unit_test_shot_records", ["GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
    d_cont = judiVector(container; segy_depth_key="RecGroupElevation")

    @test isequal(d_cont.nsrc, nsrc)
    @test isequal(typeof(d_cont.data), Array{SegyIO.SeisCon, 1})
    @test isequal(typeof(d_cont.geometry), GeometryOOC)
    @test isequal(size(d_cont), dsize)

    # contructor for single out-of-core data container and given geometry
    d_cont = judiVector(rec_geometry, container)

    @test isequal(d_cont.nsrc, nsrc)
    @test isequal(typeof(d_cont.data), Array{SegyIO.SeisCon, 1})
    @test isequal(typeof(d_cont.geometry), GeometryIC)
    @test isequal(rec_geometry, d_cont.geometry)
    @test isequal(size(d_cont), dsize)

    # contructor for out-of-core data container from cell array of containers
    container_cell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        container_cell[j] = split(container, j)
    end
    d_cont =  judiVector(container_cell; segy_depth_key="RecGroupElevation")

    @test isequal(d_cont.nsrc, nsrc)
    @test isequal(typeof(d_cont.data), Array{SegyIO.SeisCon, 1})
    @test isequal(typeof(d_cont.geometry), GeometryOOC)
    @test isequal(size(d_cont), dsize)

    # contructor for out-of-core data container from cell array of containers and given geometry
    d_cont =  judiVector(rec_geometry, container_cell)

    @test isequal(d_cont.nsrc, nsrc)
    @test isequal(typeof(d_cont.data), Array{SegyIO.SeisCon, 1})
    @test isequal(typeof(d_cont.geometry), GeometryIC)
    @test isequal(rec_geometry, d_cont.geometry)
    @test isequal(size(d_cont), dsize)

#################################################### test operations ###################################################

    # conj, transpose, adjoint
    @test isequal(size(d_obs), size(conj(d_obs)))
    @test isequal(size(d_block), size(conj(d_block)))
    @test isequal(size(d_cont), size(conj(d_cont)))

    @test isequal(reverse(size(d_obs)), size(transpose(d_obs)))
    @test isequal(reverse(size(d_block)), size(transpose(d_block)))
    @test isequal(reverse(size(d_cont)), size(transpose(d_cont)))

    @test isequal(reverse(size(d_obs)), size(adjoint(d_obs)))
    @test isequal(reverse(size(d_block)), size(adjoint(d_block)))
    @test isequal(reverse(size(d_cont)), size(adjoint(d_cont)))

    # +, -, *, /
    @test iszero(norm(2*d_obs - (d_obs + d_obs)))
    @test iszero(norm(d_obs - (d_obs + d_obs)/2))

    @test iszero(norm(2*d_block - (d_block + d_block)))
    @test iszero(norm(d_block - (d_block + d_block)/2))

    @test iszero(norm(2*d_cont - (d_cont + d_cont)))    # creates in-core judiVector
    @test iszero(norm(1*d_cont - (d_cont + d_cont)/2))

    # vcat
    d_vcat = [d_block; d_block]
    @test isequal(length(d_vcat), 2*length(d_block))
    @test isequal(d_vcat.nsrc, 2*d_block.nsrc)
    @test isequal(d_vcat.geometry.xloc[1], d_block.geometry.xloc[1])

    # dot, norm, abs
    @test isapprox(norm(d_block), sqrt(dot(d_block, d_block)))
    @test isapprox(norm(d_cont), sqrt(dot(d_cont, d_cont)))
    @test isapprox(abs.(d_block.data[1]), abs(d_block).data[1]) # need to add iterate for JUDI vector

    # vector space axioms
    u = judiVector(rec_geometry, randn(Float32, ns, nrec))
    v = judiVector(rec_geometry, randn(Float32, ns, nrec))
    w = judiVector(rec_geometry, randn(Float32, ns, nrec))
    a = randn(1)[1]
    b = randn(1)[1]

    @test isapprox(u + (v + w), (u + v) + w; rtol=eps(1f0))
    @test isapprox(u + v, v + u; rtol=eps(1f0))
    #@test isapprox(u, u + 0; rtol=eps(1f0))
    @test iszero(norm(u + u*(-1)))
    @test isapprox(a .* (b .* u), (a * b) .* u; rtol=eps(1f0))
    @test isapprox(u, u .* 1; rtol=eps(1f0))
    @test isapprox(a .* (u + v), a .* u + a .* v; rtol=eps(1f0))
    @test isapprox((a + b) .* v, a .* v + b.* v; rtol=eps(1f0))

    # subsamling
    d_block_sub = d_block[1]
    @test isequal(d_block_sub.nsrc, 1)
    @test isequal(typeof(d_block_sub.geometry), GeometryIC)
    @test isequal(typeof(d_block_sub.data), Array{Array, 1})

    d_block_sub = d_block[1:2]
    @test isequal(d_block_sub.nsrc, 2)
    @test isequal(typeof(d_block_sub.geometry), GeometryIC)
    @test isequal(typeof(d_block_sub.data), Array{Array, 1})

    d_cont =  judiVector(container_cell; segy_depth_key="RecGroupElevation")
    d_cont_sub = d_cont[1]
    @test isequal(d_cont_sub.nsrc, 1)
    @test isequal(typeof(d_cont_sub.geometry), GeometryOOC)
    @test isequal(typeof(d_cont_sub.data), Array{SegyIO.SeisCon, 1})

    d_cont_sub = d_cont[1:2]
    @test isequal(d_cont_sub.nsrc, 2)
    @test isequal(typeof(d_cont_sub.geometry), GeometryOOC)
    @test isequal(typeof(d_cont_sub.data), Array{SegyIO.SeisCon, 1})

    # Conversion to SegyIO.Block
    src_geometry = Geometry(block; key="source", segy_depth_key="SourceSurfaceElevation")
    wavelet = randn(Float32, src_geometry.nt[1])
    q = judiVector(src_geometry, wavelet)
    block_out =  judiVector_to_SeisBlock(d_block, q; source_depth_key="SourceSurfaceElevation", receiver_depth_key="RecGroupElevation")

    @test isapprox(block.data, block_out.data)
    @test isapprox(get_header(block, "SourceX"), get_header(block_out, "SourceX"); rtol=1f-6)
    @test isapprox(get_header(block, "GroupX"), get_header(block_out, "GroupX"); rtol=1f-6)
    @test isapprox(get_header(block, "RecGroupElevation"), get_header(block_out, "RecGroupElevation"); rtol=1f-6)
    @test isequal(get_header(block, "ns"), get_header(block_out, "ns"))
    @test isequal(get_header(block, "dt"), get_header(block_out, "dt"))

    # Time interpolation (inplace)
    dt_orig = 2f0
    dt_new = 1f0
    nt_orig = 501
    nt_new = 1001
    d_resample = deepcopy(d_block)
    time_resample!(d_resample, dt_new; order=2)

    @test isequal(d_resample.geometry.dt[1], dt_new)
    @test isequal(d_resample.geometry.nt[1], nt_new)
    @test isequal(size(d_resample.data[1])[1], nt_new)

    time_resample!(d_resample, dt_orig; order=2)
    @test isapprox(d_resample, d_block)

    # Time interpolation (w/ deepcopy)
    d_resample = time_resample(d_block, dt_new; order=2)

    @test isequal(d_block.geometry.dt[1], dt_orig)
    @test isequal(d_block.geometry.nt[1], nt_orig)
    @test isequal(size(d_block.data[1])[1], nt_orig)
    @test isequal(d_resample.geometry.dt[1], dt_new)
    @test isequal(d_resample.geometry.nt[1], nt_new)
    @test isequal(size(d_resample.data[1])[1], nt_new)

    d_recover = time_resample(d_resample, dt_orig; order=2)
    @test isapprox(d_recover, d_block)

    # Time interpolation (linear operator)
    I = judiTimeInterpolation(d_block.geometry, dt_orig, dt_new)
    d_resample = I*d_block

    @test isequal(d_block.geometry.dt[1], dt_orig)
    @test isequal(d_block.geometry.nt[1], nt_orig)
    @test isequal(size(d_block.data[1])[1], nt_orig)
    @test isequal(d_resample.geometry.dt[1], dt_new)
    @test isequal(d_resample.geometry.nt[1], nt_new)
    @test isequal(size(d_resample.data[1])[1], nt_new)

    d_recover = transpose(I)*d_resample
    @test isapprox(d_recover, d_block)

    # scale
    a = randn(Float32, 1)[1]
    d_scale = deepcopy(d_block)

    # broadcast multiplication
    u = judiVector(rec_geometry, randn(Float32, ns, nrec))
    v = judiVector(rec_geometry, randn(Float32, ns, nrec))
    u_scale = deepcopy(u)
    v_scale = deepcopy(v)
    a = randn(1)[1]

    # broadcast identity
    u = judiVector(rec_geometry, randn(Float32, ns, nrec))
    v = judiVector(rec_geometry, randn(Float32, ns, nrec))
    u_id = deepcopy(u)
    v_id = deepcopy(v)
    broadcast!(identity, u_id, v_id)    # copy v_id into u_id

    @test isapprox(v, v_id)
    @test isapprox(v, u_id)

    # broadcast scaling + addition
    u = judiVector(rec_geometry, randn(Float32, ns, nrec))
    v = judiVector(rec_geometry, randn(Float32, ns, nrec))
    w = judiVector(rec_geometry, randn(Float32, ns, nrec))
    u_add = deepcopy(u)
    v_add = deepcopy(v)
    w_add = deepcopy(w)
    a = randn(1)[1]

    # in-place overwrite
    u = judiVector(rec_geometry, randn(Float32, ns, nrec))
    v = judiVector(rec_geometry, randn(Float32, ns, nrec))
    u_cp = deepcopy(u)
    v_cp = deepcopy(v)
    copy!(v_cp, u_cp)

    @test isapprox(u, u_cp)
    @test isapprox(u, v_cp)

    # similar
    d_zero = similar(d_block, Float32)

    @test isequal(d_zero.geometry, d_block.geometry)
    @test isequal(size(d_zero), size(d_block))
    @test iszero(d_zero.data[1])

    # retrieve out-of-core data
    d_get = get_data(d_cont)
    @test isapprox(d_block, d_get)

end


text:
#=
Project:      MIRTdev
File:         mri_geom_compress.jl
Description:  MRI Coil Compression
Author:       Rajas Gupta
Date Created: October 8th 2019
Date Update:  October 8th 2019
=#

export ir_mri_coil_compress

# using MIRT: jim, prompt, ir_mri_sensemap_sim, ir_load_brainweb_t1_256
using LinearAlgebra: svd, norm
using Plots: plot, plot!, scatter, scatter!
using Test: @inferred

"""
    (odata, σ, Vr) = ir_mri_coil_compress(idata ; ncoil)

MRI coil compression via PCA.
Given multiple MRI surface coil images (idata), use SVD/PCA
to find a smaller number of virtual coil images (odata).

In:
* `idata` `[(N) n_in]`: noisy complex images (2D or 3D) for each coil

Option:
* `ncoil` Desired # of virtual coils (default: 1)

Out:
* `odata` `[(N) ncoil]`: virtual coil images
* `σ`     `[n_in]`: singular values.
* `Vr`    `[n_in, ncoil]`: compression matrix for reducing other data.

todo: currently ignores noise correlations

"""
function ir_mri_coil_compress(idata::AbstractArray{<:Number} ; ncoil::Int = 1)

    idim = size(idata)
    n_in = idim[end]
    idata = reshape(idata, :, n_in)

    (_, σ, V) = svd(idata)

    Vr = V[:, 1:ncoil]
    odata = idata * Vr
    odata = reshape(odata, idim[1:end-1]..., ncoil)
    return (odata, σ, Vr)
end

"""
    ir_mri_coil_compress(:test)
self test with plots
"""
function ir_mri_coil_compress(test::Symbol)
    test != :test && throw(DomainError(test, "Not valid"))
    xtrue = ir_load_brainweb_t1_256()
    (nx, ny) = size(xtrue)
    ncoil = 8
    (smap, _) = ir_mri_sensemap_sim(dims = (nx, ny), dx=192/nx, ncoil=ncoil, rcoil=100)
    idata = smap .* xtrue
    sig = snr2sigma(50, idata)
    idata = idata + sig * (randn(size(idata)) + 1im * randn(size(idata)))
    nkeep = 4
#   @inferred ir_mri_coil_compress(idata, ncoil = nkeep) # todo: fails
    (odata, σ, _) = ir_mri_coil_compress(idata, ncoil = nkeep)
    for iz=1:nkeep # normalize for display
        odata[:,:,iz] = odata[:,:,iz] / maximum(abs.(odata[:,:,iz]))
    end
    p = scatter(1:nkeep, σ[1:nkeep], marker=:circle, label="")
    scatter!(p, (nkeep+1):ncoil, σ[nkeep+1:ncoil], marker=:x, label="")
    tmp = round((norm(σ[1:nkeep]) / norm(σ))^2*100, digits=2)
    plot!(p, title="percent kept $tmp")
    plot(jim(smap, "smap"), jim(idata, "idata"), jim(odata, "odata"), p)
    prompt()
    true
end

"""
    snr2sigma(db, yb)
convert SNR in dB to noise σ for complex gaussian noise
"""
function snr2sigma(db, yb::AbstractArray{<:Complex})
    10^(-db/20) * norm(yb[:]) / sqrt(length(yb)) / sqrt(2)
end

# ir_mri_coil_compress(:test)


text:
"""
    calcdoppler(gps_orb::OrbitPropagator,
                sat_orb::OrbitPropagator,
                julian_date, eop,
                obs_ecef, sig_freq)

Use `svECItoECEF` to convert satellite state vector (r and v) from ECI to ECEF.
`orb` is a list of orbit objects initialized using `SatelliteToolbox`. The first
orbit object is the source, `gps_orb`, and the second is object the signal is
reflected off of, `sat_orb`. Only two objects are used. Anything after the
second index of `orb` is ignored.

Required Arguments:

- `gps_orb::OrbitPropagator`: GPS satellite `OrbitPropagator` struct
- `sat_orb::OrbitPropagator`: second satellite `OrbitPropagator` struct
- `julian_date::Float64`: Julian date
- `eop`: Earth Orientation Parameters data which is required to convert to ECEF
    * refer to `SatelliteToolbox` module for more info
- `obs_ecef::Vector`: receiver ECEF position in meters with format `(x, y, z)`
- `sig_freq::Float64`: signal carrier frequency in Hz

Returns:

- `fd_obs::Float64`: receiver observed Doppler frequency in Hz
"""
function calcdoppler(gps_orb::OrbitPropagator,
                     sat_orb::OrbitPropagator,
                     julian_date, eop,
                     obs_ecef, sig_freq)
    # Propagate orbits to "julian_date"
    rg, vg = propagate_to_epoch!(gps_orb, julian_date)
    rs, vs = propagate_to_epoch!(sat_orb, julian_date)
    # Calculate SV acceleration vector
    ag = calculate_accel_sv(rg)
    as = calculate_accel_sv(rs)
    # Convert orbits to state vectors
    gps_teme = orbsv(julian_date, rg, vg, ag)
    sat_teme = orbsv(julian_date, rs, vs, as)
    # Transform TEME to ECEF frame
    gps_ecef = sv_eci_to_ecef(gps_teme, TEME(), ITRF(), julian_date, eop)
    sat_teme = sv_eci_to_ecef(sat_teme, TEME(), ITRF(), julian_date, eop)
    # Get velocities and positions
    rg = gps_ecef.r
    vg = gps_ecef.v
    rs = sat_ecef.r
    vs = sat_ecef.v
    # Calculate observed radial velocities
    # Radial velocity of satellite relative to observer
    v_obs = transpose((rs-obs_ecef)/norm(rs-obs_ecef))*vs
    # Radial velocity of satellite relative to GPS satellite
    v_g2s = transpose((rs-rg)/norm(rs-rg))*vs
    # Radial velocity of GPS satellite relative to satellite
    v_s2g = transpose((rg-rs)/norm(rg-rs))*vg
    # Calculate observed Doppler frequency by observer
    fd_obs = -sig_freq*(v_obs + v_g2s + v_s2g)/c  # Hz
    return fd_obs
end

"""
    calcdoppler(r::Vector, v::Vector, obs_ecef, sig_freq)

Calculates the Doppler frequency for the direct signal case.

Required Arguments:

- `rg::Vector`: first object ECEF position vector in meters with format `(x, y, z)`
- `vg::Vector`: first object ECEF velocity vector in meters with format `(v_x, v_y, v_z)`
- `rs::Vector`: second object ECEF position vector in meters with format `(x, y, z)`
- `vs::Vector`: second object ECEF velocity vector in meters with format `(v_x, v_y, v_z)`
- `julian_date::Float64`: Julian date
- `obs_ecef::Vector`: receiver ECEF position in meters with format `(x, y, z)`
- `sig_freq::Float64`: signal carrier frequency in Hz

Returns:

- `fd_obs::Float64`: receiver observed Doppler frequency in Hz
"""
function calcdoppler(rg::Vector, vg::Vector, rs::Vector, vs::Vector,
                     obs_ecef, sig_freq)
    # Calculate observed radial velocities
    # Radial velocity of GPS satellite relative to observer
    v_obs = transpose((rs-obs_ecef)/norm(rs-obs_ecef))*vs
    # Radial velocity of satellite relative to GPS satellite
    v_g2s = transpose((rs-rg)/norm(rs-rg))*vs
    # Radial velocity of GPS satellite relative to satellite
    v_s2g = transpose((rg-rs)/norm(rg-rs))*vg
    # Calculate observed Doppler frequency by observer
    fd_obs = -sig_freq*(v_obs + v_g2s + v_s2g)/c  # Hz
    return fd_obs
end

"""
    calcdoppler(orb::OrbitPropagator, julian_date, eop,
                obs_ecef, sig_freq)

Calculates the Doppler frequency for the direct signal case.

Required Arguments:

- `orb::OrbitPropagator`: satellite `OrbitPropagator` struct
- `julian_date::Float64`: Julian date
- `eop`: Earth Orientation Parameters data which is required to convert to ECEF
    * refer to `SatelliteToolbox` module for more info
- `obs_ecef::Vector`: receiver ECEF position in meters with format `(x, y, z)`
- `sig_freq::Float64`: signal carrier frequency in Hz

Returns:

- `fd_obs::Float64`: receiver observed Doppler frequency in Hz
"""
function calcdoppler(orb::OrbitPropagator, julian_date, eop,
                     obs_ecef, sig_freq)
    # Propagate orbits to "julian_date"
    r, v = propagate_to_epoch!(orb, julian_date)
    aᵢ = calculate_accel_sv(r)
    # Convert orbits to state vectors
    orb_teme = orbsv(julian_date, r, v, aᵢ)
    # Transform TEME to ECEF frame
    orb_ecef = sv_eci_to_ecef(orb_teme, TEME(), ITRF(), julian_date, eop)
    # Get velocities and positions
    r = orb_ecef.r
    v = orb_ecef.v
    # Direct signal case
    v_obs = transpose(v)*(r-obs_ecef) / norm(r-obs_ecef)
    # Calculate observed Doppler frequency by observer
    fd_obs = -sig_freq*v_obs/c  # Hz
    return fd_obs
end

"""
    calcdoppler(r::Vector, v::Vector, obs_ecef, sig_freq)

Calculates the Doppler frequency for the direct signal case.

Required Arguments:

- `r`: object ECEF position vector in meters with format `(x, y, z)`
- `v`: object ECEF velocity vector in m/s with format `(v_x, v_y, v_z)`
- `obs_ecef`: receiver ECEF position in meters with format `(x, y, z)`
- `sig_freq`: signal carrier frequency in Hz

Returns:

- `fd_obs::Float64`: receiver observed Doppler frequency in Hz
"""
function calcdoppler(r, v, obs_ecef, sig_freq)
    # Direct signal case
    v_obs = transpose(v)*(r-obs_ecef) / norm(r-obs_ecef)
    # Calculate observed Doppler frequency by observer
    fd_obs = -sig_freq*v_obs/c  # Hz
    return fd_obs
end

"""
    calcdopplerrate(r::Vector, v::Vector, a::Vector, obs_ecef, sig_freq)

Calculates the Doppler frequency for the direct signal case.

Required Arguments:

- `r`: object ECEF position vector in meters with format `(x, y, z)`
- `v`: object ECEF velocity vector in m/s with format `(v_x, v_y, v_z)`
- `a`: object ECEF velocity vector in m/s² with format `(a_x, a_y, a_z)`
- `obs_ecef`: receiver ECEF position in meters with format `(x, y, z)`
- `sig_freq`: signal carrier frequency in Hz

Returns:

- `fd_rate_obs::Float64`: receiver observed Doppler frequency in Hz
"""
function calcdopplerrate(r, v, a, obs_ecef, sig_freq)
    # Direct signal case
    a_obs = transpose(a)*(r-obs_ecef) / norm(r-obs_ecef)  +
            norm(v)^2 / norm(r-obs_ecef) -
            (transpose(v)*(r-obs_ecef) / norm(r-obs_ecef))^2 / norm(r-obs_ecef)
    # Calculate observed Doppler frequency rate by observer
    fd_rate_obs = -sig_freq*a_obs/c  # Hz/s
    return fd_rate_obs
end


text:
"""
    use_sigdigits(str)
Evaluate input string while respecting significant digits.

# Examples
```julia-repl
julia> use_sigdigits("1.801 + 3.5")
5.3
```
"""
function use_sigdigits(str::String)
    ex = Meta.parse(str)
    return eval_sigdigits(ex)
end

"""
    use_sigdigits(expr)
Evaluate input expression while respecting significant digits.

# Examples
```julia-repl
julia> ex = Meta.parse("1.801 + 3.5")
julia> use_sigdigits(ex)
5.3
```
"""
function use_sigdigits(ex::Expr)
    return eval_sigdigits(ex)
end

"""
Output text in specified color.
"""
function print_color(text, r::Int, g::Int, b::Int)
    print("\e[1m\e[38;2;$r;$g;$b;249m",text)
end

"""
    sig_start()
Start significant digit evaluation mode with default colors.
"""
function sig_start()
    sig_start(100, 200, 255)
end

"""
    sig_start(r::Int, g::Int, b::Int)
Start significant digit evaluation mode with specified colors.
"""
function sig_start(r::Int, g::Int, b::Int)
    println("Starting input mode...")
    while true
        print_color("sig>", r, g, b)
        print_color(" ", 240, 240, 240)
        input = chomp(readline())
        if strip(input) == "exit"; break
        elseif strip(input) == "help"
            print("In input mode, any valid math expression ")
            println("you input will be handled along with its significant digits.")
        else
            try
                ex = Meta.parse(input)
                println(eval_sigdigits(ex))
            catch e
                if isa(e, ErrorException)
                    print_color("Error: invalid input!", 255, 100, 100)
                    print_color(" \n", 240, 240, 240)
                else
                    stacktrace(catch_backtrace())
                end
            end
        end
    end
end

"""
Calculate the minimum significant digits in the arguments of an expression.
"""
function find_min_sig(args, minset::Tuple, ismult)#hasmult, hasadd)
    new_minsig = minset[1]
    new_mindec = minset[2]
    tmpsets = Vector{Tuple}()
    # iterate through expression args and handle nested expressions and minimums
    for el in args
        if typeof(el) == Expr
            tmpset = calc_sig(el, minset)
            push!(tmpsets, tmpset)
        elseif typeof(el) != Symbol
            tmpsig = count_sigdigits(el)
            tmpdec = count_decimalplaces(el)
            if ismult && tmpsig < new_minsig; new_minsig = tmpsig; end
            if !ismult && tmpdec < new_mindec; new_mindec = tmpdec; end
        end
    end
    # interate through nested expression outputs to find minimum significant digits and decimal places
    for set in tmpsets
        if set[1] < new_minsig; new_minsig = set[1] end
        if set[2] < new_mindec; new_mindec = set[2] end
    end
    if new_minsig != Inf; new_minsig = floor(Int, new_minsig) end
    if new_mindec != Inf; new_mindec = floor(Int, new_mindec) end
    println("new_minset: ", (new_minsig, new_mindec))
    return (new_minsig, new_mindec)
end

"""
Calculate the minimum significant digits in an expression.
"""
function calc_sig(ex, minset::Tuple)
    operator = ex.args[1]
    ismult = false
    if operator == :* || operator == :/; ismult = true
    elseif operator == :+ || operator == :-; ismult = false
    end
    println("minset: ", minset)
    return find_min_sig(ex.args, minset, ismult)
end

"""
    eval_digdigits(ex:Expr)
Evaluate an expression while respecting significant digits.
"""
function eval_sigdigits(ex)
    minset = calc_sig(ex, (typemax(Float32), typemax(Float32)))
    ans = eval(ex)
    println(minset)
    if minset[1] != Inf; ans = round(ans, sigdigits=minset[1]) end
    nondec = count_digits(ans) - count_decimalplaces(ans)
    println("DIGITS: ", count_digits(ans), " ", count_decimalplaces(ans), " ", minset)
    if minset[2] != Inf; ans = round(ans, digits=nondec+minset[2]) end
    return ans
end

"""
    str_replace(str::String, first::Int, last::Int, new::String)
Return new string with replaced substring.
"""
function str_replace(str::String, first::Int, last::Int, new::String)
    front = str[1: first]
    back = str[last+1: end]
    return front * new * back;
end

"""
    count_digits(val)
Enumerate the number of digits in a number.
"""
function count_digits(val)
    if val == 0 return 0 end
    strval = string(decimal(val))
    len = length(strval)
    if val < 0; len -= 1 end
    if findfirst(".", strval) != nothing; len -= 1 end
    return len
end

"""
    count_sigdigits(val)
Enumerate the number of significant digits in a number.
"""
function count_sigdigits(val)
    if val == 0 return 0 end
    num = string(decimal(val))
    len = length(num)
    if val < 0; len -= 1; num = num[2:end] end
    sigdigits = 0
    unwrapped_zeros = 0 # resets every time zeros are wrapped
    prevnonzero = false
    for i in 1:len
        placeval = num[i]
        if placeval == '.'; continue
        elseif placeval == '0'; prevnonzero && (unwrapped_zeros += 1; true)
        elseif unwrapped_zeros > 0
            sigdigits += unwrapped_zeros + 1
            unwrapped_zeros = 0
            prevnonzero = true
        else
            sigdigits += 1
            prevnonzero = true
        end
    end
    return sigdigits
end

"""
    count_decimalplaces(val)
Enumerate the number of decimal places a float contains.
"""
function count_decimalplaces(val)
    num = decimal(val)
    decimalslug = num - trunc(num)
    if decimalslug == 0; return 0; end
    decimals = 0
    while decimalslug > 0
        decimals += 1
        placeval = trunc(decimalslug, digits=decimals)
        decimalslug -= placeval
    end
    return decimals
end

"""
    count_decimalplaces(val)
Enumerate the number of nonzero decimal places a float contains.
"""
function count_nonzero_decimalplaces(val)
    num = decimal(val)
    decimalslug = num - trunc(num)
    if decimalslug == 0; return 0; end
    decimals = 0
    nonzeros = 0
    while decimalslug > 0
        decimals += 1
        placeval = trunc(decimalslug, digits=decimals)
        decimalslug -= placeval
        if placeval > 0; nonzeros += 1; end
    end
    return nonzeros
end


text:
## Exponential moving average
function EWMA(x)
    ewmaVar = zeros(length(x),1)
    ewmaVar[1] = var(x)
    for i in 2:length(x)
         ewmaVar[i] = 0.94*ewmaVar[i-1]+ (1.0 - 0.94)*x[i-1]^2.0
    end
    ewmaVol = sqrt.(ewmaVar);
    return ewmaVol
end

## OLS estimate
function OLSestimator(y,x)
    estimate = (transpose(x)*x) \ (transpose(x)*y)
    return estimate
end

## Function to generate polynoms
function delta3j(js)
    alphaj =  collect(linspace(0.05, 0.95, js))
    delta0 = collect(linspace(1.0, 1.0, length(alphaj)))
    delta1 = 2.0*(alphaj - 0.5)
    delta2 = (2.0^2.0)*((alphaj - 0.5).^2.0)
    delta3 = (2.0^3.0)*((alphaj - 0.5).^3.0)
    delta4 = (2.0^4.0)*((alphaj - 0.5).^4.0)
    delta5 = (2.0^5.0)*((alphaj - 0.5).^5.0)
    return [delta0 delta1 delta2 delta3 delta4 delta5]
end

## Fucntion to obtain initial estimates
function getinitlogit(rets, js,p1,p2)
    n = length(rets)
    x2 = log.(1 + abs.(rets))[1:(n-1)];
    alphaj = collect(linspace(0.05, 0.95, js))
    coefdf = zeros(js, 3)
    estlog = 0.0
    for j in collect(1:(js))
        indx1 = (rets .< (EWMA(rets).*quantile.(Normal(0,1), alphaj[j])) ).*1.0;
        data1 = DataFrame([indx1[2:n], indx1[1:(n-1)], x2], [:y, :x1, :x2]);

        estlog = glm(@formula(y ~ x1 + x2), data1, Binomial(), LogitLink());
        coefdf[j,:] = coef(estlog)
    end

    ## Estimation of kappas from deltas
    delta = delta3j(js)
    kapa1 = OLSestimator(coefdf[:, 2], delta[:,1:(p1+1)]);
    kapa2 = OLSestimator(coefdf[:, 3], delta[:,1:(p2+1)]);
    kapainit = [coefdf[:, 1]; kapa1; kapa2]
    return kapainit
end

## LL function
function loglikeFinalstd(x, param0,p1,p2)
    ##
    rets = x
    n = length(rets)
    m = zeros(n)
    theta = zeros(n)
    n_alpha = length(param0) - (p1+p2+2)
    prob = zeros(n, n_alpha)

    alphaj = collect(linspace(0.05, 0.95, n_alpha))

    cuts=zeros(n,n_alpha);
    for i in collect(1:n_alpha)
        cuts[:,i]=EWMA(rets)*quantile.(Normal(0,1), alphaj[i])
    end
    ##
    powers = [p1; p2]
    ## looping over time=j, and quantiles=i
    for j in collect(1:n)
        for i in collect(1:n_alpha)
            ##
            theta[j] = param0[1 + (i-1)*1]
            poly_1 = param0[1 + n_alpha]
            poly_2 = param0[2 + n_alpha + powers[1]]
            ## first polynoms in delta1
            for l in collect(1:powers[1])
                poly_1 = poly_1 + param0[1 + n_alpha + l] * (2.0 * alphaj[i] - 1.0)^l
            end
            ## second polynoms in delta2
            for l in collect(1:powers[2])
                poly_2 = poly_2 + param0[2 + n_alpha + powers[1]+l] * (2.0 * alphaj[i] - 1.0)^l
            end
            ##
            if j==1
                theta[j]=theta[j] + poly_1*mean(rets.<cuts[:,i]) + poly_2*mean(log.(1.0+abs.(rets)))
            else
                theta[j]=theta[j] + poly_1*(rets[j-1]<cuts[j-1,i]) + poly_2*log.(1.0+abs.(rets[j-1]))
            end
            ##
            prob[j, i] = 1.0 / (1.0 + exp.(-theta[j]))
        end
        m[j] = (rets[j] <= cuts[j,1]) * log.(prob[j, 1])
        for i in collect(2:n_alpha)
            if (prob[j, i] - prob[j, i-1]) <= 0.0
                prob[j, i] = prob[j, i-1] + 0.00001
            end

            m[j] = m[j] + ((rets[j] <= cuts[j,i]) * (rets[j] > cuts[j,i-1])) * log.(prob[j, i] - prob[j, i-1])
        end
            if (1.0 - prob[j, n_alpha]) <= 0.0
                prob[j, n_alpha] = 0.99999
            end
        m[j] = m[j] + (rets[j] > cuts[j,n_alpha]) * log.(1.0 - prob[j, n_alpha])
    end
    ##
    return m
end

function loglikeFinal(x, param0,p1,p2)

    m = loglikeFinalstd(x, param0,p1,p2)

    return -sum(m)
end

# Estimate parameters
function OrderedLogitparameters(x,q,p1,p2)

    kapainit = getinitlogit(x,q,p1,p2)
    estone = optimize(kapainit -> loglikeFinal(x, kapainit,p1,p2), kapainit, NelderMead(), Optim.Options(g_tol=1e-5, iterations=100))

    par = estone.minimizer

    return par
end

# Parameters and inference
function OrderedLogit(x,q,p1,p2)

    kapainit = getinitlogit(x,q,p1,p2)
    estone = optimize(kapainit -> loglikeFinal(x, kapainit,p1,p2), kapainit, NelderMead(), Optim.Options(g_tol=1e-5, iterations=100))

    par = estone.minimizer

    n=length(x)
    b=par'
    k=length(b)
    frac = 0.00001;
    h=frac*b

    e = eye(k)

    DfDp=zeros(n,k)
    for i in 1:k
        DfDp[:,i] = (loglikeFinalstd(x, b'+h'.*e[i,:],p1,p2)-loglikeFinalstd(x, b'-h'.*e[i,:],p1,p2))./h[i]/2.0
    end

    DfDp2=zeros(k,k);
    for i in 1:k
        for j in 1:i
            DfDp2[i,j] = mean(DfDp[:,i].*DfDp[:,j])
            DfDp2[j,i] = DfDp2[i,j];
        end
    end

    Df2Dp=zeros(k,k);
    for i in 1:k
    for j in 1:i
        if i==j;
            Df2Dp[i,i] = mean((loglikeFinalstd(x, b'+h'.*e[i,:],p1,p2)-2.0*loglikeFinalstd(x, b,p1,p2)+loglikeFinalstd(x, b'-h'.*e[i,:],p1,p2))./(h[i]^2.0));
        else
            Df2Dp[i,j] = mean((loglikeFinalstd(x, b'+h'.*(e[i,:]+e[j,:]),p1,p2)-loglikeFinalstd(x, b'+h'.*(e[i,:]-e[j,:]),p1,p2)-loglikeFinalstd(x, b'+h'.*(e[j,:]-e[i,:]),p1,p2)+loglikeFinalstd(x, b'-h'.*(e[i,:]+e[j,:]),p1,p2))./(4.0*h[i]*h[j]));
            Df2Dp[j,i] = Df2Dp[i,j];
        end
    end
    end
    Df2Dp=-Df2Dp;
    J = inv(Df2Dp)

    sterr=sqrt.(diag(J*DfDp2*J)/n)
    t=b'./sterr
    LL=estone.minimum

    AIC = -2LL + 2k + 2k*(k+1)/(n-k-1)
    BIC = -2LL + k*log(n)

    return [par,sterr,t,LL,AIC,BIC]
end

function forecastProbs(x, y, q,p1,p2)

    kapainit = getinitlogit(x,q,p1,p2)

    estone = optimize(kapainit -> loglikeFinal(x, kapainit,p1,p2), kapainit, NelderMead(), Optim.Options(g_tol=1e-5, iterations=100))
    param0 = estone.minimizer

    rets0 = x
    rets = y
    n = length(rets)
    m = zeros(n)
    theta = zeros(n)
    n_alpha = length(param0) - (p1+p2+2)
    prob = zeros(n, n_alpha)
    alphaj = collect(linspace(0.05, 0.95, n_alpha))
    cuts=EWMA(rets)[length(rets)].*quantile.(Normal(0,1), alphaj)
    powers = [p1; p2]
    wrong=zeros(n)
    ## looping over time=j, and quantiles=i
    for j in collect(2:n)
        for i in collect(1:n_alpha)
            ##
            theta[j] = param0[1 + (i-1)*1]
            poly_1 = param0[1 + n_alpha]
            poly_2 = param0[2 + n_alpha + powers[1]]
            ## first polynoms in delta1
            for l in collect(1:powers[1])
                poly_1 = poly_1 + param0[1 + n_alpha + l] * (2.0 * alphaj[i] - 1.0)^l
            end
            ## second polynoms in delta2
            for l in collect(1:powers[2])
                poly_2 = poly_2 + param0[2 + n_alpha + powers[1]+l] * (2.0 * alphaj[i] - 1.0)^l
            end

            theta[j]=theta[j] + poly_1*(rets[j-1]<cuts[i]) + poly_2*log.(1.0+abs.(rets[j-1]))

            prob[j, i] = 1.0 / (1.0 + exp.(-theta[j]))
        end
        for i in collect(2:n_alpha)
            if (prob[j, i] - prob[j, i-1]) <= 0.0
                    prob[j, i] = prob[j, i-1] + 0.00001
                    #wrong[j] = 1
            end
        end
    end

    return prob[2:n,:]
end

## ------------------------------------------
## Statistical evaluation of forecasts
## ------------------------------------------

function Interpol(xs,ys,x)
    n = length(xs)
    dxs = xs[2:n]-xs[1:n-1]
    dys = ys[2:n]-ys[1:n-1]
    Ds = dys./dxs
    ms = [Ds[1];(Ds[1:n-2]+Ds[2:n-1])/2;Ds[n-1]]
    as = ms[1:n-1]./Ds[1:n-1]
    bs = ms[2:n]./Ds[1:n-1]
    cir = (as.^2)+(bs.^2)
    ex = cir.>9
    tau = 3./sqrt.(cir).*ex + 1-ex
    ms = ms.*[tau;1].*[1;tau]
    if findfirst(xs.>x) == 0
        ind_u = 2
    elseif findfirst(xs.>x) == 1
        ind_u = 2
    else
        ind_u = findfirst(xs.>x)
    end
    x_u = xs[ind_u]
    x_l = xs[ind_u-1]
    h = x_u-x_l
    t = (x-x_l)./h
    fx = ys[ind_u-1]*(2*t^3-3*t^2+1) + h*ms[ind_u-1]*(t^3-2*t^2+t) + ys[ind_u]*(-2*t^3+3*t^2) + h*ms[ind_u]*(t^3-t^2)

    return fx
end

function Berkowitz(pit)
    Zb = zeros(length(pit))
    pit[pit.>=1] =0.9999999
    pit[pit.<=0] =0.0000001
    for t in collect(1:length(pit))
        Zb[t] = norminvcdf(abs.(pit[t]))
    end
    naux = length(Zb)-1
    mu_hat = sum(Zb)/naux

    Xaux = Zb[1:(length(Zb)-1)]-mu_hat
    Yaux = Zb[2:length(Zb)]-mu_hat

    rho_hat = Xaux'Yaux/(Xaux'Xaux)
    eps_hat = Yaux-rho_hat*Xaux
    s2_hat = sum(eps_hat.^2)/naux

    LLu = -1/2*naux*log(s2_hat) - 1/2/s2_hat*eps_hat'eps_hat
    LLr = - 1/2*Yaux'Yaux
    LR = -2*(LLr-LLu)
    return [LR,ccdf(Chisq(3), LR)]
end

function GonzalezRivera(pit,alpha_c,alpha_main,k_main)
    # alpha_main    -- alpha-contour for k-aggregated test
    # k_main        -- lag for alpha-aggregated test

    n_a = length(alpha_c)   # collection of alpha-contours
    K_c=collect(1:n_a)
    n_k = length(K_c)                   # collection of lags for alpha-aggregated test
    Lambda_a = zeros(n_k,n_k)
    L_a =  zeros(n_k,1)
    report1=zeros(n_a,2,n_k)
    report2=zeros(1,2)
    report3=zeros(1,2)
    for k in 1:n_k
        alpha_contour_hat = zeros(length(alpha_c),1)
        PITk = [pit[1:(length(pit)-k)] pit[(1+k):length(pit)]]
        Tk = size(PITk)[1]
        Ind = zeros(Tk,n_a)
        Omega_k = zeros(n_a,n_a)

        for i in 1:n_a
            Indi = 1*(PITk.<sqrt(alpha_c[i]))
            Ind[:,i] = Indi[:,1].*Indi[:,2];
            alpha_contour_hat[i] = sum(Ind[:,i])/Tk;

            for j in 1:(i-1)
                    Omega_k[i,j] = alpha_c[j].*(1-alpha_c[i])+2*alpha_c[j]*alpha_c[i]^0.5.*(1-alpha_c[i]^0.5)
            end
            Omega_k = Omega_k + Omega_k'
            Omega_k[i,i] = alpha_c[i].*(1-alpha_c[i])+2*alpha_c[i]^1.5.*(1-alpha_c[i]^0.5)
        end

        for j in 1:n_k
            Lambda_a[j,k] = 4*alpha_main^1.5*(1-alpha_main^0.5)
        end
        Lambda_a[k,k] = alpha_main.*(1-alpha_main)+2*alpha_main^1.5.*(1-alpha_main^0.5)
        C_k = sqrt(Tk)*(alpha_contour_hat-alpha_c)
        t_c = C_k./sqrt.(alpha_c.*(1-alpha_c)+2*alpha_c.^1.5.*(1-alpha_c.^0.5))
        C_k = C_k'inv(Omega_k)*C_k

        Indi = 1*(PITk.<sqrt(alpha_main))
        Ind = Indi[:,1].*Indi[:,2]
        alpha_main_hat = sum(Ind)/Tk

        L_a[k] = sqrt(Tk)*(alpha_main_hat-alpha_main)

        # report individual k stats and p-vals
        report1[:,:,k]=[t_c ccdf.(Chisq(1),(t_c.^2))]

        # report alpha-aggregated
        if k == k_main
            report2=[C_k ccdf.(Chisq(n_a), C_k)]
        end

    end
    # report lag-aggregated
    L_a = L_a'inv(Lambda_a)*L_a;

    report3=[L_a ccdf.(Chisq(n_k), L_a)]

    return [report1,report2,report3]
end

## NUMERICAL INTEGRATION ROUTINES

function ComputeInt(r,rets,cuts,ext_probs,n)

    a = 2*minimum(rets)
    b = 2*maximum(rets)
    x = cos.(pi*(2*collect(1:1:n)-1)./2./n)
    S = 0

    for i in 1:n
        y=(x[i]+1)*(b-a)/2+a
        CDFy = Interpol([2*minimum(rets);cuts;2*maximum(rets)],ext_probs,y)
        S = S + ((CDFy-1*(y.>=r)).^2)*sqrt(1-x[i]^2)
    end

    return S*pi*(b-a)/2/n
end

max_stars_repo_name: barunik/DistributionalForecasts

max_issues_repo_name: barunik/DistributionalForecasts

max_forks_repo_name: barunik/DistributionalForecasts


text:
module InductionlessMHD

using Gridap

function dimensionless_a(X,Y,Re,N,B,dΩ, coord)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( (1/Re)*inner(∇(u),∇(v_u)) - p*(∇⋅v_u) + (N∘(coord))*((j×(B∘(coord)))⋅v_u)*(-1.0) +
     (∇⋅u)*v_p +
     j⋅v_j - φ*(∇⋅v_j) - (u×(B∘(coord)))⋅v_j +
     (∇⋅j)*v_φ ) * dΩ
end

function dimensionless_a(X,Y,Re,N::Float64,B,dΩ, coord)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( (1/Re)*inner(∇(u),∇(v_u)) - p*(∇⋅v_u) - N*((j×(B∘(coord)))⋅v_u) +
     (∇⋅u)*v_p +
     j⋅v_j - φ*(∇⋅v_j) - (u×(B∘(coord)))⋅v_j +
     (∇⋅j)*v_φ ) * dΩ
end

function dimensionless_a(X,Y,Re,N,B::VectorValue,dΩ, coord)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( (1/Re)*inner(∇(u),∇(v_u)) - p*(∇⋅v_u) + (N∘(coord))*((j×B)⋅v_u)*(-1.0) +
     (∇⋅u)*v_p +
     j⋅v_j - φ*(∇⋅v_j) - (u×B)⋅v_j +
     (∇⋅j)*v_φ ) * dΩ
end

function dimensionless_a(X,Y,Re,N::Float64,B::VectorValue,dΩ, coord)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( (1/Re)*inner(∇(u),∇(v_u)) - p*(∇⋅v_u) - N*(j×B)⋅v_u +
     (∇⋅u)*v_p +
     j⋅v_j - φ*(∇⋅v_j) - (u×B)⋅v_j +
     (∇⋅j)*v_φ ) * dΩ
end

function dimensionless_a(X,Y,Re,N::Float64,B::VectorValue,dΩ)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( (1/Re)*inner(∇(u),∇(v_u)) - p*(∇⋅v_u) - N*(j×B)⋅v_u +
     (∇⋅u)*v_p +
     j⋅v_j - φ*(∇⋅v_j) - (u×B)⋅v_j +
     (∇⋅j)*v_φ ) * dΩ
end

function dimensionless_l(Y,f_u,dΩ)
  v_u, v_p, v_j, v_φ = Y

  ∫( v_u⋅f_u ) * dΩ
end

conv(u,∇u) = (∇u')⋅u
dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)

function convective_term(X,Y,dΩ)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y

  ∫( v_u⋅conv(u,∇(u)) ) * dΩ
end

function dconvective_term(X,dX,Y,dΩ)
  u  , p  , j  , φ   = X
  du , dp , dj , dφ  = dX
  v_u, v_p, v_j, v_φ = Y

  ∫( v_u⋅dconv(du,∇(du),u,∇(u)) ) * dΩ
end

function dimensionless_residual(X,Y,Re,N,B,dΩ, coord)
  dimensionless_a(X,Y,Re,N,B,dΩ, coord) +
  convective_term(X,Y,dΩ)
end

function dimensionless_residual(X,Y,Re,N,B,f_u,dΩ, coord)
  dimensionless_a(X,Y,Re,N,B,dΩ, coord) +
  convective_term(X,Y,dΩ) -
  dimensionless_l(Y,f_u,dΩ)
end

function dimensionless_jacobian(X,dX,Y,Re,N,B,dΩ, coord)
  dimensionless_a(dX,Y,Re,N,B,dΩ, coord) +
  dconvective_term(X,dX,Y,dΩ)
end

function dimensionless_residual(X,Y,Re,N,B,dΩ)
  dimensionless_a(X,Y,Re,N,B,dΩ) +
  convective_term(X,Y,dΩ)
end

function dimensionless_residual(X,Y,Re,N,B,f_u,dΩ::Measure)
  dimensionless_a(X,Y,Re,N,B,dΩ) +
  convective_term(X,Y,dΩ) -
  dimensionless_l(Y,f_u,dΩ)
end

function dimensionless_jacobian(X,dX,Y,Re,N,B,dΩ)
  dimensionless_a(dX,Y,Re,N,B,dΩ) +
  dconvective_term(X,dX,Y,dΩ)
end

function dimensionless_darcy_l_Γ(Y,n,g_φ,dΓ)
  v_u, v_p, v_j, v_φ = Y

  ∫( -(v_j⋅n)*g_φ ) * dΓ
end

function dimensionless_conducting_wall(X,Y,n,c_w,dΓ;α=1.0)
  u  , p  , j  , φ   = X
  v_u, v_p, v_j, v_φ = Y
  ∫( α*((v_j⋅n)*(j⋅n) + c_w * (v_j⋅n)*(n⋅(∇(j)⋅n))) ) * dΓ
end

end # module


text:
using LinearAlgebra
using FermiCG
using Printf
using Test

@testset "cmf" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,0,1]))
    push!(atoms,Atom(3,"H",[0,0,2]))
    push!(atoms,Atom(4,"H",[0,0,3]))
    push!(atoms,Atom(5,"H",[0,0,4]))
    push!(atoms,Atom(6,"H",[0,0,5]))
    #basis = "6-31g"
    basis = "sto-3g"

    mol     = Molecule(0,1,atoms,basis)
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1a_fci, d1b_fci,d2_fci = FermiCG.pyscf_fci(ints,3,3)
    # @printf(" FCI Energy: %12.8f\n", e_fci)

    function compute_npairs(d2)
        npairs = 0
        for i in 1:size(d2)[1]
            for j in 1:size(d2)[1]
                npairs += d2[i,i,j,j]
            end
        end
        println(" NPairs:", round(npairs,digits=3))
    end

    compute_npairs(d2_fci)

    FermiCG.pyscf_write_molden(mol,mf.mo_coeff,filename="scf.molden")

    # Cl = FermiCG.localize(mf.mo_coeff,"boys", mf)
    # FermiCG.pyscf_write_molden(mol,Cl,filename="boys.molden")
    C = mf.mo_coeff
    rdm_mf = C[:,1:2] * C[:,1:2]'
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    clusters    = [(1:2),(3:4),(5:6)]
    #clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(1,1),(1,1),(1,1)]

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))
    rdm1a = rdm_mf
    rdm1b = rdm_mf

    # test in core method
    #f1 = FermiCG.cmf_ci_iteration(ints, clusters, rdm1a, rdm1b, init_fspace, verbose=1)
    #@test isapprox(f1[1], -2.876651063218, atol=1e-10)

    # test on the fly integral method
    #f2 = FermiCG.cmf_ci_iteration(mol, Cl, rdm1a, rdm1b, clusters, init_fspace, verbose=1)
    #@test isapprox(f2[1], -2.876651063218, atol=1e-10)

    f1 = FermiCG.cmf_ci(ints, clusters, init_fspace, rdm1a, rdm1b, verbose=1,sequential=false)
    #f2 = FermiCG.cmf_ci(mol, Cl, clusters, init_fspace, rdm1a, verbose=0)
    #@test isapprox(f1[1], f2[1], atol=1e-10)
    @test isapprox(f1[1], -2.97293813654926351, atol=1e-10)

    e_cmf, U = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1, rdm1, verbose=0, gconv=1e-6, method="cg",sequential=true)
    #e_cmf, U = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1, verbose=0, gconv=1e-6, max_iter_oo=4)
    @test isapprox(e_cmf, -3.205983033016, atol=1e-10)

    ## This is so slow, so only test a couple iterations
    ##e_cmf, U = FermiCG.cmf_oo(mol, Cl, clusters, init_fspace, rdm1, verbose=0, gconv=1e-6)
    ##@test isapprox(e_cmf, -3.205983033016, atol=1e-10)
    #e_cmf, U = FermiCG.cmf_oo(mol, Cl, clusters, init_fspace, rdm1, verbose=0, gconv=1e-6, max_iter_oo=3)
    #@test isapprox(e_cmf, -3.2036779028934497, atol=1e-10)
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
__precompile__()

module Econometrics

# Utilities
include("Utilities/stnorm.jl")
include("Utilities/prettyprint.jl")
include("Utilities/sortbyc.jl")
include("Utilities/dstats.jl")
include("Utilities/lag.jl")
include("Utilities/lags.jl")
include("Utilities/vech.jl")
include("Utilities/PrintDivider.jl")
include("Utilities/PrintEstimationResults.jl")
# linear regression
include("LinearRegression/lsfit.jl")
include("LinearRegression/ols.jl")
include("LinearRegression/tsls.jl")
include("LinearRegression/TestStatistics.jl")
# Bayesian
include("Bayesian/mcmc.jl")
# nonparametrics
include("NP/npreg.jl")
include("NP/npdensity.jl")
include("NP/kernelweights.jl")
include("NP/NeweyWest.jl")
#include("NN/TrainNet.jl")
#include("NN/AnalyzeNet.jl")
# optimization
include("Optimization/samin.jl")
include("Optimization/fminunc.jl")
include("Optimization/fmincon.jl")
# MLE
include("ML/mle.jl")
include("ML/mleresults.jl")
include("ML/Likelihoods/logit.jl")
include("ML/Likelihoods/poisson.jl")
include("ML/Likelihoods/normal.jl")
# GMM
include("GMM/gmm.jl")
include("GMM/gmmresults.jl")
export stnorm, prettyprint, sortbyc, dstats, lag, lags, vech
export PrintDivider, PrintEstimationResults
export lsfit, ols, tsls, TestStatistics, NeweyWest
export mcmc
export npreg, kernelweights
export npdensity
#export TrainNet, AnalyzeNet
export samin, fminunc, fmincon
export mle, mleresults, logit, poisson, normal
export gmm, gmmresults

end


text:
@with_kw struct LavaWorld <: MDP{GWPos, Symbol}
    gridworld::SimpleGridWorld # Underlying grdiworld
    lava::Union{Symbol, Array{GWPos}} = [GWPos(1,1)] # Define the lava or select :random
    num_lava_tiles::Int = 1 # Number of squares of lava to use
    possible_lava_squares::Union{Symbol, Array{GWPos}} = :all # Set of possible lava squares. :all means all of them can be used
    goal::Union{Symbol, GWPos} = GWPos(7,5) # can be :random
    randomize_lava_on_reset::Bool = false # Whether or not to randomize the lava on each reset
    lava_penalty::Float64 = -1. # The reward for landing in lava
    goal_reward::Float64 = 1. # The reward for landing in the goal state
    observation_type = :all # if :all, we get the full grid over 3 channels
    rng::AbstractRNG = Random.GLOBAL_RNG
    initial_state = :random
end

function LavaWorldMDP(;size = (7, 5), tprob = 1.0, discount = 0.95, rewards = nothing, num_lava_tiles::Int = 1,
                        lava::Union{Symbol, Array{GWPos}} = [GWPos(1,1)],
                        goal::Union{Symbol, GWPos} = GWPos(7,5),
                        randomize_lava_on_reset::Bool = false,
                        possible_lava_squares::Union{Symbol, Array{GWPos}} = :all,
                        lava_penalty::Float64 = -1.,
                        goal_reward::Float64 = 1.,
                        observation_type = :all,
                        initial_state = :random,
                        rng::AbstractRNG = Random.GLOBAL_RNG)
         mdp = LavaWorld(gridworld = SimpleGridWorld(size = size, tprob = tprob, discount = discount, rewards = lava == :random ? Dict() : merge(Dict(p=>lava_penalty for p in lava), Dict(goal=>goal_reward))),
            num_lava_tiles = num_lava_tiles,
            possible_lava_squares = possible_lava_squares,
            goal = goal,
            randomize_lava_on_reset = randomize_lava_on_reset,
            lava_penalty = lava_penalty,
            goal_reward = goal_reward,
            observation_type = observation_type,
            rng = rng,
            initial_state = initial_state)
    lava == :random && randomize_lava!(mdp)
    mdp
end

POMDPs.gen(mdp::LavaWorld, s, a, rng = Random.GLOBAL_RNG) = (sp = rand(rng, transition(mdp.gridworld, s, a )), r = POMDPs.reward(mdp.gridworld, s, a))

function random_lava(mdp::LavaWorld)
    size, rng = mdp.gridworld.size, mdp.rng
    g = mdp.goal == :random ? GWPos(rand(rng, 1:size[1]), rand(rng, 1:size[2])) : mdp.goal
    lava = []
    while length(lava) < mdp.num_lava_tiles
        if mdp.possible_lava_squares == :all
            p = GWPos(rand(rng, 1:size[1]), rand(rng, 1:size[2]))
        else
            p = rand(rng, mdp.possible_lava_squares)
        end
        p != g && push!(lava, p)
    end
    @assert !(g in lava)
    rewards = Dict(p => mdp.lava_penalty for p in lava)
    rewards[g] = mdp.goal_reward
    rewards
end

function randomize_lava!(mdp::LavaWorld)
    g = mdp.gridworld
    [delete!(g.rewards, k) for k in keys(g.rewards)]
    [delete!(g.terminate_from, k) for k in g.terminate_from]
    rs = random_lava(mdp)
    for k in keys(rs)
        g.rewards[k] = rs[k]
        push!(g.terminate_from, k)
    end
end

function POMDPs.initialstate(mdp::LavaWorld)
    mdp.randomize_lava_on_reset && randomize_lava!(mdp)

    # return Deterministic(GWPos(1,5))
    if mdp.initial_state == :random
        function istate(rng::AbstractRNG)
            while true
                x, y = rand(rng, 1:mdp.gridworld.size[1]), rand(rng, 1:mdp.gridworld.size[2])
                !(GWPos(x,y) in mdp.gridworld.terminate_from) && return GWPos(x,y)
            end
        end
        return ImplicitDistribution(istate)
    else
        return Deterministic(mdp.initial_state)
    end
end

POMDPs.actions(mdp::LavaWorld) = [POMDPs.actions(mdp.gridworld)...]
POMDPs.states(mdp::LavaWorld) = states(mdp.gridworld)
POMDPs.reward(mdp::LavaWorld, s) = POMDPs.reward(mdp.gridworld, s)
POMDPs.isterminal(mdp::LavaWorld, s) = isterminal(mdp.gridworld, s)
POMDPs.discount(mdp::LavaWorld) = discount(mdp.gridworld)

function POMDPs.convert_s(::Type{V}, s::GWPos, mdp::LavaWorld) where {V<:AbstractArray}
    if mdp.observation_type == :all
        svec = zeros(Float32, mdp.gridworld.size..., 3, 1)
        !isterminal(mdp, s) && (svec[s[1], s[2], 3] = 1.)
        for p in states(mdp)
            POMDPs.reward(mdp, p) < 0 && (svec[p[1], p[2], 2] = 1.)
            POMDPs.reward(mdp, p) > 0 && (svec[p[1], p[2], 1] = 1.)
        end
        return svec
    else
        return [s...]
    end
end

function POMDPs.convert_s(::Type{GWPos}, v::V, mdp::LavaWorld) where {V<:AbstractArray}
    if mdp.observation_type == :all
        return GWPos(findfirst(reshape(v, mdp.gridworld.size..., :)[:,:,3] .== 1.0).I)
    else
        return v
    end
end

goal(mdp::LavaWorld, s) = GWPos(findfirst(reshape(s, mdp.gridworld.size..., :)[:,:,1] .== 1.0).I)

function gen_occupancy(buffer, mdp)
    N = length(buffer)
    occupancy = Dict(s => 0. for s in states(mdp))
    cols = fill(:, ndims(buffer[:s])-1)
    for i=1:N
        s = convert_s(GWPos, buffer[:s][cols..., i], mdp)
        occupancy[s] += 1 / N
    end
    occupancy
end

function render(mdp::LavaWorld, s=GWPos(7,5), a=nothing; color = s->10.0*POMDPs.reward(mdp, s), policy= nothing, return_compose = false)
    img = POMDPModelTools.render(mdp.gridworld, (s = s,), color = color, policy = isnothing(policy) ? nothing : FunctionPolicy((s) ->  action(policy, convert_s(AbstractArray, s, mdp))[1]))
    return_compose && return img
    tmpfilename = "/tmp/out.png"
    img |> PNG(tmpfilename, 1cm .* mdp.gridworld.size...)
    load(tmpfilename)
end

# render_and_save(filename, g::MDP...) = hcat_and_save(filename,  [POMDPModelTools.render(gi) for gi in g]...)
#
# function hcat_and_save(filename, c::Context...)
#     set_default_graphic_size(35cm,10cm)
#     r = compose(Compose.context(0,0,1cm, 0cm), Compose.rectangle()) # spacer
#     cs = []
#     for ci in c
#         push!(cs, ci)
#         push!(cs, r)
#     end
#     hstack(cs[1:end-1]...) |> PDF(filename)
# end

max_forks_repo_name: ancorso/POMDPGym


text:
using ACRLab

str = "
xinit = -0.5
theta1init = 0.1
"

ACRLab.run3!(str, vis = false)
@test true


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -7.107954588851326 <= q <= 8.439913276188518)

                     add_NL_constraint(m, :(softplus(-0.4953615570310643 + -0.5197659840811282*softplus(0.7941458217392698 + 0.49653685598476827*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.07122185861106178*softplus(-0.06420985369672083 + 0.9131512287471262*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.004105284291951339*softplus(-0.8124442669202567 + 0.4437121866105662*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + softplus(0.7344840166166842 + -0.5219551557665656*softplus(0.7941458217392698 + 0.49653685598476827*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.12237894168399821*softplus(-0.06420985369672083 + 0.9131512287471262*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.5945370485807278*softplus(-0.8124442669202567 + 0.4437121866105662*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + softplus(-0.7546740103129919 + -0.5743271753296177*softplus(0.7941458217392698 + 0.49653685598476827*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.959961363428667*softplus(-0.06420985369672083 + 0.9131512287471262*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.0685664132974031*softplus(-0.8124442669202567 + 0.4437121866105662*softplus(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*softplus(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*softplus(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module coolFFT

export coolfft

include("fft.jl")

end # module


text:
import Base.rand
import Distributions: logpdf

ProjDir = @__DIR__
cd(ProjDir)

struct LNR{T1,T2,T3} <: ContinuousUnivariateDistribution
    μ::T1
    σ::T2
    ϕ::T3
end

Broadcast.broadcastable(x::LNR) = Ref(x)

LNR(;μ, σ, ϕ) = LNR(μ, σ, ϕ)

function rand(dist::LNR)
    @unpack μ,σ,ϕ = dist
    x = @. rand(LogNormal(μ, σ)) + ϕ
    rt,resp = findmin(x)
    return resp,rt
end

rand(dist::LNR, N::Int) = [rand(dist) for i in 1:N]

function logpdf(d::T, r::Int, t::Float64) where {T<:LNR}
    @unpack μ,σ,ϕ = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == r
            LL += logpdf(LogNormal(m, σ), t-ϕ)
        else
            LL += log(1-cdf(LogNormal(m, σ), t-ϕ))
        end
    end
    return LL
end

logpdf(d::LNR, data::Tuple) = logpdf(d, data...)


text:
n = 8
function complex_least(n)
    w = rand(Complex{T},n)
    b = rand(Complex{T})
    U = rand(Complex{T},n,n)
    M = ctranspose(U)*U

    λ = b/(ctranspose(w)*inv(M)*w)
    α = λ*conj(inv(M))*conj(w)

    return (transpose(w)*α - b, transpose(α)*M*conj(α),α)
end

function complex_least_sys(i::Int)
    n=80; m=80;
    w = rand(Complex{T},m)
    c = rand(Complex{T},n)
    b = rand(Complex{T})
    A = rand(Complex{T},n,m)
    invM = inv(ctranspose(A)*A)

    λ = (b - transpose(w)*invM*ctranspose(A)*c)/(transpose(w)*invM*conj(w))

    α = λ*invM*conj(w) + invM*ctranspose(A)*c
    # α = invM*ctranspose(A)*c

    norm(transpose(w)*α - b)/norm(b)

    return norm(A*α - c)/norm(c)
end

# the wrong solution
function g(i::Int)
    n=80; m=80;
    w = rand(Complex{T},m)
    c = rand(Complex{T},n)
    b = rand(Complex{T})
    A = rand(Complex{T},n,m)
    invM = inv(ctranspose(A)*A)

    λ = (b - transpose(w)*invM*ctranspose(A)*c)/(transpose(w)*invM*(w))

    α = λ*invM*(w) + invM*ctranspose(A)*c
    # α = invM*ctranspose(A)*c

    norm(transpose(w)*α - b)/norm(b)
    norm(A*α - c)/norm(c)
end
mean(complex_least_sys(i) for i=1:500)
mean(g(i) for i=1:500)

plot([real.(A*α),real.(c)])

return norm(A*α - c)/norm(c)

ns = 2:2:1050
fs = zeros(Float64,length(ns))

for i in eachindex(ns)
    (constr, f, a) = complex_least(ns[i])
    fs[i] = abs(f)
end
# fs[2:end] = (fs[2:end]+fs[1:(length(fs)-1)])/2
v = exp.(-(-2:0.02:2).^2)
v = v/sum(v)
fs2 = conv(fs,v)
ns2 = conv(collect(ns),v)
fs2 = fs2[length(v):(length(fs2)-length(v))]
ns2 = ns2[length(v):(length(ns2)-length(v))]
# plot([fs2, fs[4:(length(fs)-3)]])
# plot([ns2, ns[4:(length(ns)-3)]])

plot(ns2,log.(fs2))