text:
include("../simulations/McSimulation.jl")

abstract type Payoff end
abstract type Leaf <: Payoff end
abstract type BinaryNode <: Payoff end

function at(self::Payoff, p::Path)
    throw(ArgumentError("Implementation of method at() required."))
end

function obsTime(self::Payoff)
    throw(ArgumentError("Implementation of method obsTime() required."))
end

function discountedAt(self::Payoff, p::Path)
    return at(self,p) / numeraire(p,obsTime(self))
end

function discountedAt(self::Payoff, paths::Array{Path,1})
    res = zeros(size(paths,1))
    @inbounds Threads.@threads for j = 1:size(res,1)
        res[j] = at(self,paths[j]) / numeraire(paths[j],obsTime(self))
    end
    return res
end

discountedAt(self::Payoff, sim::McSimulation) = discountedAt(self, paths(sim))

function discountedAt(payoffs::Array, paths::Array{Path,1})
    res = zeros( (size(payoffs,1), size(paths,1)) )
    @inbounds Threads.@threads for j = 1:size(paths,1)
        @inbounds Threads.@threads for i = 1:size(payoffs,1)
            res[i,j] = at(payoffs[i],paths[j]) / numeraire(paths[j],obsTime(payoffs[i]))
        end
    end
    return res
end

discountedAt(payoffs::Array{}, sim::McSimulation) = discountedAt(payoffs, paths(sim))

function observationTimes(self::Payoff)
    throw(ArgumentError("Implementation of method observationTimes() required."))
end

function observationTimes(self::Leaf)
    return Set([ obsTime(self) ])
end

function observationTimes(self::BinaryNode)
    return union(observationTimes(self.x),observationTimes(self.y))
end

function observationTimes(list::Array)
    times = Set{Float64}()
    for item in list
        times = union(times,observationTimes(item))
    end
    return sort([ t for t in times ])
end

# simplify payoff scripting

import Base.+
(+)(x::Payoff,y::Payoff) = Axpy(1.0,x,y)
(+)(x::Payoff,y) = Axpy(1.0,x,Fixed(y))
(+)(x,y::Payoff) = Axpy(1.0,Fixed(x),y)
#
import Base.-
(-)(x::Payoff,y::Payoff) = Axpy(-1.0,y,x)
(-)(x::Payoff,y) = Axpy(-1.0,Fixed(y),x)
(-)(x,y::Payoff) = Axpy(-1.0,y,Fixed(x))
#
import Base.*
(*)(x::Payoff,y::Payoff) = Mult(x,y)
(*)(x::Payoff,y) = Mult(x,Fixed(y))
(*)(x,y::Payoff) = Mult(Fixed(x),y)
#
import Base./
(/)(x::Payoff,y::Payoff) = Div(x,y)
(/)(x::Payoff,y) = Div(x,Fixed(y))
(/)(x,y::Payoff) = Div(Fixed(x),y)
#
import Base.%
(%)(x::Payoff,t) = Pay(x,t)
#
import Base.<
(<)(x::Payoff,y::Payoff) = Logical(x,y,(a,b)->Float64(a<b))
(<)(x::Payoff,y) = Logical(x,Fixed(y),(a,b)->Float64(a<b))
(<)(x,y::Payoff) = Logical(Fixed(x),y,(a,b)->Float64(a<b))
import Base.<=
(<=)(x::Payoff,y::Payoff) = Logical(x,y,(a,b)->Float64(a<=b))
(<=)(x::Payoff,y) = Logical(x,Fixed(y),(a,b)->Float64(a<=b))
(<=)(x,y::Payoff) = Logical(Fixed(x),y,(a,b)->Float64(a<=b))
import Base.==
(==)(x::Payoff,y::Payoff) = Logical(x,y,(a,b)->Float64(a==b))
(==)(x::Payoff,y) = Logical(x,Fixed(y),(a,b)->Float64(a==b))
(==)(x,y::Payoff) = Logical(Fixed(x),y,(a,b)->Float64(a==b))
import Base.>=
(>=)(x::Payoff,y::Payoff) = Logical(x,y,(a,b)->Float64(a>=b))
(>=)(x::Payoff,y) = Logical(x,Fixed(y),(a,b)->Float64(a>=b))
(>=)(x,y::Payoff) = Logical(Fixed(x),y,(a,b)->Float64(a>=b))
import Base.>
(>)(x::Payoff,y::Payoff) = Logical(x,y,(a,b)->Float64(a>b))
(>)(x::Payoff,y) = Logical(x,Fixed(y),(a,b)->Float64(a>b))
(>)(x,y::Payoff) = Logical(Fixed(x),y,(a,b)->Float64(a>b))

# basic payoffs

struct Fixed{T<:AbstractFloat} <: Leaf
    x::T
end

at(self::Fixed, p::Path) = self.x
obsTime(self::Fixed) = 0.0

#

struct Pay{T<:AbstractFloat} <: Payoff
    x::Payoff
    payTime::T
end

at(self::Pay, p::Path) = at(self.x, p)
obsTime(self::Pay) = self.payTime
function observationTimes(self::Pay)
    return union(observationTimes(self.x),[self.payTime])
end

#

struct Asset{T<:AbstractFloat} <: Leaf
    obsTime::T
    alias::String
end

Asset(obsTime) = Asset(obsTime,nothing)
at(self::Asset, p::Path) = asset(p,self.obsTime,self.alias)
obsTime(self::Asset) = self.obsTime

#

# basic rates payoffs

struct ZeroBond{T<:AbstractFloat} <: Leaf
    obsTime::T
    payTime::T
    alias::String
end

ZeroBond(obsTime,payTime) = ZeroBond(obsTime,payTime,"")
at(self::ZeroBond, p::Path) = zeroBond(p,self.obsTime,self.payTime,self.alias)
obsTime(self::ZeroBond) = self.obsTime

#

struct LiborRate{T<:AbstractFloat} <: Leaf
    obsTime::T
    startTime::T
    endTime::T
    yearFraction::T
    tenorBasis::T
    alias::String
    _dummy_::Bool  # avoid stack overflow during constructor
end

function LiborRate(obsTime,startTime,endTime;yearFraction=nothing,tenorBasis=1.0,alias=nothing)
    if isnothing(yearFraction)
        yearFraction = endTime - startTime
    end
    if isnothing(alias)
        alias = ""
    end
    return LiborRate(obsTime,startTime,endTime,yearFraction,tenorBasis,alias,true)
end

function at(self::LiborRate, p::Path)
    pStart = zeroBond(p, self.obsTime, self.startTime, self.alias)
    pEnd   = zeroBond(p, self.obsTime, self.endTime, self.alias)
    return (pStart/pEnd*self.tenorBasis - 1.0)/self.yearFraction
end
obsTime(self::LiborRate) = self.obsTime

#

struct SwapRate{T<:AbstractFloat} <: Leaf
    obsTime::T
    floatTimes::T
    floatWeights::T
    fixedTimes::T
    fixedWeights::T
    alias::String
end

SwapRate(obsTime,floatTimes,floatWeights,fixedTimes,fixedWeights) =
    SwapRate(obsTime,floatTimes,floatWeights,fixedTimes,fixedWeights,"")
function at(self::SwapRate, p::Path)
    num = sum([ w*zeroBond(p,self.obsTime,T,self.alias) for (w,T) in zip(self.floatWeights,self.floatTimes) ])
    den = sum([ w*zeroBond(p,self.obsTime,T,self.alias) for (w,T) in zip(self.fixedWeights,self.fixedTimes) ])
    return num / den
end
obsTime(self::SwapRate) = self.obsTime

#

struct FixedLeg{T<:AbstractFloat} <: Leaf
    obsTime::T
    payTimes::T
    payWeights::T
    alias::String
end

FixedLeg(obsTime,payTime,payWeights) = FixedLeg(obsTime,payTime,payWeights,"")
function at(self::FixedLeg, p::Path)
    return sum([ w*zeroBond(p,self.obsTime,T,self.alias) for (w,T) in zip(self.payWeights,self.payTimes) ])
end
obsTime(self::FixedLeg) = self.obsTime

# arithmetic operations

struct Axpy{T<:AbstractFloat} <: BinaryNode
    a::T
    x::Payoff
    y::Payoff
end

at(self::Axpy, p::Path) = self.a * at(self.x,p) + at(self.y,p)
obsTime(self::Axpy) = 0.0

#

struct Mult <: BinaryNode
    x::Payoff
    y::Payoff
end

at(self::Mult, p::Path) = at(self.x,p) * at(self.y,p)
obsTime(self::Mult) = 0.0

#

struct Div <: BinaryNode
    x::Payoff
    y::Payoff
end

at(self::Div, p::Path) = at(self.x,p) / at(self.y,p)
obsTime(self::Div) = 0.0

#

struct Max <: BinaryNode
    x::Payoff
    y::Payoff
end

at(self::Max, p::Path) = max(at(self.x,p),at(self.y,p))
obsTime(self::Div) = 0.0

#

struct Min <: BinaryNode
    x::Payoff
    y::Payoff
end

at(self::Min, p::Path) = min(at(self.x,p),at(self.y,p))
obsTime(self::Div) = 0.0

#

struct Logical <: BinaryNode
    x::Payoff
    y::Payoff
    op
end

at(self::Logical, p::Path) = self.op(at(self.x,p),at(self.y,p))
obsTime(self::Div) = 0.0

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
@doc raw"""
`SODE`: Split Ordinary Differential Equation

Defines an initial value problem
```math
\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,
```
with vector field ``v``, initial condition ``q_{0}`` and the solution
``q`` taking values in ``\mathbb{R}^{d}``. Here, the vector field ``v``
is given as a sum of vector fields
```math
v (t) = v_1 (t) + ... + v_r (t) .
```

### Fields

* `d`: dimension of dynamical variable ``q`` and the vector field ``v``
* `v`: tuple of functions computing the vector field
* `t₀`: initial time
* `q₀`: initial condition

The functions `v_i` providing the vector field must have the interface
```julia
    function v_i(t, q₀, q₁, h)
        q₁[1] = q₀[1] + ...
        q₁[2] = q₀[2] + ...
        ...
    end
```
where `t` is the current time, `q₀` is the current solution vector, `q₁` is the
new solution vector which holds the result of computing one substep with the
vector field ``v_i`` on `t` and `q₀`, and `h` is the (sub-)timestep to compute
the update for.

The fact that the function `v` returns the solution and not just the vector
field for each substep increases the flexibility for the use of splitting
methods, e.g., it allows to use another integrator for solving substeps.

"""
struct SODE{dType <: Number, tType <: Number,
            vType <: Tuple, pType <: Union{Tuple,Nothing}, N} <: AbstractEquationODE{dType, tType}

    d::Int
    n::Int
    v::vType
    t₀::tType
    q₀::Array{dType,N}
    parameters::pType
    periodicity::Vector{dType}

    function SODE(DT::DataType, N::Int, d::Int, n::Int, v::vType, t₀::tType, q₀::AbstractArray{dType};
                 parameters=nothing, periodicity=zeros(DT,d)) where {
                        dType <: Number, tType <: Number, vType <: Tuple}

        @assert d == size(q₀,1)
        @assert n == size(q₀,2)
        @assert ndims(q₀) == N ∈ (1,2)

        new{DT, tType, vType, typeof(parameters), N}(d, n, v, t₀,
                convert(Array{DT}, q₀), parameters, periodicity)
    end
end

function SODE(v, t₀, q₀::AbstractArray{DT}; kwargs...) where {DT}
    SODE(DT, ndims(q₀), size(q₀,1), size(q₀,2), v, t₀, q₀; kwargs...)
end

function SODE(v, q₀; kwargs...)
    SODE(v, zero(eltype(q₀)), q₀; kwargs...)
end

Base.hash(ode::SODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.v, hash(ode.t₀,
        hash(ode.q₀, hash(ode.periodicity, hash(ode.parameters, h)))))))

Base.:(==)(ode1::SODE, ode2::SODE) = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.v == ode2.v
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.parameters == ode2.parameters
                             && ode1.periodicity == ode2.periodicity)

function Base.similar(ode::SODE, q₀; kwargs...)
    similar(ode, ode.t₀, q₀; kwargs...)
end

function Base.similar(ode::SODE, t₀::TT, q₀::AbstractArray{DT};
                      parameters=ode.parameters, periodicity=ode.periodicity) where {DT  <: Number, TT <: Number}
    @assert ode.d == size(q₀,1)
    SODE(ode.v, t₀, q₀; parameters=parameters, periodicity=periodicity)
end

@inline Base.ndims(ode::SODE) = ode.d

@inline periodicity(equation::SODE) = equation.periodicity

ext: jl


text:
using Base: @propagate_inbounds

using OffsetArrays
using JLD2

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Fields

using Oceananigans.Grids: topology, total_size, interior_parent_indices
using Oceananigans.Fields: show_location

import Oceananigans: short_show
import Oceananigans.Fields: Field, set!, interior

struct FieldTimeSeries{X, Y, Z, K, A, T, D, G, B, χ} <: AbstractDataField{X, Y, Z, A, G, T, 4}
                   data :: D
           architecture :: A
                   grid :: G
    boundary_conditions :: B
                  times :: χ

    function FieldTimeSeries{X, Y, Z, K}(data::D, arch::A, grid::G, bcs::B, times::χ) where {X, Y, Z, K, D, A, G, B, χ}
        T = eltype(grid)
        return new{X, Y, Z, K, A, T, D, G, B, χ}(data, arch, grid, bcs, times)
    end
end

#####
##### Constructors
#####

"""
    FieldTimeSeries{LX, LY, LZ}([architecture = CPU()], grid, times, boundary_conditions=nothing)

Return `FieldTimeSeries` at location `(LX, LY, LZ)`, on `grid`, at `times`, with
`boundary_conditions`, and initialized with zeros of `eltype(grid)`.
"""
function FieldTimeSeries{LX, LY, LZ}(architecture, grid, times, boundary_conditions=nothing) where {LX, LY, LZ}
    location = (LX, LY, LZ)
    Nt = length(times)
    data_size = total_size(location, grid)
    raw_data = zeros(architecture, grid, data_size..., Nt)
    data = offset_data(raw_data, grid, location)
    return FieldTimeSeries{LX, LY, LZ, InMemory}(data, architecture, grid, boundary_conditions, times)
end

# CPU() default
FieldTimeSeries{LX, LY, LZ}(grid::AbstractGrid, times, bcs=nothing) where {LX, LY, LZ} =
    FieldTimeSeries{LX, LY, LZ}(CPU(), grid, times, bcs)

"""
    FieldTimeSeries(path, name;
                    architecture = CPU(),
                    backend = InMemory(),
                    grid = nothing,
                    iterations = nothing,
                    times = nothing)

Returns a `FieldTimeSeries` for the field `name` describing a field's time history from a JLD2 file
located at `path`.

Keyword arguments
=================

- `architecture`: The architecture on which to store time series data. `CPU()` by default.

- `backend`: `InMemory()` to load data into a 4D array or `OnDisk()` to lazily load data from disk
             when indexing into `FieldTimeSeries`.

- `grid`: A grid to associated with data, in the case that the native grid
          was not serialized properly.

- `iterations`: Iterations to load. Defaults to all iterations found in the file.

- `times`: Save times to load, as determined through an approximate floating point
           comparison to recorded save times. Defaults to times associated with `iterations`.
           Takes precedence over `iterations` if `times` is specified.
"""
FieldTimeSeries(path, name; architecture=CPU(), backend=InMemory(), kwargs...) =
    FieldTimeSeries(path, name, architecture, backend; kwargs...)

#####
##### InMemory time serieses
#####

const InMemoryFieldTimeSeries{X, Y, Z} = FieldTimeSeries{X, Y, Z, InMemory}

struct UnspecifiedBoundaryConditions end

function FieldTimeSeries(path, name, architecture, backend::InMemory;
                         grid = nothing,
                         location = nothing,
                         boundary_conditions = UnspecifiedBoundaryConditions(),
                         iterations = nothing,
                         times = nothing)

    file = jldopen(path)

    # Non-defaults
    isnothing(grid)       && (grid       =  file["serialized/grid"])
    isnothing(iterations) && (iterations =  parse.(Int, keys(file["timeseries/t"])))
    isnothing(times)      && (times      =  [file["timeseries/t/$i"] for i in iterations])
    isnothing(location)   && (location   =  file["timeseries/$name/serialized/location"])

    if boundary_conditions isa UnspecifiedBoundaryConditions
        boundary_conditions = file["timeseries/$name/serialized/boundary_conditions"]
    end

    LX, LY, LZ = location

    time_series = FieldTimeSeries{LX, LY, LZ}(architecture, grid, times, boundary_conditions)

    set!(time_series, path, name)

    return time_series
end

Base.getindex(fts::InMemoryFieldTimeSeries{LX, LY, LZ}, n::Int) where {LX, LY, LZ} =
    Field(LX, LY, LZ, fts.architecture, fts.grid, fts.boundary_conditions, view(fts.data, :, :, :, n))

#####
##### set!
#####

"""
    Field(path::String, name::String, iter; architecture=GPU(), grid=nothing)

Load a Field saved in JLD2 file at `path`, with `name` and at `iter`ation.
`architecture = CPU()` by default, and `grid` is loaded from `path` if not specified.
"""
function Field(location, path::String, name::String, iter;
               architecture = CPU(),
               grid = nothing,
               boundary_conditions = nothing)

    file = jldopen(path)

    raw_data = arch_array(architecture, file["timeseries/$name/$iter"])
    isnothing(grid) && (grid = file["serialized/grid"])

    close(file)

    try
        data = offset_data(raw_data, grid, location)
        return Field(location, architecture, grid, boundary_conditions, data)
    catch
        field = Field(location, architecture, grid, boundary_conditions)
        interior(field) .= raw_data
        return field
    end
end

function set!(time_series::InMemoryFieldTimeSeries, path::String, name::String)

    file = jldopen(path)
    file_iterations = parse.(Int, keys(file["timeseries/t"]))
    file_times = [file["timeseries/t/$i"] for i in file_iterations]
    close(file)

    for (n, time) in enumerate(time_series.times)
        file_index = findfirst(t -> t ≈ time, file_times)
        file_iter = file_iterations[file_index]

        field_n = Field(location(time_series), path, name, file_iter,
                        boundary_conditions = time_series.boundary_conditions,
                        grid = time_series.grid)

        set!(time_series[n], field_n)
    end

    return nothing
end

function set!(time_series::FieldTimeSeries, fields_vector::AbstractVector{<:AbstractField})
    raw_data = parent(time_series.data)
    ArrayType = array_type(time_series.architecture)

    file = jldopen(path)

    for (n, field) in enumerate(fields_vector)
        raw_data[:, :, :, n] .= parent(field)
    end

    close(file)

    return nothing
end

# TODO: this is a bit of type-piracy (with respect to the Oceananigans.Fields module)...
# is there a better way?

# FieldTimeSeries[i] returns ViewField
const ViewField = Field{<:Any, <:Any, <:Any, <:Any, <:SubArray}

using OffsetArrays: IdOffsetRange

parent_indices(idx::Int) = idx
parent_indices(idx::Base.Slice{<:IdOffsetRange}) = Colon()

# Is this too surprising?
Base.parent(vf::ViewField) = view(parent(parent(vf.data)), parent_indices.(vf.data.indices)...)

"Returns a view of `f` that excludes halo points."
@inline interior(f::FieldTimeSeries{X, Y, Z}) where {X, Y, Z} =
    view(parent(f.data),
         interior_parent_indices(X, topology(f, 1), f.grid.Nx, f.grid.Hx),
         interior_parent_indices(Y, topology(f, 2), f.grid.Ny, f.grid.Hy),
         interior_parent_indices(Z, topology(f, 3), f.grid.Nz, f.grid.Hz),
         :)

#####
##### OnDisk time serieses
#####

struct OnDiskData
    path :: String
    name :: String
end

function FieldTimeSeries(path, name, architecture, backend::OnDisk; grid=nothing)
    file = jldopen(path)

    if isnothing(grid)
        grid = file["serialized/grid"]
    end

    iterations = parse.(Int, keys(file["timeseries/t"]))
    times = [file["timeseries/t/$i"] for i in iterations]

    data = OnDiskData(path, name)
    LX, LY, LZ = file["timeseries/$name/serialized/location"]
    bcs = file["timeseries/$name/serialized/boundary_conditions"]

    close(file)

    return FieldTimeSeries{LX, LY, LZ, OnDisk}(data, architecture, grid, bcs, times)
end

# For creating an empty `FieldTimeSeries`.
function FieldTimeSeries(grid, location, times; architecture=CPU(), ArrayType=array_type(architecture), name="", filepath="", bcs=nothing)
    LX, LY, LZ = location

    Nt = length(times)
    data_size = total_size(location, grid)

    raw_data = zeros(data_size..., Nt) |> ArrayType
    data = offset_data(raw_data, grid, location)

    return FieldTimeSeries{LX, LY, LZ}(InMemory(), data, architecture, grid, bcs, times, name, filepath, 4)
end

#####
##### Methods
#####

# Include the time dimension.
@inline Base.size(fts::FieldTimeSeries) = (size(location(fts), fts.grid)..., length(fts.times))

@propagate_inbounds Base.getindex(f::FieldTimeSeries{X, Y, Z, InMemory}, i, j, k, n) where {X, Y, Z} = f.data[i, j, k, n]

function Base.getindex(fts::FieldTimeSeries{X, Y, Z, OnDisk}, n::Int) where {X, Y, Z}
    # Load data
    file = jldopen(fts.data.path)
    iter = keys(file["timeseries/t"])[n]
    raw_data = file["timeseries/$(fts.data.name)/$iter"] |> array_type(fts.architecture)
    close(file)

    # Wrap Field
    loc = (X, Y, Z)
    field_data = offset_data(raw_data, fts.grid, loc)

    return Field(loc..., fts.architecture, fts.grid, fts.boundary_conditions, field_data)
end

Base.setindex!(fts::FieldTimeSeries, val, inds...) = Base.setindex!(fts.data, val, inds...)

#####
##### Show methods
#####

backend_str(::InMemory) = "InMemory"
backend_str(::OnDisk) = "OnDisk"

#####
##### show
#####

short_show(fts::FieldTimeSeries{X, Y, Z, K}) where {X, Y, Z, K} =
    string("$(join(size(fts), "×")) FieldTimeSeries{$(backend_str(K()))} located at $(show_location(fts))")

Base.show(io::IO, fts::FieldTimeSeries{X, Y, Z, K, A}) where {X, Y, Z, K, A} =
    print(io, "$(short_show(fts))\n",
          "├── architecture: $A\n",
          "└── grid: $(short_show(fts.grid))")


text:
# Read in NexSci planets with non-zero lower limits
# and both masses & radii measured.

if ~@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end

using DelimitedFiles; using PyPlot

clf()
figure(figsize=(6,6))
datadai = readdlm("../../../data/apjab3a3bt2_ascii.txt",',',comments=true,comment_char='#')
datadai[:,10] .*= -1
datadai[:,13] .*= -1
#datadaiflt = convert(Array{Float64,2},datadai[:,9:14])
#errorbar(datadai[:,9],datadai[:,12],xerr=permutedims(datadai[:,10:11]),yerr=permutedims(datadai[:,13:14]),fmt="o")
offsetx = zeros(13) .+ 0.35
#offsetx[13] *= 0.6  # K-93b
#offsetx[7] *= 1.1  # K2-141b
xerrs = reverse(permutedims(datadai[:,10:11]),dims=1)
yerrs = reverse(permutedims(datadai[:,13:14]),dims=1)
for i=1:13
  errorbar(datadai[i,9],datadai[i,12],xerr=reshape(xerrs[:,i],2,1),yerr=reshape(yerrs[:,i],2,1),fmt="o",alpha=0.18/minimum(xerrs[:,i]),color="C3",linewidth=2.0)
#  text([datadai[i,9]] .+ offsetx[i],[datadai[i,12]-0.01],datadai[i,1],fontsize=8)
#  text([1.0],[datadai[i,12]-0.01+offsety[i]],datadai[i,1],fontsize=4)
#  text([datadai[i,12]^3.7*offsetx[i]],[datadai[i,12]-0.01],datadai[i,1],fontsize=8)
#  text([datadai[i,12]^3.7*offsetx[i]],[datadai[i,12]-0.01],datadai[i,1],fontsize=5,rotation=90)
  println(i," ",offsetx[i]," ",datadai[i,1])
end
#dataeu = readdlm("exoplanet.eu_catalog.csv",',',skipstart=1,comments=true,comment_char='#')

#errorbar(dataeu[:,3].*(MJUPITER/MEARTH),dataeu[:,9].*(RJUPITER/REARTH),xerr=permutedims(dataeu[:,4:5].*(MJUPITER/MEARTH)),
#  yerr=permutedims(data[:,10:11].*(RJUPITER/REARTH)),fmt="o")

data = readdlm("../../../data/planets_2020.02.26_11.12.22.csv",',',skipstart=34,comments=true,comment_char='#')
data[:,14] .*= -1
data[:,18] .*= -1
#igood = (data[:,12] .> (5data[:,14])) .& (data[:,16] .> 12.5data[:,17]) .& (data[:,16] .> 12.5data[:,18])
igood = data[:,12] .> 5data[:,14]
#xerrs = reverse(permutedims(data[igood,13:14]),dims=1)
xerrs = reverse(permutedims(data[:,13:14]),dims=1)
#yerrs = reverse(permutedims(data[igood,17:18]),dims=1)
yerrs = reverse(permutedims(data[:,17:18]),dims=1)
#errorbar(data[igood,12],data[igood,16]./data[igood,12].^(1/3.7),xerr=xerrs,yerr=yerrs,fmt=".")
offset = zeros(31) .+0.35
#offset[15] *= 0.7 # Kepler-105c
#offset[13] *= 1.55 # GJ 357b
#offset[23] *= 0.75  # Kepler-36b
#offset[22] *= 0.8  # Kepler-10b
#offset[24] *= 1.2  # HD 219134 c
#offset[13] = -0.8; offset[18] = 0.6; offset[23] = -2.3; offset[25] = -2.4;  offset[24] = 0.5; offset[30] = -2.6; offset[27] = -2.4; offset[15] = -2.8; offset[31] = 0.7
for i=1:31
  if igood[i]
    if i != 23
      errorbar(data[i,12],data[i,16],xerr=reshape(xerrs[:,i],2,1),yerr=reshape(yerrs[:,i],2,1),fmt="o",alpha=0.18/minimum(xerrs[:,i]),color="C3",linewidth=2.0)
    else
      errorbar(data[i,12],data[i,16],xerr=reshape(xerrs[:,i],2,1),yerr=reshape(yerrs[:,i],2,1),fmt="o",label="NexSci/Dai/Dressing",alpha=0.18/minimum(xerrs[:,i]),color="C3",linewidth=2.0)
    end
#    text([data[i,12]] .+ offset[i],[data[i,16]-0.01],string(data[i,1]," ",data[i,2]),fontsize=8)
#    text([1.0],[data[i,16]-0.01],string(data[i,1]," ",data[i,2]),fontsize=4)
#    text([data[i,16]^3.7*offset[i]],[data[i,16]-0.01],string(data[i,1]," ",data[i,2]),fontsize=5,rotation=90)
    println(i," ",offset[i]," ",data[i,1]," ",data[i,2])
  end
end
m = collect(0.05:0.01:10.0)
#semilogx(m,m.^(1/3.7)./m.^(1/3.7))
#semilogx(m,m.^(1/3.7),label="CMF=0.33")
mrearth = readdlm("../../../data/massradiusEarthlikeRocky.txt")
semilogx(mrearth[9:31,1],mrearth[9:31,2],label="Earth-like (Zeng et al. 2019)")
#mrearth_case5 = readdlm("zeng2013_tab02_case2.txt")
#plot(mrearth_case5[10:29,1],mrearth_case5[10:29,2],label="Case 2 CMF=0.1787 (Zeng & Sasselov 2013)")
# Read in the ternary diagram (Table 3) from Zeng & Sasselov (2013), but
# without any of the water planet cases:
#mrearth_cmf = readdlm("zeng2013_tab03_dry.txt",skipstart=1)
#mrearth_cmf = reshape(mrearth_cmf,5,41,8)
## Find a 17% CMF case:
#mrearth_17pct = zeros(2,41)
#mrearth_17pct[1,:] = mrearth_cmf[1,:,1]
#pctfe = 0.17
#for i=1:41
#  for j=1:4
#    if (mrearth_cmf[j,i,6] > pctfe) && (mrearth_cmf[j+1,i,6] < pctfe)
#      x = (mrearth_cmf[j,i,6] - pctfe)/(mrearth_cmf[j,i,6] - mrearth_cmf[j+1,i,6])
#      mrearth_17pct[2,i] = mrearth_cmf[j,i,2] + (mrearth_cmf[j+1,i,2] - mrearth_cmf[j,i,2])*x
#    end
#  end
#end
#plot(mrearth_17pct[1,1:30],mrearth_17pct[2,1:30],label="Zeng & Sasselov (2013); 17% CMF")
#semilogx(m,m.^(1/3.7)*(1.07-0.21*0.25)./m.^(1/3.7))
#semilogx(m,m.^(1/3.7)*(1.07-0.21*0.25),label="CMF=0.25")

# Now plot trappist:

mrt1 = readdlm("../../../data/t1_mass_radius.txt",',')
#errorbar(mrt1[:,2],mrt1[:,4]./mrt1[:,2].^(1/3.7),xerr=mrt1[:,3],yerr=mrt1[:,5],fmt=".")
errorbar(mrt1[:,2],mrt1[:,4],xerr=mrt1[:,3],yerr=mrt1[:,5],fmt=".",label="Trappist-1",linewidth=2.0)

# Plot Earth, Venus, Mars, Mercury(?):
mss = [5.9736e27,0.3302e27,4.8685e27,0.64185e27]./5.9736e27
rss = [6371.00e5,2439.7e5,6051.8e5,3389.5e5]./6371.0e5

#plot(mss,rss./mss.^(1/3.7),"o")
plot(mss,rss,"o",label="Solar System")

#dmr_carbon = readdlm("../Figures/Dorn/M_R_carbon-rich_Eric.ddat")
dmr_nocore_solar = readdlm("../../../data/MR_trappist_corefree_Solar.ddat",skipstart=1)
#dmr_nocore_unter = readdlm("../Figures/Dorn/MR_trappist_corefree_Unterborn.ddat",skipstart=1)

#plot(dmr_carbon[:,1],dmr_carbon[:,2],label="Carbon-rich (Miozzi et al. 2018)")
plot(dmr_nocore_solar[:,1],dmr_nocore_solar[:,2],label="No-core, Solar",color="C4")
#plot(dmr_nocore_unter[:,1],dmr_nocore_unter[:,2],label="No-core, U17")

axis([0.04,10.0,0.35,2.0])

xlabel(L"$M/M_\oplus$")
#ylabel(L"$(R/R_\oplus)(M_\oplus/M)^{0.27}$")
ylabel(L"$R/R_\oplus$")
legend()
tight_layout()
text([0.28],[0.77],"h")
text([0.35],[0.81],"d")
text([0.64],[0.94],"e")
text([0.96],[1.05],"f")
text([1.22],[1.15],"g")
text([1.15],[1.08],"c")
text([1.40],[1.125],"b")

text([1.07],[0.98],"Earth")
text([0.88],[0.93],"Venus")
text([0.12],[0.52],"Mars")
text([0.064],[0.38],"Mercury")
savefig("../mass_radius_relation_comparison.pdf",bbox_inches="tight")


text:
######################################################################
# weierstrass.jl : curves given by weierstrass equations
######################################################################

export Weierstrass, ShortWeierstrass, SeparatedWeierstrass, EllipticCurve, from_j_invariant

export curvetype, a_invariants, b_invariants, c_invariants, discriminant, j_invariant, equation, projective_equation, divisionpolynomial, coordinate_ring

######################################################################
# Type definitions
######################################################################

"""
Concrete type for elliptic curves in long Weierstrass form over a ring.
"""
struct Weierstrass{T} <: AbstractWeierstrass{T}
    a1::T
    a2::T
    a3::T
    a4::T
    a6::T
end

"""
Concrete types for elliptic curves in short Weierstrass form over a ring.
"""
struct ShortWeierstrass{T} <: AbstractWeierstrass{T}
    a::T
    b::T
end

struct SeparatedWeierstrass{T} <: AbstractWeierstrass{T}
	a2::T
	a4::T
	a6::T
end

curvetype(E::SeparatedWeierstrass) = SeparatedWeierstrass

curvetype(E::ShortWeierstrass) = ShortWeierstrass

curvetype(E::Weierstrass) = Weierstrass

######################################################################
# Calling with six arguments
######################################################################

function ShortWeierstrass(a1::T, a2::T, a3::T, a4::T, a6::T) where T
	a1 == 0 && a2 == 0 && a3 == 0 || throw(ArgumentError("Tried to call ShortWeierstrass with nonzero coefficients"))
	return ShortWeierstrass(a4, a6)
end

function SeparatedWeierstrass(a1::T, a2::T, a3::T, a4::T, a6::T) where T
	a1 == 0 && a3 == 0 || throw(ArgumentError("Tried to call SeparatedWeierstrass with nonzero coefficients"))
	return SeparatedWeierstrass(a2, a4, a6)
end

######################################################################
# Calling with the 'EllipticCurve' constructor
######################################################################

function EllipticCurve(a4::T, a6::T) where T
	return ShortWeierstrass(a4, a6)
end

function EllipticCurve(a2::T, a4::T, a6::T) where T
	return SeparatedWeierstrass(a2, a4, a6)
end

function EllipticCurve(a1::T, a2::T, a3::T, a4::T, a6::T) where T
	return Weierstrass(a1, a2, a3, a4, a6)
end

function from_j_invariant(j::T) where T
    a = -j * (j - 1728)
    b = -a * (j - 1728)
    return ShortWeierstrass(3*a, 2*b)
end

######################################################################
# Basic methods
######################################################################

base_ring(E::Weierstrass) = Nemo.parent(E.a1)

base_ring(E::ShortWeierstrass) = Nemo.parent(E.a)

base_ring(E::SeparatedWeierstrass) = Nemo.parent(E.a2)

"""
Get a description of an elliptic curve in long Weierstrass form.
"""
function show(io::IO, E::Weierstrass)
    print(io, "Elliptic Curve in long Weierstrass form y² + $(E.a1) xy + $(E.a3) y = x³ + $(E.a2) x² + $(E.a4) x + $(E.a6)  over ")
    show(io, base_ring(E))
end

"""
Get a description of an elliptic curve in short Weierstrass form.
"""
function show(io::IO, E::ShortWeierstrass)
    print(io, "Elliptic Curve in short Weierstrass form y² = x³ + $(E.a) x + $(E.b)  over ")
    show(io, base_ring(E))
end

"""
Get a description of an elliptic curve in short Weierstrass form.
"""
function show(io::IO, E::SeparatedWeierstrass)
    print(io, "Elliptic Curve in separated Weierstrass form y² = x³ + $(E.a2) x² + $(E.a4) x + $(E.a6)  over ")
    show(io, base_ring(E))
end

######################################################################
# Invariants and equations
######################################################################

"""
Get the a-invariants (a1, a2, a3, a4, a6) of an elliptic curve in long Weierstrass form.
"""
a_invariants(E::Weierstrass) = (E.a1, E.a2, E.a3, E.a4, E.a6)

"""
Get the a-invariants (a1, a2, a3, a4, a6) of an elliptic curve in short Weierstrass form.
"""
function a_invariants(E::ShortWeierstrass)
    zero = Nemo.zero(base_ring(E))
    return (zero, zero, zero, E.a, E.b)
end

"""
Get the a-invariants (a1, a2, a3, a4, a6) of an elliptic curve in separated Weierstrass form.
"""
function a_invariants(E::SeparatedWeierstrass)
    zero = Nemo.zero(base_ring(E))
    return (zero, E.a2, zero, E.a4, E.a6)
end

==(E1::AbstractWeierstrass, E2::AbstractWeierstrass) = a_invariants(E1) == a_invariants(E2)

"""
Get the b-invariants (b2, b4, b6, b8) of an elliptic curve in Weierstrass form.
"""
function b_invariants(E::AbstractWeierstrass)
    a1, a2, a3, a4, a6 = a_invariants(E)
    return (a1*a1 + 4*a2,
             a1*a3 + 2*a4,
             a3^2 + 4*a6,
             a1^2 * a6 + 4*a2*a6 - a1*a3*a4 + a2*a3^2 - a4^2)
end

"""
Get the c-invariants (c4, c6) of an elliptic curve in Weierstrass form.
"""
function c_invariants(E::AbstractWeierstrass)
    b2, b4, b6, b8 = b_invariants(E)
    return (b2^2 - 24*b4, -b2^3 + 36*b2*b4 - 216*b6)
end

"""
Get the discriminant of an elliptic curve in Weierstrass form.
"""
function discriminant(E::AbstractWeierstrass)
    b2, b4, b6, b8 = b_invariants(E)
    return -b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6
end

"""
Check if a given instance of AbstractWeierstrass is indeed a nonsingular curve.
"""
isvalid(E::AbstractWeierstrass) = !iszero(discriminant(E))

"""
Get the j-invariant of an elliptic curve in Weierstrass form.
This may fail if the base ring is not a field.
"""
function j_invariant(E::AbstractWeierstrass)
    c4, _ = c_invariants(E)
    disc = discriminant(E)
	return c4^3 // disc
end

function equation(E::AbstractWeierstrass) #use nested polynomial rings for efficient reduction
	K = base_ring(E)
	RX, X = PolynomialRing(K, "X")
	RXY, Y = PolynomialRing(RX, "Y")
	a1, a2, a3, a4, a6 = a_invariants(E)
	return Y^2 + a1 * X * Y + a3 * Y - X^3 - a2 * X^2 - a4 * X - a6
end

function equation(E::AbstractWeierstrass, X, Y)
	a1, a2, a3, a4, a6 = a_invariants(E)
	return Y^2 + a1 * X * Y + a3 * Y - X^3 - a2 * X^2 - a4 * X - a6
end

function projective_equation(E::AbstractWeierstrass)
	K = base_ring(E)
	R, (X, Y, Z) = PolynomialRing(K, ["X", "Y", "Z"])
	a1, a2, a3, a4, a6 = a_invariants(E)
	return Z * Y^2 + a1 * X * Y * Z + a3 * Y * Z^2 - X^3 - a2 * X^2 * Z - a4 * X * Z^2 - a6 * Z^3
end

######################################################################
# Division polynomials for Weierstrass curves
######################################################################

"""
Internal function to compute the 2-division polynomial of an elliptic curve in Weierstrass form.
"""
_2divpoly(E::AbstractWeierstrass{T}, x::PolyElem{T}, D::Dict{Int, PolyElem{T}}) where T = _divpoly(E, -1, x, D)

"""
Internal function to compute the mth division polynomial of an elliptic curve in Weierstrass form, using the well-known induction formulas. If the index is even, *the usual factor 2*y is dropped*.
"""
function _divpoly(E::AbstractWeierstrass{T}, m::Integer, x::PolyElem{T}, D::Dict{Int, PolyElem{T}}) where T
	b2, b4, b6, b8 = b_invariants(E)
	R = parent(x)
	if m in keys(D)
		return D[m]
	elseif m == -2
		D[-2] = _divpoly(E, - 1, x, D)^2
		return D[-2]
	elseif m == -1
		D[-1] = 4 * x^3 + b2 * x^2 + 2 * b4  *x + b6
		return D[-1]
	elseif m <= 0
		throw(ArgumentError("Non-positive integers in _divpoly must be either -1 or -2, got $m"))
	elseif m == 1
		D[1] = R(1)
		return D[1]
	elseif m == 2
		D[2] = R(1)
		return D[2]
	elseif m == 3
		D[3] = 3 * x^4 + b2 * x^3 + 3 * b4 * x^2 + 3 * b6 * x + b8
		return D[3]
	elseif m == 4
		pm2 = _divpoly(E, -2, x, D)
		p3 = _divpoly(E, 3, x, D)
		D[4] = - pm2 + (6*x^2 + b2*x + b4) * p3
		return D[4]
	elseif isodd(m)
		k = div(m - 1, 2)
		m1 = _divpoly(E, k - 1, x, D)
		eq = _divpoly(E, k, x, D)
		p1 = _divpoly(E, k + 1, x, D)
		p2 = _divpoly(E, k + 2, x, D)
		if iseven(k)
			D[m] = _divpoly(E, -2, x, D) * p2 * eq^3 - m1 * p1^3
			return D[m]
		else #k is odd
			D[m] = p2 * eq^3 - _divpoly(E, -2, x, D) * m1 * p1^3
			return D[m]
		end
	else  #m is even
		k = div(m, 2)
		m2 = _divpoly(E, k - 2, x, D)
		m1 = _divpoly(E, k - 1, x, D)
		eq = _divpoly(E, k, x, D)
		p1 = _divpoly(E, k + 1, x, D)
		p2 = _divpoly(E, k + 2, x, D)
		D[m] = eq * (p2 * m1^2 - m2 * p1^2)
		return D[m]
	end
end

"""
Compute the division polynomial of an elliptic curve in Weierstrass form.

If m is even, you may specify an additional parameter two_torsion:
* 0 simply ignores the points of order two
* 1 returns the exact bivariate division polynomial
* 2 (default) adds extra roots at points of two-torsion in order to get a univariate polynomial.
"""

function divisionpolynomial(E::AbstractWeierstrass{T}, m::Nemo.Integer, two_torsion::Integer = 2) where T
	m <= -3 && throw(ArgumentError("m must be positive or -1 or -2"))

	if m<=0
		K = base_ring(E)
		R, x = Nemo.PolynomialRing(K, "x")
		D = Dict{Int, Nemo.PolyElem{T}}()
		poly = _divpoly(E, m, x, D)
		return poly
	elseif m%2 == 0
		if two_torsion == 1
			K = base_ring(E)
			a1, _, a3, _, _ = a_invariants(E)
			R, (x, y) = Nemo.PolynomialRing(K, ["x", "y"])
			D = Dict{Int, Nemo.PolyElem{T}}()
			poly = _divpoly(E, m, x, D)
			return (2 * y + a1 * x + a3) * poly
		elseif two_torsion == 2
			K = base_ring(E)
			R, x = Nemo.PolynomialRing(K, "x")
			D = Dict{Int, Nemo.PolyElem{T}}()
			poly = _divpoly(E, m, x, D)
			pm1 = _divpoly(E, -1, x, D)
			return pm1 * poly
		else #two_torsion == 0
			K = base_ring(E)
			R, x = Nemo.PolynomialRing(K, "x")
			D = Dict{Int, Nemo.PolyElem{T}}()
			poly = _divpoly(E, m, x, D)
			return poly
		end
	else # m is odd
		K = base_ring(E)
		R, x = Nemo.PolynomialRing(K, "x")
		D = Dict{Int, Nemo.PolyElem{T}}()
		poly = _divpoly(E, m, x, D)
		return poly
	end
end

######################################################################
# Coordinate ring
######################################################################

function coordinate_ring(E::AbstractWeierstrass)
	Eq = equation(E)
	Y = gen(parent(Eq))
	X = gen(base_ring(Eq))
	RXY = parent(Eq)
	return ResidueRing(RXY, Eq), X, Y
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DouglasPeucker(ϵ)

The Douglas-Peucker algorithm for simplifying polygonal chains given
a deviation tolerance `ϵ`.

## References

* Douglas, D. and Peucker, T. 1973. [Algorithms for the Reduction of
  the Number of Points Required to Represent a Digitized Line or its
  Caricature](https://www.sciencedirect.com/science/article/abs/pii/0167839691900198)
"""
struct DouglasPeucker{T} <: SimplificationMethod
  ϵ::T
end

function simplify(polyarea::PolyArea, method::DouglasPeucker)
  c = [simplify(chain, method) for chain in chains(polyarea)]
  PolyArea(c[1], c[2:end])
end

function simplify(chain::Chain, method::DouglasPeucker)
  v = simplify(vertices(chain), method)
  isclosed(chain) ? Chain([v; first(v)]) : Chain(v)
end

# simplify chain assuming it is open
function simplify(v::AbstractVector{Point{Dim,T}},
                  method::DouglasPeucker) where {Dim,T}
  # find vertex with maximum distance
  imax, dmax = 0, zero(T)
  for i in 2:length(v)-1
    d = evaluate(Euclidean(), v[i], Line(first(v), last(v)))
    if d > dmax
      imax = i
      dmax = d
    end
  end

  if dmax < method.ϵ
    [first(v), last(v)]
  else
    v₁ = simplify(v[begin:imax], method)
    v₂ = simplify(v[imax:end],   method)
    [v₁[begin:end-1]; v₂]
  end
end


text:
min = 0

function getMin(arr)
#
  global min
  if length(arr) == 0
    return Nothing
  end

  min = arr[1] #
  for ele in arr
    if ele < min
      min = ele #
    end
  end

  return min
end

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
# ------- IPOPT ------------

"""
    IPOPT(options=Dict())

Use Ipopt as the optimizer

# Arguments
- `options::Dict`: options for Ipopt.  see Ipopt docs.
"""
struct IPOPT{TD} <: AbstractSolver
    options::TD
end

IPOPT() = IPOPT(Dict())

# wrapper for IPOPT
function optimize(solver::IPOPT, cache, x0, lx, ux, lg, ug, rows, cols)

    # initialize
    nx = length(x0)
    ng = length(lg)
    options = solver.options

    # initialize data for caching results since ipopt separates functions out
    xlast = 2*x0
    f = 0.0
    g = zeros(ng)
    df = zeros(nx)
    dg = zeros(length(rows))

    function ipcon(x, g)  # ipopt calls this function first

        f = evaluate!(g, df, dg, x, cache)
        xlast = x

        return nothing
    end

    function ipobj(x)
        if !isequal(x, xlast)
            f = func!(g, x)
        end
        return f
    end

    function ipgrad_obj(x, grad)
        if !isequal(x, xlast)
            f = evaluate!(g, df, dg, x, cache)
            xlast = x
        end
        grad[:] = df
        return nothing
    end

    function ipjac_con(x, mode, r, c, values)
        if mode == :Structure
            r[:] = rows
            c[:] = cols
        else
            if !isequal(x, xlast)
                f = evaluate!(g, df, dg, x, cache)
                xlast = x
            end
            values[:] = dg
        end
        return nothing
    end

    nzJ = length(rows)
    nH = 1  # irrelevant for quasi-newton
    prob = Ipopt.createProblem(nx, lx, ux, ng, lg, ug, nzJ, nH,
        ipobj, ipcon, ipgrad_obj, ipjac_con)

    # set options
    Ipopt.addOption(prob, "hessian_approximation", "limited-memory")
    for (key, value) in options
        Ipopt.addOption(prob, key, value)
    end

    # open output file
    filename = "ipopt.out"
    print_level = 5
    if haskey(options, "output_file")
        filename = options["output_file"]
    end
    if haskey(options, "print_level")
        print_level = options["print_level"]
    end
    Ipopt.openOutputFile(prob, filename, print_level)

    # solve problem
    prob.x = x0
    status = solveProblem(prob)

    return prob.x, prob.obj_val, Ipopt.ApplicationReturnStatus[status], nothing
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

# # Learning with a single neuron using Flux.jl
#
# In this notebook, we'll use `Flux` to create a single neuron and teach it to learn, as we did by hand in notebook 10!

#-

# ### Read in data and process it
#
# Let's start by reading in our data

using CSV, DataFrames, Flux

apples = DataFrame(CSV.File(datapath("data/apples.dat"), delim='\t', allowmissing=:none, normalizenames=true))
bananas = DataFrame(CSV.File(datapath("data/bananas.dat"), delim='\t', allowmissing=:none, normalizenames=true))

# and processing it to extract information about the red and green coloring in our images:

x_apples  = [ [row.red, row.green] for row in eachrow(apples)]
x_bananas = [ [row.red, row.green] for row in eachrow(bananas)]

xs = vcat(x_apples, x_bananas)

ys = vcat(fill(0, size(x_apples)), fill(1, size(x_bananas)))

# The input data is in `xs` and the labels (true classifications as bananas or apples) in `ys`.

#-

# ### Using `Flux.jl`

#-

# Now we can load `Flux` to really get going!

## using Pkg; Pkg.add("Flux")
using Flux

# We saw in the last notebook that σ is a built-in function in `Flux`.
#
# Another function that is used a lot in neural networks is called `ReLU`; in Julia, the function is called `relu`.

#-

# #### Exercise 1
#
# Use the docs to discover what `ReLU` is all about.
#
# `relu.([-3, 3])` returns
#
# A) [-3, 3]

# B) [0, 3]

# C) [0, 0]

# D) [3, 3]

#-

# #### Solution
#
#
# If you run
#
# ```julia
# ?relu
# ```
#
# you'll see that `relu` returns `0` for nonpositive input values and the input value itself for positive input values. Therefore `relu.([-3, 3])` will return (B) `[0, 3]`.

using Plots
plot(relu, -5, 5)

# ### Making a single neuron in Flux

#-

# Let's use `Flux` to build our neuron with 2 inputs and 1 output:
#
#

#-

# We previously put the two weights in a vector, $\mathbf{w}$. Flux instead puts weights in a $1 \times 2$ matrix (i.e. a matrix with 1 *row* and 2 *columns*).
#
# Previously, to compute the dot product of $\mathbf{w}$ and $\mathbf{x}$ we had to use either the `dot` function, or we had to transpose the vector $\mathbf{w}$:
#
# ```julia
# # transpose w
# b = w' * x
# # or use dot!
# b = dot(w, x)
# ```
# If the weights are instead stored in a $1 \times 2$ matrix, `W`, then we can simply multiply `W` and `x` together instead!
#
# We start off with random values for our parameters and data:

W = rand(1, 2)

#-

x = rand(2)

# Note that the product of `W` and `x` will now be an array (vector) with a single element, rather than a single number:

W * x

# This means that our bias `b` is treated as an array when we're using `Flux`:

b = rand(1)

# #### Exercise 2
#
# Write a function `mypredict` that will take a single input, array `x` and use `W`, `b`, and built-in `σ` to generate an output prediction (stored as an array). This function defines our neural network!
#
# Hint: This function will look very similar to $f_{\mathbf{w},\mathbf{b}}$ from the last notebook but has changed since our data structures to store our parameters have changed!

#-

# #### Solution

mypredict(x) = σ.(W*x + b)

# **Test**

W = rand(1, 2)
x = rand(2)
b = rand(1)

isapprox(mypredict(x), σ.(W*x + b))

# #### Exercise 3
#
# Define a loss function called `loss`.
#
# `loss` should take two inputs: a vector storing data, `x`, and a vector storing the correct "labels" for that data. `loss` should return the sum of the squares of differences between the predictions and the correct labels.

#-

# #### Solution

loss(x, y) = Flux.mse(mypredict(x), y)

# **Tests**

x, y = rand(2), rand(1)
isapprox( loss(x, y), sum((mypredict(x) .- y).^2) )

# ## Calculating gradients using Flux: backpropagation

#-

# For learning, we know that what we need is a way to calculate derivatives of the `loss` function with respect to the parameters `W` and `b`. So far, we have been doing that using finite differences.
#
# `Flux.jl` instead implements a numerical method called **backpropagation** that calculates gradients (essentially) exactly, in an automatic way, by indirectly applying the rules of calculus.
# To do so, it provides a new type of object called "tracked" arrays. These are arrays that store not only their current value, but also information about gradients, which is used by the backpropagation method.
#
# [If you want to understand the maths behind backpropagation, we recommend e.g. [this lecture](https://www.youtube.com/watch?v=i94OvYb6noo).]

#-

# To do so, `Flux` provides a function `param` to define such objects that will contain the information for a *param*eter.

#-

# Let's start, as usual, by setting up some random initial values for the parameters:

W_data = rand(1, 2)
b_data = rand(1)

W_data, b_data

# We now set up `Flux.jl` objects that will contain these values *and* their derivatives, and allow to propagate
# this information around:

W = param(W_data)
b = param(b_data)

# Here, `param` is a function that `Flux` provides to create an object that represents a parameter of a machine learning model, i.e. an object which has both a value and derivative information, and such that other objects know how to *keep track* of when it is used in an expression.

#-

# #### Exercise 4
#
# What type does `W` have?
#
# A) Array (1D)

# B) Array (2D)

# C) TrackedArray (1D)

# D) TrackedArray (2D)

# E) Parameter (1D)

# F) Parameter (2D)

#-

# #### Solution
#
# D) `TrackedArray` (2D)
#
# ```julia
# typeof(W)
# ```
# gives `TrackedArray{…,Array{Float64,2}}`. This can be interpreted as `TrackedArray{…,Array{T,N}}` where `T` is the type and `N` is the dimension.

#-

# #### Exercise 5
#
# `W` stores not only its current value, but also has space to store gradient information. You can access the values and gradient of the weights as follows:
#
# ```julia
# W.data
# W.grad
# ```
#
# At this point, are the values of the weights or the gradient of the weights larger?
#
# A) the values of the weights

# B) the gradient of the weights

#-

# #### Solution
#
# A) the values of the weights.
#
# We've randomly initialized the data in `W` with values between `0` and `1`. When we use `params` to create a `TrackedArray`, the gradient is initialized to `0`.

#-

# #### Exercise 6
#
# For data `x` and `y` where
#
# ```julia
# x, y = [0.413759, 0.692204], [0.845677]
# ```
# apply the loss function to `x` and `y` to give a new variable `l`. What is the type of `l`? (How many dimensions does it have?)
#
# A) Array (0D)

# B) Array (1D)

# C) TrackedArray (0D)

# D) TrackedArray (1D)

# E) Float64

# F) Int64

#-

# #### Solution
#
# C) `TrackedArray` (0D) (Note that this means that `l` is just a single number!)

mypredict(x) = σ.(W*x + b)
loss(x, y) = sum((mypredict(x) .- y).^2)
x, y = [0.413759, 0.692204], [0.845677]
l = loss(x, y)
typeof(l)

# Having set up the structure, we can now **propagate information about derivatives backwards ** from the `loss` function to all of the objects that are used to calculate it:

using Flux.Tracker

back!(l)   # backpropagate derivatives of the loss function

# and now we can look at the derivatives again:

W.grad

#-

b.grad

# What are these results? They are the components of the **gradient of the loss function** with respect to each component of the object `W`, and with respect to `b`! So as promised, `Flux` has done the hard work of calculating derivatives for us!
#
# *Bonus info*:
#
# To do so, internally Flux sets up a "computational graph" and propagates the information on derivatives backwards through the graph. Each node of the graph knows which objects feed into that node, so it tells them to also update their gradients, etc.

#-

# ### Stochastic gradient descent

#-

# We can now use these features to reimplement stochastic gradient descent, following the method we used in the previous notebook, but now using backpropagation!

#-

# #### Exercise 7
#
# Modify the code from the previous notebook for stochastic gradient descent to use Flux instead.

#-

# #### Solution

x, y = [0.413759, 0.692204], [0.845677]
loss(x, y)

#-

function stochastic_gradient_descent(loss, W, b, xs, ys, N=1000)

    η = 0.01

    for step in 1:N

        i = rand(1:length(xs))  # choose a data point

        x = xs[i]
        y = ys[i]

        l = loss(x, y)
        back!(l)
        b.data .-= η * b.grad
        W.data .-= η * W.grad
    end

    return W, b

end

# ### Investigating stochastic gradient descent

#-

# Let's look at the values stored in `b` before we run stochastic gradient descent:

b

# After running `stochastic_gradient_descent`, we find the following:

W_final, b_final = stochastic_gradient_descent(loss, W, b, xs, ys, 1000)

# we can look at the values of `W_final` and `b_final`, which our machine learned to generate our desired classification.

W_final

#-

b_final

# #### Exercise 8
#
# Plot the data and the learned function.

#-

# #### Solution

## using Pkg; Pkg.add("Plots")
using Plots; gr()

# Let's draw the function that the network has learned, together with the data:

heatmap(0:0.01:1, 0:0.01:1, (x,y)->mypredict([x, y]).data[1])

scatter!(first.(x_apples), last.(x_apples), m=:cross)
scatter!(first.(x_bananas), last.(x_bananas))

# #### Exercise 9
#
# Do this plot every so often as the learning process is proceeding in order to have an animation of the process.

#-

# ### Automation with Flux.jl

#-

# We will need to repeat the above process for a lot of different systems.
# Fortunately, `Flux.jl` provides us with tools to automate this!

#-

# Flux allows to create a neuron in a simple way:

## using Pkg; Pkg.add("Flux")
using Flux

#-

model = Dense(2, 1, σ)

# The `2` and `1` refer to the number of inputs and outputs, and the neuron is defined using the $\sigma$ function.

typeof(model)

# We have made an object of type `Dense`, defined by `Flux`, with the name `model`. This represents a "dense neural network layer" (see later for more on neural network layers).
# The parameters that will be modified during the learning process live *inside* the `model` object.

#-

# #### Exercise 10
#
# Investigate which variables live inside the `model` object and what type they are. How does that compare to the call to create the `Dense` object that we started with?

#-

# #### Solution

#-

# Use `model.<TAB>` (i.e. write `model.` and then press the `TAB` key to do autocompletion) to check interactively what is inside the `model` object, or

model.W

#-

model.b

# The fact that `model.W` and `model.b` are of size $1 \times 2$ and $1$, respectively, comes from the `(2, 1)` pair in the call to the `Dense` constructor when we created `model`.
#
# `model.W` will be multiplied by a vector `x` of length 2, which it is why it needs to be of size $1 \times 2$.
#
# Again, these are tracked arrays so that Flux can calculate their gradients.

#-

# ### Model object as a function

#-

# We can apply the `model` object to data just as if it were a standard function:

model(rand(2))

# #### Exercise 11
#
# Prove to yourself that you understand what is going on when we call `model`. Create two arrays `W` and `b` with the same elements as `model.W` and `model.b`. Use `W` and `b` to generate the same answer that you get when we call `model([.5, .5])`.

#-

# ### Using Flux

#-

# We now need to provide Flux with three pieces of information:
#
# 1. A loss function
# 2. Some training data
# 3. An optimization method

#-

# ### Data

#-

# The data can take a couple of different forms.
# One form is a single **iterator**, consisting of pairs $(x, y)$ of data and labels.
# We can achieve this with `zip`.

#-

# #### Exercise 12
#
# Use `zip` to "zip together" `xs` and `ys`. Then use the `collect` function to check what `zip` actually does.

#-

# #### Solution

data = zip(xs, ys)

#-

collect(data)

# ### Optimization routine

#-

# Now we need to tell Flux what kind of optimization routine to use. It has several built in; the standard stochastic gradient descent algorithm that we have been using is called `SGD`. We must pass it two things: a list of parameter objects which will be modified by the optimization routine, and a step size:

opt = SGD([model.W, model.b], 0.01)
## give a list of the parameters that will be modified

# The gradient calculations and parameter updates will be carried out by this optimizer function; we do not see those details, but if you are curious, you can, of course, look at the `Flux.jl` source code!

#-

# ### Training

#-

# We now have all the pieces in place to actually **train** our model (a single neuron) on the data.
# "Training" refers to using pre-labeled data to learn the function that relates the input data to the desired output data given by the labels.
#
# `Flux` provides the function `train!`, which performs a single pass through the data and does a single step of optimization using the partial cost function for each data point:

Flux.train!(loss, data, opt)

# We can then just repeat this several times to train the network more and coax it towards the minimum of the cost function:

for i in 1:100
    Flux.train!(loss, data, opt)
end

# Now let's look at the parameters after training:

model.W

#-

model.b

# Instead of writing out a list of parameters to modify, `Flux` provides the function `params`, which extracts all available parameters from a model:

opt = SGD(params(model), 0.01)

#-

params(model)

# ## Adding more features

#-

# #### Exercise 13
#
# So far we have just used two features, red and green.
#
# (i) Add a third feature, blue. Plot the new data.
#
# (ii) Train a neuron with 3 inputs and 1 output on the data.
#
# (iii) Can you find a good way to visualize the result?

max_stars_repo_path: Courses/Deep learning with Flux/1300.Intro-to-Flux.jl-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
#. PREAMBLE

#= This file defines many term orders. A term order is basically a way to
transform an exponent vector into an "order vector". The ordering is then the
lexicographic ordering on the order vector. The specific application of symbolic
integration requires very baroque term order. Thanks to Julia's JIT compilation,
even baroque monomial order can be compiled rather efficiently.

The type of the exponents is not made explicit, we can choose any signed integer type.
=#

using StaticArrays

#. Some helpers for SVectors

"""
    srange(v :: StaticVector, ::Val{first}, ::Val{last})

Extract consecutive coefficients from a static vector.
"""
@generated function srange(v :: StaticVector, ::Val{first}, ::Val{last}) where {first, last}
    :(SVector( $([:(v[$i]) for i in first:last]...) ))
end

@generated function _concat(::Size{p}, va :: StaticVector, ::Size{q}, vb :: StaticVector) where {p,q}
    :(SVector{p[1]+q[1]}($([:(va[$i]) for i in 1:p[1]]...), $([:(vb[$i]) for i in 1:q[1]]...)))
end

"""
    concat(va :: StaticVector, vb :: StaticVector)

Concatenate two static vectors into a static vector.
"""
concat(va, vb) = _concat(Size(va), va, Size(vb), vb)

#. Abstract interface

""" Abstract type for term order on monomials with N exponents. Includes term
order on modules (where one or several exponents are used to indicate the
position).

Interface for T <: TermOrder{N}:
    - `ordervector(::T, ::SVector) :: SVector`
        return a static vector on which to perform lexicographic ordering.

    - iscompatible(::T, e::SVector, f::SVector)
        indicates if `e`  and `f` have a common multiple. Only useful for module orders.
"""
abstract type TermOrder{N} end

nvars(t::TermOrder{N}) where N = N

"""
    divides(t::TermOrder{N}, e::SVector{N}, f::SVector{N})

Returns true if and only if the two terms are compatible (that is, in the case of a module order, they live at the same position)
and `f` is a multiple of `e`.
"""
divides(t::TermOrder{N}, e::SVector{N}, f::SVector{N}) where N =
    iscompatible(t, e, f) && all(e .<= f)

""" Abstract type for monomial order in a ring. """
abstract type MonomialOrder{N} <: TermOrder{N} end

iscompatible(::MonomialOrder{N}, e::SVector{N}, f::SVector{N}) where N = true

#. Concrete types

#.. Lexicographic order
struct Lex{N} <: MonomialOrder{N} end
Lex(N) = Lex{N}()

ordervector(::Lex{N}, e::SVector{N}) where N = e

#.. Graded reverse lexicographic order
struct Grevlex{N} <: MonomialOrder{N} end
Grevlex(N) = Grevlex{N}()

# should be compiled statically
ordervector(::Grevlex{N}, e::SVector{N}) where N = insert(deleteat(reverse(-e), N), 1, sum(e))

#.. Weighted grevlex
struct WGrevlex{N} <: MonomialOrder{N}
    weights :: NTuple{N, Int}
end

ordervector(o::WGrevlex{N}, e::SVector{N}) where N = ordervector(Grevlex(N), o.weights .* e)

#.. Block order

""" Lexicographic on two blocks """
struct Block{N, O1 <: TermOrder, O2 <: TermOrder} <: TermOrder{N}
    o1 :: O1
    o2 :: O2
end

Block(o1, o2) = Block{nvars(o1)+nvars(o2), typeof(o1), typeof(o2)}(o1, o2)

Base.:*(o1::TermOrder, o2::TermOrder) = Block(o1, o2)

split(b::Block{N}, e::SVector{N}) where {N} =
    (srange(e, Val(1), Val(nvars(b.o1))), srange(e, Val(nvars(b.o1)+1), Val(N)))

ordervector(b::Block{N}, e::SVector{N}) where N =
    concat(ordervector(b.o1, split(b, e)[1]), ordervector(b.o2, split(b, e)[2]))

iscompatible(b::Block{N}, e::SVector{N}, f::SVector{N}) where N =
    iscompatible(b.o1, split(b, e)[1], split(b, f)[1]) &&
    iscompatible(b.o2, split(b, e)[2], split(b, f)[2])

#. Sort interface

@generated function lt(t::TermOrder{N}, e::SVector{N}, f::SVector{N}) where N
    quote
        oe = ordervector(t, e)
        of = ordervector(t, f)

        $([:(oe[$i] < of[$i] && return true ; oe[$i] > of[$i] && return false) for i in 1:N]...)

        return false
    end
end

# @generated function divides(t::TermOrder{N}, e::SVector{N}, f::SVector{N}) where N
#     quote
#         $([:(e[$i] > f[$i] && return false) for i in 1:N]...)

#         return true
#     end
# end


text:
a = rand(Float32, 5, 5)

function test(x, inc)
    return x + inc
end

function invalidate(x, threshold)::Bool
    x > threshold ? true : false
end

b = test.(a, 2)
c = map(x -> test(x, 2), a)

println(c)

d = invalidate.(a, 0.9)
e = map(x -> invalidate(x, 0.9), a)

println(a)
println(d)
println(e)

function invalidate2(x, datamin, datamax, ignrval)::Bool
    val = Float32(x)

    !isfinite(val) || (val < datamin) || (val > datamax) || (val <= ignrval)
end

f = map(x -> invalidate2(x, 0.3, 0.7, 0.0), a)

println(f)


text:
## Example:
## DGP - Bivariate VAR(1) Model from Kilian, RESTAT, 1998
# B11 set to 0.01
using VectorAutoregressions
using Plots
plotly()
using Random, LinearAlgebra
using Statistics, GrowableArrays

const T,K = 1000,2
const H = 24
const nrep = 1000
const p = 1
const intercept = false
const burnin = 100

# Random.seed!(1234)
mR  = [0.9 0.01; 0.5 0.5]             #  Coefficients of the 1st lag B11 set to 0.01
mΣ_true = [1 0.3; 0.3 1]              #  vcv of u
#mP = cholfact(mΣ_true)[:L]            #  Cholesky PP'= Σ
mP = cholesky(mΣ_true, Val(false)).L
mI = Matrix(1.0I,K,K)
mIRF_true = (mP'*mI)[:]       # Initialize structural IRF
# True IRF
for h=1:H
    vIRF = (mP*mR^h)[:]
    global  mIRF_true = hcat(mIRF_true,vIRF)
end

# Generate data from the VAR(1) DGP
y = zeros(T+burnin,K)
gen_var1_data!(y,mR,mP,burnin)

# Fit VAR(1) to data
V = VAR(y,p,intercept)
mIRFb = IRFs_b(V,H,nrep,intercept)
mIRFa = IRFs_a(V,H,intercept)

# Plot true irf vs model irf + asy ci
pIRF_asy = plot(layout = grid(K,K));
[plot!(pIRF_asy, [mIRF_true[i,:] mIRFa.IRF[i,:] mIRFa.CI.CIl[i,:] mIRFa.CI.CIh[i,:]], color = ["blue" "red" "red" "red"],
line = [:solid :solid :dash :dash], legend = false, label = ["true" "model"], subplot = i) for i in 1:K^2]
gui(pIRF_asy)

# Plot true irf vs model irf + bootstraped ci
pIRF_boot = plot(layout = grid(K,K));
[plot!(pIRF_boot, [mIRF_true[i,:] mIRFb.CI.CIl[i,:] mIRFb.IRF[i,:] mIRFb.CI.CIh[i,:]], color = ["blue" "black" "black" "black"],
line = [:solid :solid :dash :dash], legend = false, label = ["true" "model"], subplot = i) for i in 1:K^2]
gui(pIRF_boot)

max_stars_repo_name: lucabrugnolini/VAR

max_issues_repo_name: lucabrugnolini/VAR

max_forks_repo_name: lucabrugnolini/VAR


text:
using Random, Distributions, Turing, Plots

σ_true = 1 # true standard deviation of measurement device
temperature = 21.0 # true temperature

# we take 100 measurements of the temperature
measurements = rand(Normal(temperature, σ_true), 100)

histogram(measurements)

#=
prior on temperature: I'm about 95% certain that the temperature will
lie between -40 and +40, i.e. 95% of the distribution lies between μ-2σ and  μ+2σ
=#
temperature_prior = Normal(0, 20)
plot(x -> pdf(temperature_prior, x), -100, 100)
xlabel!("Prior belief about temperature")

# inference model
@model infer_temperture(y) = begin
    N = length(y)

    # prior on measurement noise
    σ ~ Truncated(Cauchy(0, 2), 0, Inf)
    # prior on temperature
    μ ~ Normal(0, 20)

    # observations
    for i ∈ 1:N
        y[i] ~ Normal(μ, σ)
    end
    return μ, σ
end

# P(μ | y) ∝ P(y | μ) * P(μ)
# algorithm1 = NUTS(1000,  0.65)
algorithm = SMC(3000)

temp_samples = sample(infer_temperture(measurements), algorithm)

histogram(temp_samples[:σ])
plot(temp_samples[:σ])
plot(temp_samples[:μ])

mean(temp_samples[:μ])
mean(temp_samples[:σ])

max_stars_repo_name: kogpsy/vpompd

max_issues_repo_name: kogpsy/vpompd

max_forks_repo_name: kogpsy/vpompd


text:
using Flux, LinearAlgebra, Polynomials, Plots

const n = 100
# target polynomial
c = ChebyshevT([-1,0,-2,0,1,0,1,2,3])
target = convert(Polynomial, c)
plot(target, (-1.,1.)...,label="target")

length(target.coeffs)

# mapping polynomial to circulant matrix
S = circshift(Matrix{Float64}(I, n, n),(1,0))
param = zeros(n)
param[1:9] = target.coeffs
Circulant = param
for k in 1:n-1
    Circulant = hcat(Circulant, S^k*param)
end

# creating dataset
bs = 3000
x = randn(Float32,n,1,bs)
y = convert(Array{Float32},reshape(transpose(Circulant)*dropdims(x;dims=2) ,(n,1,bs)))
data = [(x,y)]

# padding function to work modulo n
function pad_cycl(x;l=1,r=1)
    last = size(x,1)
    xl = selectdim(x,1,last-l+1:last)
    xr = selectdim(x,1,1:r)
    cat(xl, x, xr, dims=1)
end

# neural network with 7 convolution layers
model = Chain(
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros()),
    x -> pad_cycl(x,l=0,r=2),
    CrossCor((3,),1=>1,bias=Flux.Zeros())
)

loss(x, y) = Flux.Losses.mse(model(x), y)
ps = Flux.params(model)
loss_vector = Vector{Float32}()
logging_loss() = push!(loss_vector, loss(x, y))
opt = ADAM(0.2)
n_epochs = 1700
for epochs in 1:n_epochs
    Flux.train!(loss, ps, data, opt, cb=logging_loss)
    if epochs % 50 == 0
        println("Epoch: ", epochs, " | Loss: ", loss(x,y))
    end
end

plot(loss_vector)

pred = Polynomial([1])
for p in ps
    if typeof(p) <: Array
        pred *= Polynomial([p...])
    end
end
plot(target, (-1.,1.)...,label="target")
ylims!((-10,10))
plot!(pred, (-1.,1.)...,label="pred")


text:
# This file was generated, do not modify it. # hide
nnregressor.batch_size = 2
MLJ.fit!(mach, rows=train, verbosity=3)


text:
import PeriodicGraphs

function export_dataline(f, x)
    inbetween = '\t'
    if any(isspace, x)
        if any(==('\n'), x)
            x = ';'*x*"\n;"
            inbetween = '\n'
        else
            @toggleassert x[1]!= '\''
            x = '\''*x*'\''
        end
    end
    println(f, inbetween*x)
end

_representative_atom(ty, i) = last(representative_atom(ty, i))
function PeriodicGraphEmbeddings.export_vtf(file, c::Union{CrystalNet,Crystal}, repeatedges=6, colorname=false)
    if c isa CrystalNet
        net = CrystalNet3D(c)
        return export_vtf(file, net.pge, net.types, repeatedges, colorname, string_atomtype, _representative_atom)
    end
    return export_vtf(file, c.pge, c.types, repeatedges, colorname, string_atomtype, _representative_atom)
end

function export_cif(file, __c::Union{Crystal, CIF})
    mkpath(splitdir(file)[1])
    c, info = if __c isa CIF
        ___c = expand_symmetry(__c)
        ___c, copy(__c.cifinfo)
    else
        __c, Dict{String,Union{String,Vector{String}}}("data"=>last(splitdir(file)))
    end
    loops = Dict{Int, Tuple{Vector{String}, Vector{Vector{String}}}}()
    open(file, write=true) do f
        print(f, "data_"); println(f, popstring!(info, "data")*"\n\n")
        for (id, data) in info
            if data isa String
                print(f, '_'); print(f, id)
                export_dataline(f, data)
            else
                @toggleassert data isa Vector{String}
                n = length(data)
                if haskey(loops, n)
                    push!(loops[n][1], id)
                    push!(loops[n][2], data)
                else
                    loops[n] = (String[id], Vector{String}[data])
                end
            end
        end

        cell::Cell{Rational{Int}} = c isa CIF ? c.cell : c.pge.cell
        #scale_factor::Float64 = unique!(sort(c.types)) == [:Si] && c.pge.pos[:,1] != [0,0,0] ? 1.2 : 1.0
        ((__a, __b, __c), (__α, __β, __γ)), _ = cell_parameters(cell)
        _a, _b, _c, _α, _β, _γ = Float64.((__a, __b, __c, __α, __β, __γ))

        println(f, """

        _cell_length_a\t\t$_a
        _cell_length_b\t\t$_b
        _cell_length_c\t\t$_c
        _cell_angle_alpha\t\t$_α
        _cell_angle_beta\t\t$_β
        _cell_angle_gamma\t\t$_γ

        _symmetry_space_group_name_H-M\t'P 1'
        _symmetry_Int_Tables_number\t1
        _symmetry_cell_setting\ttriclinic
        """)

        println(f, """
        loop_
        _symmetry_equiv_pos_as_xyz
        x,y,z""")

        for eq in cell.equivalents
            println(f, eq)
        end
        println(f)

        n::Int = c isa CIF ? length(c.ids) : length(c.types)
        if !haskey(loops, n)
            loops[n] = (String[], Vector{String}[])
        end
        append!(loops[n][1],
                ["atom_site_label", "atom_site_type_symbol",
                 "atom_site_fract_x", "atom_site_fract_y", "atom_site_fract_z"])
        labels = String[string_atomtype(c.types[c isa CIF ? c.ids[i] : i]::Symbol)*string(i) for i in 1:n]
        pos = string.(round.(c isa CIF ? c.pos::Matrix{Float64} : reduce(hcat, c.pge.pos)::Matrix{Float64}; sigdigits=6))
        append!(loops[n][2], [labels, string_atomtype.((c.types[x])::Symbol for x in (c isa CIF ? c.ids : collect(1:n))::Vector{Int}),
                pos[1,:], pos[2,:], pos[3,:]])

        if c isa Union{Crystal{Clusters}, Crystal{Nothing}}
            bonds = edges((c isa CIF ? c.graph : c.pge.g)::PeriodicGraph3D)
            src = String[]
            dst = String[]
            mult = Int[]
            last_bond = (0, 0)
            n = 0
            for e in bonds
                if (e.src, e.dst.v) == last_bond
                    mult[end] += 1
                else
                    push!(src, labels[e.src])
                    push!(dst, labels[e.dst.v])
                    push!(mult, 1)
                    last_bond = (e.src, e.dst.v)
                    n += 1
                end
            end
            if !haskey(loops, n)
                loops[n] = (String[], Vector{String}[])
            end
            append!(loops[n][1], ["geom_bond_atom_site_label_1", "geom_bond_atom_site_label_2", "geom_bond_multiplicity"])
            append!(loops[n][2], [src, dst, string.(mult)])
        end

        for (n, (ids, datas)) in loops
            println(f, "loop_")
            for id in ids
                print(f, '_'); println(f, id)
            end
            for i in 1:n
                for j in 1:length(datas)
                    print(f, datas[j][i]); print(f, '\t')
                end
                println(f)
            end
            println(f)
        end
        println(f, "\n#END")
    end
    nothing
end

PeriodicGraphEmbeddings.export_cgd(file, c::Crystal) = export_cgd(file, c.pge)
PeriodicGraphEmbeddings.export_cgd(file, c::CrystalNet) = export_cgd(file, PeriodicGraph3D(c.pge.g))

function export_attributions(crystal::Crystal{Clusters}, path=joinpath(tempdir(),tempname()))
    frame = Chemfiles.Frame()
    m = length(crystal.clusters.classes)
    # residues = [Chemfiles.Residue(string(i)) for i in 1:m]
    # for i in 1:m
    #     Chemfiles.set_property!(residues[i], "chainname", string(i))
    #     Chemfiles.set_property!(residues[i], "chainid", string(i))
    #     Chemfiles.add_residue!(frame, residues[i])
    # end
    ((_a, _b, _c), (_α, _β, _γ)), mat = cell_parameters(c.pge.cell)
    Chemfiles.set_cell!(frame, Chemfiles.UnitCell(Float64[_a, _b, _c], Float64[_α, _β, _γ]))
    recenter::SVector{3,Float64} = minimum(reduce(hcat, crystal.pge.pos); dims=2)
    for i in 1:length(crystal.types)
        # resid = crystal.clusters.attributions[i]
        atom = Chemfiles.Atom(string(crystal.clusters.classes[crystal.clusters.attributions[i]]))
        Chemfiles.set_type!(atom, string_atomtype(crystal.types[i]))
        Chemfiles.add_atom!(frame, atom, collect(Float64.(mat * (crystal.pge.pos[i] .- recenter))))
        # Chemfiles.add_atom!(residues[resid], i)
    end
    for e in edges(crystal.pge.g)
        iszero(last(last(e))) || continue
        Chemfiles.add_bond!(frame, first(e)-1, first(last(e))-1)
    end
    target = isempty(last(splitext(path))) ? path*".pdb" : path
    output = Chemfiles.Trajectory(target, 'w')
    write(output, frame)
    close(output)
end

"""
    export_arc(path, arc=CRYSTAL_NETS_ARCHIVE)

Export archive `arc` to the specified `path`. If unspecified, the exported archive is the
current one.
"""
function export_arc(path, arc=CRYSTAL_NETS_ARCHIVE)
    mkpath(splitdir(path)[1])
    open(path, "w") do f
        println(f, "Made by CrystalNets.jl v$CRYSTAL_NETS_VERSION")
        if !isnothing(arc)
            println(f)
            pairs = sort(collect(arc); by=last)
            for (genome, id) in pairs
                println(f, "key\t", genome)
                println(f, "id\t", id)
                println(f)
            end
        end
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    ln_mean(x, y)

Compute the logarithmic mean

    ln_mean(x, y) = (y - x) / (log(y) - log(x)) = (y - x) / log(y / x)

Problem: The formula above has a removable singularity at `x == y`. Thus,
some care must be taken to implement it correctly without problems or loss
of accuracy when `x ≈ y`. Here, we use the approach proposed by
Ismail and Roe (2009).
Set ξ = y / x. Then, we have

    (y - x) / log(y / x) = (x + y) / log(ξ) * (ξ - 1) / (ξ + 1)

Set f = (ξ - 1) / (ξ + 1) = (y - x) / (x + y). Then, we use the expansion

    log(ξ) = 2 * f * (1 + f^2 / 3 + f^4 / 5 + f^6 / 7) + O(ξ^9)

Inserting the first few terms of this expansion yields

    (y - x) / log(ξ) ≈ (x + y) * f / (2 * f * (1 + f^2 / 3 + f^4 / 5 + f^6 / 7))
                     = (x + y) / (2 + 2/3 * f^2 + 2/5 * f^4 + 2/7 * f^6)

Since divisions are usually more expensive on modern hardware than
multiplications (Agner Fog), we try to avoid computing two divisions. Thus,
we use

    f^2 = (y - x)^2 / (x + y)^2
        = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y)

Given ε = 1.0e-4, we use the following algorithm.

    if f^2 < ε
      # use the expansion above
    else
      # use the direct formula (y - x) / log(y / x)
    end

# References
- Ismail, Roe (2009).
  Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks.
  [DOI: 10.1016/j.jcp.2009.04.021](https://doi.org/10.1016/j.jcp.2009.04.021)
- Agner Fog.
  Lists of instruction latencies, throughputs and micro-operation breakdowns
  for Intel, AMD, and VIA CPUs.
  https://www.agner.org/optimize/instruction_tables.pdf
"""
@inline function ln_mean(x, y)
  epsilon_f2 = 1.0e-4
  f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2
  if f2 < epsilon_f2
    return (x + y) / @evalpoly(f2, 2, 2/3, 2/5, 2/7)
  else
    return (y - x) / log(y / x)
  end
end

"""
    inv_ln_mean(x, y)

Compute the inverse `1 / ln_mean(x, y)` of the logarithmic mean
[`ln_mean`](@ref).

This function may be used to increase performance where the inverse of the
logarithmic mean is needed, by replacing a (slow) division by a (fast)
multiplication.
"""
@inline function inv_ln_mean(x, y)
  epsilon_f2 = 1.0e-4
  f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2
  if f2 < epsilon_f2
    return @evalpoly(f2, 2, 2/3, 2/5, 2/7) / (x + y)
  else
    return log(y / x) / (y - x)
  end
end

end # @muladd


text:
interpolation_keywords = ["displacement","epsilon","pressure_interpolated","sigma","temperature_interpolated","\"displacement\"","\"epsilon\"","\"pressure_interpolated\"","\"sigma\"","\"temperature_interpolated\""]
uncompress_keywords = ["connectivity","offsets","bulk_node_ids","bulk_element_ids","Points","MaterialIDs","\"connectivity\"","\"offsets\"","\"bulk_node_ids\"","\"bulk_element_ids\"","\"Points\"","\"MaterialIDs\"","types","\"types\""]

"""
    VTUKeyWords

Container for interpolation and uncompress keywords. Will be created as a module-global constant value `vtukeywords`.

# Fields
- `interpolation_keywords::Vector{String}` : Vector with names of VTU fields onto which the math operators should apply
- `uncompress_keywords` : Vector with names of VTU field that should be uncompressed and stored in memory

# Standard values

By default, keywords were set that match results of the simulation software [OpenGeoSys](https://www.opengeosys.org/).

```julia
interpolation_keywords = ["displacement","epsilon","pressure_interpolated","sigma","temperature_interpolated"]
uncompress_keywords = ["connectivity","offsets","bulk_node_ids","bulk_element_ids","Points","MaterialIDs"]
```

# Example

If a VTU result file e.g. owns a field named `field1` to which the math operators should be applied, it can be added via:
```julia
add_interpolation_keywords(["field1"]])
```
"""
mutable struct VTUKeyWords
	interpolation_keywords::Vector{String}
	uncompress_keywords::Vector{String}
end

const vtukeywords = VTUKeyWords(interpolation_keywords,uncompress_keywords)

"""
	set_uncompress_keywords(uk::Vector{String})

Sets all fields that should be accessible and therefore needs to be uncompressed

# Arguments
- `uk::Vector{String}`: fieldnames
"""
function set_uncompress_keywords(uk::Vector{String})
	vtukeywords.uncompress_keywords = uk
	for k in copy(uk)
		push!(vtukeywords.uncompress_keywords, "\""*k*"\"") #hack for downwards compatibility with older ogs6-versions
	end
	return nothing
end

"""
	add_uncompress_keywords(uk::Vector{String})

Adds fields that should be accessible and therefore needs to be uncompressed

# Arguments
- `uk::Vector{String}`: fieldnames
"""
function add_uncompress_keywords(uk::Vector{String})
	for k in copy(uk)
		push!(vtukeywords.uncompress_keywords, "\""*k*"\"") #hack for downwards compatibility with older ogs6-versions
	end
	for k in copy(uk)
		push!(vtukeywords.uncompress_keywords, k) #hack for downwards compatibility with older ogs6-versions
	end
	return nothing
end

"""
	set_interpolation_keywords(ik::Vector{String})

Sets all fields onto which the math operators should be applied

# Arguments
- `ik::Vector{String}`: fieldnames
"""
function set_interpolation_keywords(ik::Vector{String})
	vtukeywords.interpolation_keywords = ik
	for k in copy(ik)
		push!(vtukeywords.interpolation_keywords, "\""*k*"\"") #hack for downwards compatibility with older ogs6-versions
	end
	return nothing
end

"""
	add_interpolation_keywords(ik::Vector{String})

Adds fields onto which the math operators should be applied

# Arguments
- `ik::Vector{String}`: fieldnames
"""
function add_interpolation_keywords(ik::Vector{String})
	for k in copy(ik)
		push!(vtukeywords.interpolation_keywords, "\""*k*"\"") #hack for downwards compatibility with older ogs6-versions
	end
	for k in copy(ik)
		push!(vtukeywords.interpolation_keywords, k) #hack for downwards compatibility with older ogs6-versions
	end
	return nothing
end


text:
using GeoFormatTypes, Test
using GeoFormatTypes: Geom, CRS, Mixed

@testset "Test construcors" begin
    @test_throws ArgumentError ProjString("+lat_ts=56.5 +ellps=GRS80")
    @test_throws ArgumentError EPSG("ERROR:4326")
    @test EPSG("EPSG:4326") == EPSG(4326)
end

@testset "Test conversion to string or int" begin
    @test convert(String, ProjString("+proj=test")) == "+proj=test"
    @test convert(String, EPSG(4326)) == "EPSG:4326"
    @test convert(Int, EPSG(4326)) == 4326
    @test convert(String, WellKnownText("test")) == "test"
    @test convert(String, WellKnownText2("test")) == "test"
    @test convert(String, GML("test")) == "test"
    @test convert(String, KML("test")) == "test"
    @test convert(String, GeoJSON("test")) == "test"
end

# `convert` placeholder methods
Base.convert(target::Type{<:GeoFormat}, mode::Union{Geom,Type{Geom}}, source::GeoFormat; kwargs...) =
    (:geom, kwargs...)
Base.convert(target::Type{<:GeoFormat}, mode::Union{CRS,Type{CRS}}, source::GeoFormat; kwargs...) =
    (:crs, kwargs...)

@testset "Test convert mode allocation" begin
    @testset "Test identical type is passed through unchanged" begin
        @test convert(WellKnownText, WellKnownText(Mixed(), "test")) == WellKnownText(Mixed(), "test")
        @test convert(ProjString, ProjString("+proj=test")) == ProjString("+proj=test")
    end
    @testset "Test conversions are assigned to crs or geom correctly" begin
        @test convert(WellKnownText, WellKnownText2(CRS(), "test")) == (:crs,)
        @test convert(WellKnownText2, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(WellKnownBinary, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(ProjString, WellKnownText(CRS(), "test")) == (:crs,)
        @test convert(EPSG, ProjString("+proj=test")) == (:crs,)
        @test convert(CoordSys, ProjString("+proj=test")) == (:crs,)

        @test convert(GeoJSON, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(KML, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(GML, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(ESRIWellKnownText, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownBinary, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Geom(), "test")) == (:geom,)
        @test convert(WellKnownText, WellKnownText2(Geom(), "test")) == (:geom,)

        @test convert(GeoJSON, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(KML, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(GML, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(ESRIWellKnownText, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(WellKnownBinary, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(WellKnownText2, WellKnownText(Mixed(), "test")) == (:geom,)
        @test convert(WellKnownText, WellKnownText2(Mixed(), "test")) == (:geom,)
    end
    @testset "Test kargs pass through convert" begin
        @test convert(WellKnownText, WellKnownText2(CRS(), "test"); order=:trad) == (:crs, :order=>:trad,)
        @test convert(GML, WellKnownText(Mixed(), "test"); order=:custom) == (:geom, :order=>:custom)
    end
    @testset "Test conversions that are not possible throw an error" begin
        @test_throws ArgumentError convert(KML, ProjString("+proj=test"))
        @test_throws ArgumentError convert(GeoJSON, ProjString("+proj=test"))
        @test_throws ArgumentError convert(ProjString, WellKnownText(Geom(), "test"))
        @test_throws ArgumentError convert(CoordSys, WellKnownText(Geom(), "test"))
        @test_throws ArgumentError convert(EPSG, WellKnownText(Geom(), "test"))
    end
end


text:
using DataFrames, Statistics, LinearAlgebra
using Random, Distributions, Plots, StatsBase

ProjDir = @__DIR__
#Random.seed!(12591)

include(joinpath(ProjDir,
    "./optimal_price_probabilities/optimal_price_probabilities.jl"))

# parameters
prices = [1.99, 2.49, 2.99, 3.49, 3.99, 4.49]

# parameters of the true (unknown) demand model
true_intc = 50
true_slope = 7

# prior distribution for each price
θ = [DataFrame() for i in 1:length(prices)]
for (i,p) in enumerate(prices)
    append!(θ[i],
        DataFrame(
            :price => p,
            :α => 30.0,
            :β => 1.0,
            :mean => 30.0
        )
    )
end

demands = true_intc .- true_slope * prices
probs = optimal_price_probabilities(prices, demands, 60)

α_demand = 50
β_demand = 7

# distribution for each price
θ = DataFrame()
for p in prices
    append!(θ,
        DataFrame(
            :price => p,
            :α => 30.0,
            :β => 1.0,
            :mean => 30.0
        )
    )
end

function sample_demands_from_model(theta)
    [rand(Gamma(v[:α], 1/v[:β]), 1)[1] for v in eachrow(θ)]
end

function sample_demand(price, demand_a=α_demand, demand_b=β_demand)
    demand = demand_a - demand_b * price
    rand(Poisson(demand), 1)[1]
end

function simulate(;T=T, verbose=false, report_final=true)
    # history for each price
    history = DataFrame()

    for t in 1:T
        demands = sample_demands_from_model(θ)
        price_probs = optimal_price_probabilities(prices, demands, 60)
        price_index_t = sample(1:length(prices), price_probs)
        price_t = prices[price_index_t]
        demand_t = sample_demand(price_t)
        revenue = round(demand_t * price_t, digits=2)

        θ[price_index_t, :α] += demand_t
        θ[price_index_t, :β] += 1
        θ[price_index_t, :mean] = θ[price_index_t, :α] / θ[price_index_t, :β]

        if verbose || (report_final && t == T)
            println("\nPrices = $prices")
            println("Demands (step=$t) = $(round.(demands, digits=2))")
            println("Price probs (step=$t) = $(price_probs)")
            println("Price index (step=$t) = $(price_index_t)")
            print("Selected price (step=$t) = $price_t =>")
            println(" demand = $demand_t, revenue = $revenue")
        end

        v = θ[price_index_t, :]
        append!(history,
            DataFrame(
                :step => t,
                :index => price_index_t,
                :price => v[:price],
                :α => v[:α],
                :β => v[:β]
            )
        )
    end

    (θ, history)
end

(θ, history) = simulate(T=100)
#history[history.step .== T, :] |> display
println()
θ |> display

x = 0.0:0.01:60.0
figs = Vector{Plots.Plot{Plots.GRBackend}}(undef, 1)
figs[1] = plot(xlabel="Price", ylabel="Demand pdf")
for v in eachrow(θ)
    figs[1] = plot!(figs[1], x, pdf.(Gamma(v[:α], 1/v[:β]), x),
        label="price=$(v[:price])")
end
plot(figs..., layout=(1,1))
savefig(joinpath(@__DIR__, "thomson.png"))


text:
module rectangular_plate_examples
using Statistics
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using PGFPlotsX
using Test

E = 210e3*phun("MPa");
nu = 0.3;
rho = 7850*phun("KG/M^3");
a=4.00*phun("M"); b=1.00*phun("M"); h= 0.1*phun("M");
na= 10; nb=  5; nh =2;
na= 2*10; nb=  2*5; nh =8;
neigvs = 20                   # how many eigenvalues
OmegaShift = (10.0*2*pi)^2;
# Fundamental frequency
f_analytical =[0 0 0  0 0 0 33.78 82.28  92.99 170.06];

function rectangular_plate_esnice()
    global a,b,h,na,nb,nh
    MR = DeforModelRed3D
    fens, fes = T4block(a,b,h,na,nb,nh);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    @show minimum(vec(femm.nphis)), maximum(vec(femm.nphis))
    @pgf a = Axis({
            xlabel = "Entity",
            ylabel = "Stabilization factor",
            grid="major",
            legend_pos  = "north east"
        },
        Plot({"only marks", mark="+"}, Table([:x => vec(1:count(fes)), :y => vec(femm.ephis)]))
        )
    display(a)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")
    println("f/f_analytical = $(fs[7:10]  ./ f_analytical[7:10] .* 100) %")

    vectors = []
    for i = 7:length(fs)
        scattersysvec!(u, v[:, i])
        push!(vectors, ("Mode_$i", deepcopy(u.values)))
    end
    File  =   "rectangular_plate_esnice.vtk"
    vtkexportmesh(File, connasarray(fes), fens.xyz, FinEtools.MeshExportModule.T4; vectors = vectors)
    @async run(`"paraview.exe" $File`)

    true
end # function

function allrun()
    println("#####################################################")
    println("# rectangular_plate_esnice ")
    rectangular_plate_esnice()
    return true
end # function allrun

end # module LE11NAFEMS_examples


text:
#Reference https://math.stackexchange.com/questions/1435448/negacyclic-fft-multiplication

module MulFFT
    using FFTW
    struct TwistFFTPlan
        fftplan::FFTW.Plan
        ifftplan::FFTW.Plan
        twist::Array{ComplexF64,1}
    end
    function TwistFFTPlan(N::Int32)
        TwistFFTPlan(FFTW.plan_fft(Float64.(rand(-2^9:2^9,div(N,2)))+im*rand(-2^9:2^9,div(N,2))),FFTW.plan_ifft(Float64.(rand(-2^9:2^9,div(N,2)))+im*rand(-2^9:2^9,div(N,2))),TwistGen(N))
    end
    function TwistGen(N::Int32)
        return [exp(im*k*pi/N) for k in 0:div(N,2)-1]
    end

    function TwistFFT(a::Array{Int32,1},fftplan::FFTW.Plan,twist::Array{Complex{Float64}})
        return fftplan*(complex.(Float64.(a[1:div(length(a),2)]),Float64.(a[div(length(a),2)+1:length(a)])).* twist)
    end
    function TwistIFFT(a::Array{Complex{Float64},1},ifftplan::FFTW.Plan,twist::Array{Complex{Float64}})
        b = (ifftplan*a).* conj(twist)
        return append!(real(b),imag(b))
    end

    function PolyMul(a::Array{Int32,1},b::Array{Int32,1},TwistFFTPlan::TwistFFTPlan)
        return map(x->UInt32(mod(round(Int64,x),Int64(2)^32)),TwistIFFT(TwistFFT(a,TwistFFTPlan.fftplan,TwistFFTPlan.twist).*TwistFFT(b,TwistFFTPlan.fftplan,TwistFFTPlan.twist),TwistFFTPlan.ifftplan,TwistFFTPlan.twist))
    end
end


text:
using Test
using IntervalArithmetic

@testset "widened interval is superset" begin
    a = unit_ball(Float64)
    @test a ⊂ widened(a)
    @test !(widened(a) ⊂ widened(a))
    @test a ⪽ widened(a)
    @test !(widened(a) ⪽ widened(a))
    @test a ⊆ widened(a)
    @test widened(a) ⊆ widened(a)
end

@testset "sign in unit ball" begin
    a = unit_ball(Float64)
    b = nonnegatives(Float64)
    c = interval(2)
    @test sign(c) ⊆ a
    @test sign(-c) ⊆ a
    @test sign(b) ⊆ a
    @test sign(a) ⊆ a
end


text:
using IntervalArithmetic
using Test

using LinearAlgebra

@testset "Complex interval operations" begin
    a = @interval 1im
    b = @interval 4im + 3
    c = (@interval -1 4) + (@interval 0 2)*im

    @test typeof(a) == Complex{IntervalArithmetic.Interval{Float64}}
    @test a ==  Interval(0) + Interval(1)*im
    @test a * a == Interval(-1)
    @test a + a == Interval(2)*im
    @test a - a == 0
    @test a / a == 1

    @test 3+2im ∈ c
    @test a ∪ b == (@interval 0 3) + (@interval 1 4)*im
    @test c ∩ (a ∪ b) == (@interval 0 3) + (@interval 1 2)*im
    @test a ∩ b == ∅ + ∅*im
    @test isdisjoint(a,b) == true
end

@testset "Complex functions" begin
    Z = (3 ± 1e-7) + (4 ± 1e-7)*im
    @test sin(Z) == complex(sin(real(Z)) * cosh(imag(Z)), sinh(imag(Z)) * cos(real(Z)))
    @test exp(-im * Interval(π)) == Interval(-1.0, -0.9999999999999999) - Interval(1.224646799147353e-16, 1.2246467991473532e-16)*im

    sZ = sqrt(Z)
    @test sZ == Interval(1.99999996999999951619,2.00000003000000070585) + Interval(0.99999996999999984926,1.00000003000000048381)*im
    @test sqrt(-Z) == imag(sZ) - real(sZ)*im

    @test sqrt((@interval -1 0) + (@interval 0)*im) .== (@interval 0 1)*im
    @test sqrt((@interval -1 1) + (@interval 0)*im) .== (@interval 0 1) + (@interval 0 1)*im
    @test sqrt((@interval -9//32 Inf)*im + (@interval 0)) .== (@interval 0 Inf) + (@interval -3//8 Inf)*im
end

@testset "Complex powers" begin
    x = (3..3) + 4im
    @test x^2 == -7 + 24im
    @test sqrt(x) ⊆ x^0.5
    @test x^-2 == inv(x)^2

    a = -3.1
    @test x ⊆ (x^a)^(1/a)
end

@testset "abs2 and abs" begin
    x = (0..3) + (0..4)*im
    @test abs2(x) == 0..25
    @test abs(x) == norm(x) == 0..5

    y = (-1..1) + (-2..2)*im
    @test abs(y).lo == 0.0
    @test abs2(y).lo == 0.0
end


text:
using Pkg.TOML, Libdl

config_toml = joinpath(first(DEPOT_PATH), "prefs", "MPI.toml")
mkpath(dirname(config_toml))

if !isfile(config_toml)
    touch(config_toml)
end

# This file contains code that runs at three different occasions:
#
# 1. When `build MPI` is called. At this time, `MPI.jl` has not yet
#    been loaded. This code runs essentially as a script. However, it
#    calls the MPI function `MPI_Get_library_version`, which means
#    that the MPI library must be callable. (`MPI_Init` is not
#    called.) The result of this call is saved, to ensure that the MPI
#    library doesn't accidentally change at a later time (e.g. when
#    the wrong set of modules is loaded), which could lead to
#    segfaults that are difficult to debug.
#
# 2. When `MPI.jl` is precompiled. Most things should be initialized
#    at that time.
#
# 3. When `MPI.jl` is loaded, i.e. when its `__init__` function is
#    run. Some information is only available at load time, e.g.
#    certatin MPI constants such as `MPI_COMM_DUP_FN`. This might be
#    function pointers whose value is only known at load time. Another
#    reason to run code at this time would be to set environment
#    variables; their values are not captured while precompiling.
#
# Generally, this file `build.jl` is run at stage (1). It creates a
# file `deps.jl` that is run at stage (2), and which also defines a
# function `__init__deps` that is run at stage (3). Certain actions
# need to happen at two or all three stages.

config = TOML.parsefile(config_toml)
update_config = false

# MPI.toml has 4 keys
#  binary   = "" (default) | "system" | "MPICH_jll" | "MPItrampoline_jll" | "OpenMPI_jll" | "MicrosoftMPI_jll"
#  path     = "" (default) | top-level directory location
#  library  = "" (default) | library name/path | list of library names/paths
#  mpiexec  = "" (default) | executable name/path | [executable name/path, extra args...]

# 1. check if environment variables have changed
if haskey(ENV, "JULIA_MPI_BINARY")
    config["binary"] = ENV["JULIA_MPI_BINARY"]
    update_config = true
elseif haskey(ENV, "JULIA_MPI_LIBRARY") || haskey(ENV, "JULIA_MPI_PATH")
    config["binary"] = "system"
end

if haskey(ENV, "JULIA_MPI_PATH")
    config["path"] = ENV["JULIA_MPI_PATH"]
    update_config = true
end

if haskey(ENV, "JULIA_MPI_LIBRARY")
    config["library"] = ENV["JULIA_MPI_LIBRARY"]
    update_config = true
end

if haskey(ENV, "JULIA_MPIEXEC")
    if ENV["JULIA_MPIEXEC"] == ""
        delete!(config, "mpiexec")
    else
        config["mpiexec"] = ENV["JULIA_MPIEXEC"]
    end
    update_config = true
end

if update_config
    open(config_toml, write=true) do f
        TOML.print(f, config)
    end
end

binary = get(config, "binary", "")

# 2. generate deps.jl
if binary == "system"

    @info "using system MPI"
    library = get(config, "library", "")
    if library == ""
        library = ["libmpi", "libmpi_ibm", "msmpi", "libmpich", "libmpitrampoline"]
    end
    path    = get(config, "path", "")
    mpiexec = get(config, "mpiexec", path == "" ? "mpiexec" : joinpath(path, "bin", "mpiexec"))

    const libmpi = find_library(library, path == "" ? [] : [joinpath(path, "lib"), joinpath(path, "lib64")])
    if libmpi == ""
        error("libmpi could not be found")
    end

    const mpiexec_cmd = Cmd(mpiexec isa String ? [mpiexec] : mpiexec)

    _doc_external(fname) = ""

    include(joinpath("..", "src", "implementations.jl"))

    @info "using implementation" libmpi mpiexec_cmd MPI_LIBRARY_VERSION_STRING

    include("prepare_mpi_constants.jl")

    deps = quote
        const libmpi = $libmpi
        const libmpiconstants = joinpath(dirname(@__DIR__), "deps", "libload_time_mpi_constants.so")
        const mpiexec_cmd = $mpiexec_cmd
        const mpiexec_path = mpiexec_cmd[1]
        _mpiexec(fn) = fn(mpiexec_cmd)

        using Requires

        function __init__deps()
            if Get_library_version() != $MPI_LIBRARY_VERSION_STRING
                error("MPI library has changed, please re-run `Pkg.build(\"MPI\")`.")
            end

            # error if a jll is loaded.
            @require(MPICH_jll         = "7cb0a576-ebde-5e09-9194-50597f1243b4",
                     error("MPICH_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
            @require(MPItrampoline_jll = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748",
                     error("MPItrampoline_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
            @require(MicrosoftMPI_jll  = "9237b28f-5490-5468-be7b-bb81f5f5e6cf",
                     error("MicrosoftMPI_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
            @require(OpenMPI_jll       = "fe0851c0-eecd-5654-98d4-656369965a5c",
                     error("OpenMPI_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
        end
    end

elseif binary == ""

    @info "using default MPI jll"

    if Sys.iswindows()
        using MicrosoftMPI_jll
        # run(MicrosoftMPI_jll.generate_compile_time_mpi_constants())
        cp("compile_time_mpi_constants_msmpi_$(Sys.ARCH).jl", "compile_time_mpi_constants.jl"; force=true)
    else
        using MPICH_jll
        run(MPICH_jll.generate_compile_time_mpi_constants())
    end

    deps = quote
        if Sys.iswindows()
            using MicrosoftMPI_jll
            const libmpiconstants = MicrosoftMPI_jll.libload_time_mpi_constants
            const _mpiexec = MicrosoftMPI_jll.mpiexec
            const mpiexec_path = MicrosoftMPI_jll.mpiexec_path
        else
            using MPICH_jll
            const libmpiconstants = MPICH_jll.libload_time_mpi_constants
            const _mpiexec = MPICH_jll.mpiexec
            const mpiexec_path = MPICH_jll.mpiexec_path
        end

        function __init__deps()
            if (haskey(ENV, "SLURM_JOBID") || haskey(ENV, "PBS_JOBID") || haskey(ENV, "LSB_JOBID")) &&
                get(ENV, "JULIA_MPI_CLUSTER_WARN", "") != "n"
                @warn """
                    You appear to be using MPI.jl with the default MPI binary on a cluster.
                    We recommend using the system-provided MPI, see the Configuration section of the MPI.jl docs.
                    """
            end
        end
    end

elseif binary == "MPICH_jll"

    @info "using MPICH_jll"

    using MPICH_jll
    run(MPICH_jll.generate_compile_time_mpi_constants())

    deps = quote
        using MPICH_jll
        const libmpiconstants = MPICH_jll.libload_time_mpi_constants
        const _mpiexec = MPICH_jll.mpiexec
        const mpiexec_path = MPICH_jll.mpiexec_path

        __init__deps() = nothing
    end

elseif binary == "MPItrampoline_jll"

    @info "using MPItrampoline_jll"
    using MPItrampoline_jll
    @assert MPItrampoline_jll.is_available()
    run(MPItrampoline_jll.generate_compile_time_mpi_constants())

    deps = quote
        @info "using MPItrampoline_jll"
        using MPItrampoline_jll
        @assert MPItrampoline_jll.is_available()
        mpiexec_path = get(ENV, "MPITRAMPOLINE_MPIEXEC", nothing)
        if mpiexec_path === nothing
            if "MPITRAMPOLINE_LIB" ∉ keys(ENV)
                @info "[MPI] Using built-in MPICH with MPItrampoline"
                # In principle, MPItrampoline_jll.mpiwrapperexec should
                # forward to MPItrampoline_jll.mpich_mpiexec. However,
                # there is no mechanism to update the path to where the
                # Julia artifact is installed, so we forward manually
                # here.
                const _mpiexec = MPItrampoline_jll.mpich_mpiexec
                mpiexec_path = MPItrampoline_jll.mpich_mpiexec_path
            else
                const _mpiexec = MPItrampoline_jll.mpiwrapperexec
                mpiexec_path = MPItrampoline_jll.mpiwrapperexec_path
            end
        else
            const mpiexec_cmd = Cmd([mpiexec_path])
            _mpiexec(fn) = fn(mpiexec_cmd)
        end
        const libmpiconstants = MPItrampoline_jll.libload_time_mpi_constants_path

        __init__deps() = nothing
    end

elseif binary == "OpenMPI_jll"

    @info "using OpenMPI_jll"

    using OpenMPI_jll
    run(OpenMPI_jll.generate_compile_time_mpi_constants())

    deps = quote
        using OpenMPI_jll
        const _mpiexec = OpenMPI_jll.mpiexec
        const mpiexec_path = OpenMPI_jll.mpiexec_path

        __init__deps() = nothing
    end

elseif binary == "MicrosoftMPI_jll"

    @info "using MicrosoftMPI_jll"

    using MicrosoftMPI_jll
    # run(MicrosoftMPI_jll.generate_compile_time_mpi_constants())
    cp("compile_time_mpi_constants_msmpi_$(Sys.ARCH).jl", "compile_time_mpi_constants.jl"; force=true)

    deps = quote
        using MicrosoftMPI_jll
        const libmpiconstants = MicrosoftMPI_jll.libload_time_mpi_constants
        const _mpiexec = MicrosoftMPI_jll.mpiexec
        const mpiexec_path = MicrosoftMPI_jll.mpiexec_path

        __init__deps() = nothing
    end

else
    error("Unknown binary $binary")
end

remove_line_numbers(x) = x
function remove_line_numbers(ex::Expr)
    if ex.head == :macrocall
        ex.args[2] = nothing
    else
        ex.args = [remove_line_numbers(arg) for arg in ex.args if !(arg isa LineNumberNode)]
    end
    return ex
end

# Only update deps.jl if it has changed.
# This allows users to call Pkg.build("MPI") without triggering another round of precompilation.
deps_str =
    """
    # This file has been generated automatically.
    # It will be overwritten the next time `Pkg.build("MPI")` is called.
    """ *
    string(remove_line_numbers(deps)) *
    """
    """

if !isfile("deps.jl") || deps_str != read("deps.jl", String)
    write("deps.jl", deps_str)
end