text:
module Axes
export CubeAxis, QuantileAxis, TimeAxis, TimeHAxis, VariableAxis, LonAxis, LatAxis, CountryAxis,
SpatialPointAxis,Axes,YearStepRange,CategoricalAxis,RangeAxis,axVal2Index,MSCAxis,
ScaleAxis, axname, @caxis_str, findAxis, AxisDescriptor, get_descriptor, ByName, ByType, ByValue, ByFunction, getAxis,
getOutAxis, ByInference, renameaxis!, axsym
import NetCDF.NcDim
using ..Cubes
import ..Cubes: caxes, _read
using Dates

macro defineCatAxis(axname,eltype)
  newname=esc(Symbol(string(axname,"Axis")))
  quote
    const $newname = CategoricalAxis{T,$(QuoteNode(axname)),R} where R<:AbstractVector{T} where T<:$(eltype)
    $(newname)(v::AbstractVector{T}) where T = CategoricalAxis{T,$(QuoteNode(axname)),typeof(v)}(v)
  end
end

macro defineRanAxis(axname,eltype)
  newname=esc(Symbol(string(axname,"Axis")))
  quote
    const $newname = RangeAxis{T,$(QuoteNode(axname)),R} where R<:AbstractVector{T} where T<:$(eltype)
    $(newname)(v::AbstractVector{T}) where T = RangeAxis{T,$(QuoteNode(axname)),typeof(v)}(v)
  end
end

"""
    abstract CubeAxis{T} <: AbstractCubeData{T,1}

Supertype of all axes. Every `CubeAxis` is an 1D Cube itself and can be passed
to mapCube operations. Although all cube axes are instances of the parametric typealias
`CategoricalAxis` and `RangeAxis`, there are some typealiases defined
to provide shorter and more convenient names for commonly used cube axes. Here is a list
of the aliases:

### Categorical Axes

* `VariableAxis` represents different variables
* `SpatialPointAxis` represents a list of coordinates
* `CountryAxis` countries
* `ScaleAxis` time scales after time series decomposition
* `QuantileAxis` represents different quantiles

### Continuous Axes

* `LonAxis` longitudes
* `LatAxis` latitudes
* `TimeAxis` time
* `MSCAxis` time step inside a year (for seasonal statistics)

"""
abstract type CubeAxis{T,S} <: AbstractCubeMem{T,1} end

Base.size(x::CubeAxis)=(length(x.values),)
Base.size(x::CubeAxis,i)=i==1 ? length(x.values) : error("Axis has only a single dimension")
Base.ndims(x::CubeAxis)=1

"""
    CategoricalAxis{T,S}

To represent axes that are categorical, where `T` is the element type.
The type parameter `S` denotes the axis name (a symbol).
The default constructor is:

    CategoricalAxis(axname::String,values::Vector{T})

"""
struct CategoricalAxis{T,S,RT} <: CubeAxis{T,S}
  values::RT
end

CategoricalAxis(s::Symbol,v)=CategoricalAxis{eltype(v),s,typeof(v)}(v)
CategoricalAxis(s::AbstractString,v)=CategoricalAxis(Symbol(s),v)

@defineCatAxis Variable String
@defineCatAxis SpatialPoint Tuple{Number,Number}
@defineCatAxis Scale String
@defineCatAxis Quantile AbstractFloat

struct RangeAxis{T,S,R<:AbstractVector{T}} <: CubeAxis{T,S}
  values::R
end

"""
    RangeAxis{T,S,R}

To represent axes that are categorical, where `T` is the element type.
The type parameter `S` denotes the axis name (a symbol) and `R` the type of the
range which is used to represent the axis values.
The default constructor is:

    RangeAxis(axname::String,values::Range{T})

"""
RangeAxis(s::Symbol,v::AbstractVector{T}) where T = RangeAxis{T,s,typeof(v)}(v)
RangeAxis(s::AbstractString,v)=RangeAxis(Symbol(s),v)

@defineRanAxis MSC TimeType
@defineRanAxis Lon Number
@defineRanAxis Lat Number
@defineRanAxis Time TimeType

Base.length(a::CubeAxis)=length(a.values)

MSCAxis(n::Int)=MSCAxis(DateTime(1900):Day(ceil(Int,366/n)):DateTime(1900,12,31,23,59,59))

caxes(x::CubeAxis)=CubeAxis[x]

axname(::CategoricalAxis{T,S}) where {T,S}=string(S)
axname(::RangeAxis{T,S}) where {T,S}=string(S)
axname(::Type{T}) where T<:CubeAxis{S,U} where {S,U} = U
axsym(ax::CubeAxis{<:Any,S}) where S = S
axunits(::CubeAxis)="unknown"
axunits(::LonAxis)="degrees_east"
axunits(::LatAxis)="degrees_north"

get_step(r::AbstractRange)=step(r)
get_step(r::AbstractVector)=length(r)==0 ? zero(eltype(r)) : r[2]-r[1]

axVal2Index_ub(a::RangeAxis, v; fuzzy=false)=axVal2Index(a,v-abshalf(get_step(a.values)),fuzzy=fuzzy)
axVal2Index_lb(a::RangeAxis, v; fuzzy=false)=axVal2Index(a,v+abshalf(get_step(a.values)),fuzzy=fuzzy)

axVal2Index_ub(a::RangeAxis, v::Date; fuzzy=false)=axVal2Index(a,DateTime(v)-abshalf(get_step(a.values)),fuzzy=fuzzy)
axVal2Index_lb(a::RangeAxis, v::Date; fuzzy=false)=axVal2Index(a,DateTime(v)+abshalf(get_step(a.values)),fuzzy=fuzzy)

abshalf(a) = abs(a/2)
abshalf(a::Day) = abs(Millisecond(a)/2)
abshalf(a::Month) = iseven(Dates.value(a)) ? a/2 : Month(a÷2) + Day(15)

get_bb(ax::RangeAxis) = first(ax.values)-abshalf(get_step(ax.values)), last(ax.values)+abshalf(get_step(ax.values))
function axisfrombb(name,bb,n)
  offs = (bb[2]-bb[1])/(2*n)
  RangeAxis(name,range(bb[1]+offs,bb[2]-offs,length=n))
end

function axVal2Index(a::RangeAxis{<:Any,<:Any,<:AbstractRange},v;fuzzy=false)
  dt = v-first(a.values)
  r = round(Int,dt/step(a.values))+1
  return max(1,min(length(a.values),r))
end

convert_time(T::Type{<:TimeType}, v::TimeType) = T(year(v), month(v), day(v), hour(v), minute(v), second(v))
convert_time(T::Type{<:TimeType}, v::Date) = T(year(v), month(v), day(v), 0, 0, 0)
convert_time(::Type{Date},v::TimeType) = Date(year(v),month(v),day(v))
convert_time(::Type{Date},v::Date) = Date(year(v),month(v),day(v))

function axVal2Index(a::RangeAxis{T},v;fuzzy=false) where T<:TimeType
  vconverted = convert_time(T,v)
  dd = map(i->abs((i-vconverted)),a.values)
  mi,ind = findmin(dd)
  return ind
end
function axVal2Index(a::RangeAxis{T,<:Any,<:AbstractRange},v;fuzzy=false) where T<:TimeType
  vconverted = convert_time(T,v)
  dd = map(i->abs((i-vconverted)),a.values)
  mi,ind = findmin(dd)
  return ind
end
function axVal2Index(axis::CategoricalAxis{String},v::String;fuzzy::Bool=false)
  r=findfirst(isequal(v),axis.values)
  if r==nothing
    if fuzzy
      r=findall(axis.values) do i
        startswith(lowercase(i),lowercase(v[1:min(length(i),length(v))]))
      end
      if length(r)==1
        return(r[1])
      else
        error("Could not find unique value of $v in $axis")
      end
    else
      error("$v not found in $axis")
    end
  end
  r
end
axVal2Index(x,v::CartesianIndex{1};fuzzy::Bool=false)=min(max(v.I[1],1),length(x))
function axVal2Index(x,v;fuzzy::Bool=false)
  i = findfirst(isequal(v),x.values)
  isa(i,Nothing) && error("Value $v not found in x")
  return i
end
abstract type AxisDescriptor end
getAxis(d::Any,v::Any)=error("getAxis not defined for $d $v")
struct ByName <: AxisDescriptor
  name::String
end
struct ByInference <: AxisDescriptor end
struct ByType{T} <: AxisDescriptor
  t::Type{T}
end
struct ByValue <: AxisDescriptor
  v::CubeAxis
end
struct ByFunction <: AxisDescriptor
  f::Function
end

findAxis(a,c::AbstractCubeData)=findAxis(a,caxes(c))
get_descriptor(a::String)=ByName(a)
get_descriptor(a::Type{T}) where {T<:CubeAxis}=ByType(a)
get_descriptor(a::CubeAxis)=ByValue(a)
get_descriptor(a::Function)=ByFunction(a)
get_descriptor(a)=error("$a is not a valid axis description")
get_descriptor(a::AxisDescriptor)=a

const VecOrTuple{S} = Union{Vector{<:S},NTuple{<:Any,S}} where S

"Find a certain axis type in a vector of Cube axes and returns the index"
function findAxis(bt::ByType,v::VecOrTuple{S}) where S<:CubeAxis
  a=bt.t
  for i=1:length(v)
    isa(v[i],a) && return i
  end
  return nothing
end
function findAxis(bs::ByName,axlist::VecOrTuple{T}) where T<:CubeAxis
  matchstr=bs.name
  ism=map(i->startswith(lowercase(axname(i)),lowercase(matchstr)),axlist)
  sism=sum(ism)
  sism==0 && return nothing
  sism>1 && error("Multiple axes found matching string $matchstr")
  i=findfirst(ism)
end
function findAxis(bv::ByValue,axlist::VecOrTuple{T}) where T<:CubeAxis
  v=bv.v
  return findfirst(i->i==v,axlist)
end
function getAxis(desc,axlist::Vector{T}) where T<:CubeAxis
  i = findAxis(desc,axlist)
  if isa(i,Nothing)
    return nothing
  else
    return axlist[i]
  end
end
getOutAxis(desc,axlist,incubes,pargs,f) = getAxis(desc,unique(axlist))
function getOutAxis(desc::ByFunction,axlist,incubes,pargs,f)
  outax = desc.f(incubes,pargs)
  isa(outax,CubeAxis) || error("Axis Generation function $(desc.f) did not return an axis")
  outax
end
import DataStructures: counter
function getOutAxis(desc::Tuple{ByInference},axlist,incubes,pargs,f)
  inAxes = map(caxes,incubes)
  inAxSmall = map(i->filter(j->in(j,axlist),i) |>collect,inAxes)
  inSizes = map(i->(map(length,i)...,),inAxSmall)
  intypes = map(eltype, incubes)
  testars = map((s,it)->zeros(it,s...),inSizes, intypes)
  map(testars) do ta
    ta .= rand(Base.nonmissingtype(eltype(ta)),size(ta)...)
    if eltype(ta) >: Missing
      # Add some missings
      randind = rand(1:length(ta),length(ta)÷10)
      ta[randind] .= missing
    end
  end
  resu = f(testars...,pargs...)
  isa(resu,AbstractArray) || isa(resu,Number) || error("Function must return an array or a number")
  isa(resu,Number) && return ()
  outsizes = size(resu)
  outaxes = map(outsizes,1:length(outsizes)) do s,il
    if s>2
      i = findall(i->i==s,length.(axlist))
      if length(i)==1
        return axlist[i[1]]
      elseif length(i)>1
        @info "Found multiple matching axes for output dimension $il"
      end
    end
    return RangeAxis("OutAxis$(il)",1:s)
  end
  if !allunique(outaxes)
    #TODO: fallback with axis renaming in this case
    error("Could not determine unique output axes from output shape")
  end
  return (outaxes...,)
end
"""
    getAxis(desc::String, c::AbstractCubeData)

Given the string of an axis name and a cube, returns this axis of the cube.
"""
getAxis(desc,c::AbstractCubeData)=getAxis(desc,caxes(c))
getAxis(desc::ByValue,axlist::Vector{T}) where {T<:CubeAxis}=desc.v

"Fallback method"
findAxis(a,axlist)=findAxis(get_descriptor(a),axlist)

getSubRange(x::CubeAxis,i)=x.values[i],nothing
getSubRange(x::TimeAxis,i)=view(x,i),nothing

renameaxis(r::RangeAxis{T,<:Any,V}, newname) where {T,V} = RangeAxis{T,Symbol(newname),V}(r.values)
renameaxis(r::CategoricalAxis{T,<:Any,V}, newname) where {T,V} = CategoricalAxis{T,Symbol(newname),V}(r.values)
function renameaxis!(c::AbstractCubeData,p::Pair)
  i = findAxis(p[1],c.axes)
  c.axes[i]=renameaxis(c.axes[i],p[2])
  c
end
function renameaxis!(c::AbstractCubeData,p::Pair{<:Any,<:CubeAxis})
  i = findAxis(p[1],c.axes)
  i === nothing && throw(ArgumentError("Axis not found"))
  length(c.axes[i].values) == length(p[2].values) || throw(ArgumentError("Length of replacement axis must equal length of old axis"))
  c.axes[i]=p[2]
  c
end

function _read(ax::CubeAxis, ar::AbstractArray, I::CartesianIndices)
  ar[:] .= ax.values[I.indices[1]]
end

macro caxis_str(s)
  :(CategoricalAxis{String,$(QuoteNode(Symbol(s)))})
end

import Base.==
import Base.isequal
==(a::CubeAxis,b::CubeAxis)=(a.values==b.values) && (axname(a)==axname(b))
isequal(a::CubeAxis, b::CubeAxis) = a==b

import CFTime: timeencode
function NcDim(a::CubeAxis{T},start::Integer,count::Integer) where T<:TimeType
  if start + count - 1 > length(a.values)
    count = oftype(count,length(a.values) - start + 1)
  end
  tv=a.values[start:(start+count-1)]
  startyear=Dates.year(first(a.values))
  starttime=T(startyear)
  timeunits = "days since $startyear-01-01"
  atts=Dict{Any,Any}("units"=>timeunits, "calendar"=>"standard")
  d = timeencode(tv, timeunits)
  NcDim(axname(a),length(d),values=d,atts=atts)
end

#Default constructor
NcDim(a::CubeAxis{T},start::Integer,count::Integer) where {T<:Real}=NcDim(axname(a),count,values=collect(a.values[start:(start+count-1)]),atts=Dict{Any,Any}("units"=>axunits(a)))
NcDim(a::CubeAxis,start::Integer,count::Integer)=NcDim(axname(a),count,values=Float64[start:(start+count-1);],atts=Dict{Any,Any}("units"=>axunits(a)))
NcDim(a::CubeAxis)=NcDim(a,1,length(a))

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize.jl
# Load subproblem type optimization routines and define general optimize!
# throw_optimize_hook!.
#############################################################################

include("optimize_lp.jl")
include("optimize_conic.jl")
include("optimize_convex.jl")
include("optimize_nonconvex.jl")

function throw_optimize_hook!(m::Optimizer{Q,S,T}) where {Q,S,T}
    optimize_hook!(_ext(m), m)
end

function MOI.optimize!(m::Optimizer{Q,S,T}) where {Q,S,T}

    m._global_optimizer._start_time = time()
    # Runs the branch and bound routine
    if !m.enable_optimize_hook

        # Need to eliminate fixed variables prior to checking
        # for convex quadratic constraints
        initial_parse!(m)

        # Determines if the problem is an LP, MILP, SOCP, MISCOP,
        # CONVEX, OF MINCVX PROBLEM TYPE
        parse_classify_problem!(m._global_optimizer)

        m._global_optimizer._parse_time = m._global_optimizer._start_time - time()

        # Throws the problem to the appropriate solution routine
        optimize!(m._working_problem._problem_type, m)
    else

        # throws to user-defined optimization hook
        throw_optimize_hook!(m)
    end

    return
end

ext: jl

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using Rhapsodie
using DelimitedFiles
using EasyFITS

#include("test_separable_reconstruction.jl")

tau=0.03;

par=readdlm("test_results/Parameters.txt")
DSIZE=Int64(par[1]);
NTOT=Int64(par[2]);
Nframe=Int64(par[3]);
Nrot=Int64(par[4]);
Nangle=NTOT÷(Nframe*4)
Center=par[5:6];
DerotAng=Vector{Float64}();
Epsilon=Vector{Tuple{Array{Float64,1},Array{Float64,1}}}();

for iter=1:NTOT
    ind=div(iter-1, NTOT/4)+1
    push!(Epsilon,([0. ,0. ],par[end-1:end]));
    push!(DerotAng, deg2rad(-par[7]));
end

psf_center=readdlm("../data/PSF_centers_Airy.txt");

load_parameters((DSIZE, 2*DSIZE, NTOT), Nframe, Nrot, Nangle, Center, (psf_center[1:2], psf_center[3:4]), Epsilon, DerotAng)

Load_Data("test_results/DATA_$tau-$DSIZE.fits",
                    "test_results/WEIGHT_$tau-$DSIZE.fits")

PSF=read(FitsArray, "../data/PSF_parametered_Airy.fits");
const A=set_fft_op(PSF[1:end÷2,:]'[:,:],psf_center[1:2]);

X0 = PolarimetricMap("mixed", zeros(Rhapsodie.get_par().cols));
#X0 = read("mixed", "test_results/Results_Separable_Linear_$tau-$DSIZE.fits");

regularisation_parameters = [0.5 , -1. , -1., -3.]; #(in log10)

@time x = apply_rhapsodie(X0, A, Rhapsodie.dataset, regularisation_parameters,
                          maxeval=1000, maxiter=1000);
crop!(x)
write(x, "test_results/RHAPSODIE_nonlinearresults_$tau-$DSIZE.fits")

max_stars_repo_name: LaurenceDenneulin/Rhapsodie

max_issues_repo_name: LaurenceDenneulin/Rhapsodie

max_forks_repo_name: LaurenceDenneulin/Rhapsodie


text:
module LazyExpressionTest

using Test
using Random
using LinearAlgebra
using Parametron
using StaticArrays

import Random
import Parametron: setdirty!, mock_model

@testset "basics" begin
    model = mock_model()
    a = 2
    b = Parameter{Int}(() -> 3, model)
    c = 4
    expr = @expression(a + b * c)
    @test expr() == a + b() * c
end

@testset "bad expression" begin
    @test_throws(ArgumentError, @expression x ? y : z)
end

@testset "parameter" begin
    model = mock_model()
    a = 3
    b = 4.0
    cval = Ref(0)
    c = Parameter{Int}(() -> cval[], model)
    expr = @expression(a + b * c)

    @test expr() == a + b * cval[]
    cval[] = 4
    setdirty!(c)
    @test expr() == a + b * cval[]
end

@testset "nested" begin
    model = mock_model()
    a = 3
    b = 4.0
    cval = Ref(5)
    c = Parameter{Int}(() -> cval[], model)
    expr1 = @expression(a + b * c)
    expr2 = @expression(4 * expr1)
    @test expr2() == 4 * expr1()
    show(devnull, expr1)
end

module M
export SpatialMat, angular, linear
struct SpatialMat
    angular::Matrix{Float64}
    linear::Matrix{Float64}
end
angular(mat::SpatialMat) = mat.angular
linear(mat::SpatialMat) = mat.linear
end

using .M

@testset "user functions" begin
    mat = SpatialMat(rand(3, 4), rand(3, 4))
    scalar = Ref(1.0)
    updatemat! = let scalar = scalar # https://github.com/JuliaLang/julia/issues/15276
        mat -> (mat.angular .= scalar[]; mat.linear .= scalar[]; mat)
    end
    model = mock_model()
    pmat = Parameter(updatemat!, mat, model)
    pmat_angular = @expression angular(pmat)
    result = pmat_angular()
    @test result === angular(mat)
    @test all(result .== scalar[])

    setdirty!(model)
    allocs = @allocated begin
        setdirty!(model)
        pmat_angular()
    end
    @test allocs == 0
end

@testset "matvecmul!" begin
    m = mock_model()
    A = Parameter(rand!, zeros(3, 4), m)
    x = Variable.(1 : 4)
    expr = @expression A * x
    @test expr() == A() * x
    setdirty!(m)
    allocs = @allocated begin
        setdirty!(m)
        expr()
    end
    @test allocs == 0

    wrapped = Parametron.wrap(expr)
    @test wrapped isa Parametron.WrappedExpression{Vector{AffineFunction{Float64}}}
    setdirty!(m)
    @test wrapped() == expr()
    allocs = @allocated begin
        setdirty!(m)
        wrapped()
    end
    @test allocs == 0
    wrapped² = Parametron.wrap(wrapped)
    @test wrapped² === wrapped
end

@testset "StaticArrays" begin
    m = mock_model()
    A = Parameter{SMatrix{3, 3, Int, 9}}(m) do
        @SMatrix ones(Int, 3, 3)
    end
    x = Variable.(1 : 3)

    expr1 = @expression A * x
    @test expr1() == A() * x
    setdirty!(A)
    allocs = @allocated expr1()
    @test allocs == 0
    @test expr1() isa SVector{3, AffineFunction{Int}}

    y = Parameter(m, val=SVector{3}(x))
    expr2 = @expression y + y
    @test expr2() == y() + y()
    allocs = @allocated expr2()
    @test allocs == 0
    @test expr2() isa SVector{3, AffineFunction{Int}}

    expr3 = @expression y - y
    @test expr3() == y() - y()
    allocs = @allocated expr3()
    @test allocs == 0
    @test expr3() isa SVector{3, AffineFunction{Int}}
end

@testset "mul! optimization" begin
    m = mock_model()
    weight = Parameter(() -> 3, m)
    x = Variable.(1 : 3)
    expr = @expression weight * (x ⋅ x)
    vals = Dict(zip(x, [1, 2, 3]))
    xvals = map(var -> vals[var], x)
    @test expr()(vals) == 3 * xvals ⋅ xvals
    allocs = @allocated expr()
    @test allocs == 0
end

@testset "scale! optimization" begin
    model = mock_model()
    x = [Variable(model) for _ in 1 : 3]
    dt = Parameter{Float64}(model) do
        2.0
    end
    @testset "Number" begin
        expr = @expression dt * [1, 2, 3, 4, 5]
        @test expr() == dt() * [1, 2, 3, 4, 5]
        @test @allocated(expr()) == 0

        expr = @expression [1, 2, 3, 4, 5] * dt
        @test expr() == [1, 2, 3, 4, 5] * dt()
        @test @allocated(expr()) == 0
    end
    @testset "Variable" begin
        expr = @expression dt * x
        @test expr() == dt() * x
        @test @allocated(expr()) == 0

        expr = @expression x * dt
        @test expr() == x * dt()
        @test @allocated(expr()) == 0
    end

    @testset "AffineFunction" begin
        Ax = ones(3, 3) * x
        expr = @expression dt * Ax
        @test expr() == dt() * Ax
        @test @allocated(expr()) == 0

        expr = @expression Ax * dt
        @test expr() == Ax * dt()
        @test @allocated(expr()) == 0
    end
end

@testset "vcat optimization" begin
    Random.seed!(42)
    m = mock_model()
    A = Parameter(rand!, zeros(3, 4), m)
    B = Parameter(rand!, zeros(3, 3), m)
    x = Variable.(1 : 4)
    y = Variable.(5 : 7)
    f1 = @expression A * x
    f2 = @expression B * y

    v1 = @expression vcat(f1, f2)
    @test v1() == vcat(f1(), f2())
    setdirty!(m)
    @test (@allocated begin
        setdirty!(m)
        v1()
    end) == 0

    # Make sure we expand vcat expressions
    v2 = @expression [f1; f2]
    @test v2() == vcat(f1(), f2())
    @test (@allocated begin
        setdirty!(m)
        v2()
    end) == 0

    # Make sure static arrays still work
    C = Parameter{SMatrix{2, 2, Int, 4}}(m) do
        @SMatrix [1 2; 3 4]
    end
    z = SVector(Variable(8), Variable(9))
    f3 = @expression C * z
    v3 = @expression [f3; f3]
    @test v3() == vcat(f3(), f3())
    @test v3() isa SVector{4, AffineFunction{Int}}
    @test (@allocated begin
        setdirty!(m)
        v3()
    end) == 0

    # Other numbers of arguments
    v4 = @expression vcat(f3)
    @test v4() == f3()
    @test (@allocated begin
        setdirty!(m)
        v4()
    end) == 0

    v5 = @expression vcat(f1, f2, f3)
    @test v5() == vcat(f1(), f2(), f3())
    @test (@allocated begin
        setdirty!(m)
        v5()
    end) == 0

    # Generic fallbacks that allocate memory but should still give the right answer
    @test (@expression vcat(f1, 3))() == vcat(f1(), 3)
end

@testset "vect optimization" begin
    m = mock_model()
    x = Variable.(1:2)
    p = Parameter(m) do
        SVector(1., 2.)
    end
    expr = @expression [dot(p, x)]
    @test expr() == [dot(p(), x)]
    @test @allocated(expr()) == 0

    p2 = Parameter(m) do
        2.0
    end
    ex2 = @expression [p2]
    @test ex2() == [2.0]
    @test @allocated(ex2()) == 0
end

@testset "convert optimization" begin
    m = mock_model()
    A = Parameter{SMatrix{3, 3, Int, 9}}(m) do
        @SMatrix ones(Int, 3, 3)
    end
    x = Variable.(1 : 3)
    expr = @expression convert(Vector, A * x)

    @test expr() == A() * x
    @test expr() isa Vector{AffineFunction{Int}}
    allocs = @allocated expr()
    @test allocs == 0
end

@testset "vecdot optimization" begin
    m = mock_model()
    x = Variable.(1:3)
    p = Parameter(m) do
        SVector(1., 2., 3.)
    end
    ex1 = @expression dot(x, p)
    ex2 = @expression dot(p, x)
    @test ex1() == ex2() == dot(p(), x)
    @test ex1() isa AffineFunction
    @test @allocated(ex1()) == 0
    @test @allocated(ex2()) == 0

    ex3 = @expression dot(x + p, p)
    ex4 = @expression dot(p, x + p)
    @test ex3() == ex4() == dot(p(), x + p())
    @test ex3() isa AffineFunction
    @test @allocated(ex3()) == 0
    @test @allocated(ex4()) == 0
end

@testset "adjoint optimization" begin
    m = mock_model()
    x = Variable.(1:2)
    p = Parameter(m) do
        @SMatrix [1. 2.; 3. 4.]
    end
    ex = @expression adjoint(p) * x
    @test ex() == adjoint(p()) * x
    @test @allocated(ex()) == 0

    p2 = Parameter(rand!, zeros(2, 2), m)
    ex2 = @expression adjoint(p2) * x
    @test ex2() == adjoint(p2()) * x
    @test @allocated(ex2()) == 0
end

@testset "issue 26" begin
    model = mock_model()
    n = 2
    x = Variable.(1 : n)
    vals = Dict(zip(x, rand(n)))
    q = Parameter(zeros(n), model) do q
        q[1] = 1
        q[2] = 2
    end
    expr1 = @expression transpose(x) * Matrix(I, n, n) * x + q' * x
    expr2 = @expression q' * x + transpose(x) * Matrix(I, 2, 2) * x
    @test expr1()(vals) ≈ expr2()(vals) atol = 1e-14
end

@testset "issue 32" begin
    model = mock_model()
    v = [Variable(model) for _ in 1:2]
    v0 = zeros(2)
    Δt = 0.01
    u = [Variable(model) for _ in 1:2]
    H = Parameter(zeros(2, 2), model) do H
        H[1, 1] = rand()
        H[2, 2] = rand()
    end
    c = Parameter(zeros(2), model) do c
        c[1] = rand()
        c[2] = rand()
    end
    expr = @expression(H * (v - v0) - Δt * (u - c))
    @test expr() == H() * (v - v0) - Δt * (u - c())
    @test @allocated(expr()) == 0
end

@testset "issue 23, fully qualified names" begin
    model = mock_model()
    p = Parameter{Int}(() -> 2, model)
    expr = @expression Base.:+(2, p)
    @test expr() == 4
end

struct MyWrapper{T}
    x::T
end

function allocations_in_local_scope(x)
    @allocated x()
end

@testset "getfield optimization" begin
    model = mock_model()
    p = Parameter(model) do
        MyWrapper(rand())
    end
    ex1 = @expression p.x
    @test ex1() == p().x
    setdirty!(p)
    @test ex1() == p().x
    setdirty!(p)
    @test allocations_in_local_scope(ex1) == 0

    ex2 = @expression p.x + 1
    @test ex2() == p().x + 1
    setdirty!(p)
    @test ex2() == p().x + 1
    setdirty!(p)
    @test allocations_in_local_scope(ex2) == 0
end

@testset "hcat" begin
    model = mock_model()
    p = Parameter(model) do
        rand(3, 3)
    end
    hcat_expr = @expression [p p]
    @test hcat_expr() == [p() p()]
    setdirty!(p)
    @test hcat_expr() == [p() p()]
end

@testset "getindex" begin
    model = mock_model()
    p = Parameter(model) do
        rand(3, 3)
    end
    getindex_expr = @expression p[:,2]
    @test getindex_expr() == p()[:,2]
    setdirty!(p)
    @test getindex_expr() == p()[:,2]
end

@testset "A' * x (https://github.com/tkoolen/Parametron.jl/pull/79#issuecomment-423732633)" begin
    model = mock_model()
    n, m = 5, 15
    Adata = randn(n, m)
    A = Parameter(model, val=Adata)
    x = [Variable(model) for _ = 1:n]
    resid = @expression A' * x
    @test resid() == Adata' * x
end

@testset "tuple (https://github.com/tkoolen/Parametron.jl/issues/88)" begin
    A = [1, 2, 3, 4]
    @test(@expression(reshape(A, (2, 2))) == reshape(A, (2, 2)))
end

@testset "scale SVector" begin
    model = mock_model()
    x = Parameter(model, val=SVector(Variable(1), Variable(2)))
    expr = @expression x * 3
    @test expr() == SVector(3 * x.val[][1], 3 * x.val[][2])
    expr = @expression 3 * x
    @test expr() == SVector(x.val[][1] * 3, x.val[][2] * 3)
end

end # module


text:
function HarmonicOscillator(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return HarmonicOscillator((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function value(obj::HarmonicOscillator, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::HarmonicOscillator, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
# new Core

include("kpreso.jl") # for eename_clause() means a cycle dependency happend

# clause format
# Meta.parse("[x,y].[+P(ha=x,desu=f(x),-Q(ga=x,ni=y)]")
# as readclausesfromfile()

struct CForm2
 cid
 vars
 body
end

struct VForm2
 cid
 vars
end

struct LForm2
 lid
 body
end

struct STEPC
 rid
 leftp
 rightp
end

struct STEP
 rid
 leftp
 rightp
 sigma
 rename
end

mutable struct CORE
 name
 maxcid
 maxrid
 ldb
 cdb
 clmap
 lcmap
 allpsym
 level0
 proof
 proc
 trycnt
 succnt
end

function stringtoclause(cid, cls)
 vars = cls.args[1].args
 body = lit2kplit(cls.args[2].args[1].args)
# body = cls.args[2].args[1].args
 rcls = rename_clause(cid, vars, body)
 return CForm2(cid, rcls.vars, rcls.body)
end

function cform2ofclause(clss, initn=1)
 nclss=Dict()
 maxcid = cidof(1+initn-1)
 for cno in 1:length(clss)
  cid = cidof(cno+initn-1)
  maxcid = cid
  cls = stringtoclause(cid, clss[cno])
#  -F(x,y),-F(y,z),+F(x,z) is satisfiable... but not good to remove it.
#   so, may satisfible clause be removed or not ?
#  if satisfiable(cls.vars, cls.body); continue end
  nclss[cid] = cls
 end
 (maxcid, nclss)
end

function numberingliterals(ncls, lno0, ldb, lcmap, clmap, allpsym)
 cid = ncls.cid
 lids = []
 nlno0=0
 for lno in 1:length(ncls.body)
  cls=ncls.body[lno]
  if isa(cls, KPLiteral)
    lit = cls
  else
    lit = lit2kplit(ncls.body[lno])
  end
  lno = lno+lno0-1
  nlno0 = lno+1
  lid = lidof(lno)
  lcmap[lid] = cid
  push!(allpsym, lit.atom.Psym)
#  push!(allpsym, lit.args[2].args[1])
  push!(lids, lid)
  ldb[lid] = LForm2(lid, lit)
 end
 clmap[cid] = lids

 return nlno0, ldb, lcmap, clmap, allpsym
end

function createLDB(clss)
 ldb = Dict()
 lno = 1
 ltoc = Dict()
 ctol = Dict()
 allpsym = Set()
 for ncid in keys(clss)
  ncls=clss[ncid]
  lno, ldb, ltoc, ctol, allpsym = numberingliterals(ncls, lno, ldb, ltoc, ctol, allpsym)
 end
 return ldb, ltoc, ctol, allpsym
end

function vform2ofclause(cdb)
 vcl=Dict()
 for cid in keys(cdb)
   vcl[cid] = VForm2(cid, cdb[cid].vars)
 end
 vcl
end

function createcore(name, clss, proc)
 (maxcid, cdb) = cform2ofclause(clss)
 ldb, lcmap, clmap, allpsym =createLDB(cdb)
 cvdb = vform2ofclause(cdb)
 graph= Dict()
 CORE(name, [numof(maxcid)], [0], ldb, cvdb, clmap, lcmap, sort(collect(allpsym)), graph, Dict(), proc, [0], [0])
end

# proof operation
stepof(rid, core) = core.proof[rid]

function proofcof(rid, core)
 proof = []
 rids  = [rid]

 while !isempty(rids)
  rid = pop!(rids)

  if isrid(rid)
    astep=stepof(rid, core)
    push!(rids, cidof(astep.leftp, core))
    push!(rids, cidof(astep.rightp, core))
    push!(proof, STEPC(rid, cidof(astep.leftp, core), cidof(astep.rightp, core)))
  else
    continue
  end
 end
 proof
end

#core operation

cidof(num::Int) = Symbol(:C, num)
lidof(num::Int) = Symbol(:L, num)
ridof(num::Int) = Symbol(:R, num)

iscid(sym) = ('C' == string(sym)[1])
isrid(sym) = ('R' == string(sym)[1])
islid(sym) = ('L' == string(sym)[1])

function newlid(xid, lid)
  Symbol(lid,xid)
end

function newrid(core)
 core.maxrid[1] += 1
 ridof(core.maxrid[1])
end

function origof(xid)
 pix = findfirst(isequal(SEPSYM), string(xid))
 pix == nothing && return xid
 ix = findfirst(isequal('_'), string(xid))
 ix == nothing && return xid
 return Symbol(string(xid)[1:(ix -1 )])
end

function origtermof(term)
 if isa(term, Symbol); return origof(term) end
 if isa(term, Expr)
  for aix in 2:length(term.args)
   term.args[aix] = origtermof(term.args[aix])
  end
 end
 term
end

numof(xid) = Meta.parse(string(xid)[2:end])

function varsof(cid, core)
 core.cdb[cid].vars
end

function bodyof(cid, core)
 clmap=core.clmap[cid]
 body = []
 for lid in core.clmap[cid]
  push!(body, core.ldb[lid].body)
 end
 body
end

function clause2of(cid, core)
 CForm2(cid, varsof(cid, core), bodyof(cid, core))
end

function lidsof(cid, core)
 core.clmap[cid]
end

function cidof(lid, core)
 core.lcmap[lid]
end

function literalof(lid, core)
 core.ldb[lid]
end

function literalsof(lids, core)
 map(lid->literalof(lid, core).body, lids)
end

function evalproc(proc)
  for p in proc
    eval(Meta.parse(p))
  end
end

function evalcore(core)
  evalproc(core.proc)
end

#### analyzer of cnf
struct COREINFO
 vsyms
 csyms
 fsyms
 psyms

end

function analyze_term(vars, term)
 v = []
 c = []
 f = []
 if issym(term)
  isvar(term, vars) && push!(v, term)
  isconst(term, vars) && push!(c, term)
 else
  push!(f, term.args[1])
  nv, nc, nf = analyze_args(vars, term.args[2:end])
  append!(v, nv)
  append!(c, nc)
  append!(f, nf)
 end
 return v,c,f
end

function analyze_args(vars, args)
 v = []
 c = []
 f = []
 for term in args
   nv, nc, nf = analyze_term(vars, term)
   append!(v, nv)
   append!(c, nc)
   append!(f, nf)
 end
 return v,c,f
end

function analyze_lit(vars, atom)
 (v,c,f) = analyze_args(vars, atom.args[2:end])
 return v, c, f, atom.args[1]
end

function analyze_sym(core)
 ldb = core.ldb
 vpool = []
 cpool = []
 fpool = []
 ppool = []

 for lid in keys(ldb)
  atom = ldb[lid].body.args[2] ## remove sign of lit
  if @isdefined(atom.args[1]); continue end
  cid = cidof(lid, core)
  vars = varsof(cid, core)

  (v,c,f,p) = analyze_lit(vars, atom)
  append!(vpool, v)
  append!(cpool, c)
  append!(fpool, f)
  push!(ppool, p)
 end
 return COREINFO(Set(vpool), Set(cpool), Set(fpool), Set(ppool))
end

#### resolvent
## proof information in resolvent

struct RForm2
 rid
 vars
 left
 right
 sigma
 body
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# # Van der Pol oscillator
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/VanDerPol.ipynb)
#
#md # !!! note "Overview"
#md #     System type: polynomial continuous system\
#md #     State dimension: 2\
#md #     Application domain: Nonlinear physics
#
# The Van der Pol oscillator was introduced by the Dutch physicist Balthasar van der Pol.
# This is a famous model, typically investigated in the study of nonlinear dynamics.
# The model presents non-conservative oscillations with non-linear damping.
# In the past, it has been of relevance in several practical problems of engineering such as
# circuits containing vacuum tubes. For more information on the model see the wikipedia article
# [Van der Pol oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator).

# In this notebook we compute the evolution of the limit cycle in the phase plane
# for a set of initial conditions. For such set, we consider a *safety condition*
# originally from [^ARCHCOMP19LIN], that aims at verifying that there is no solution
# starting from the initial set that exceeds a prescribed upper bound on the
# velocity $y(t) := x'(t)$ at all times in the given time span. Moreover, we also
# illustrate the computation of an invariant set using the obtained flowipe.
# Finally, we study the limit cycle à la Poincaré, by defining a function that computes
# the cross section of the flowpipe at each revolution. Interestingly, this method
# gives an algorithmic proof that the safety condition obtained previously is actually
# verified at all times, i.e. over the unbounded time horizon $[0, \infty)$.

# ## Dynamics
#
# The dynamics of the Van der Pol oscillator are described by the following ODE with
# two variables:
#
# ```math
# \begin{aligned}
#   \dot{x} &= y \\
#   \dot{y} &= \mu (1 - x^2) y - x
# \end{aligned}
# ```
# The system has a stable limit cycle. Such limit cycle becomes increasingly sharp
# for higher values of ``μ``. Here we consider the parameter ``μ = 1``.

#!jl using ReachabilityAnalysis, Plots

@taylorize function vanderpol!(dx, x, params, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

#md # !!! tip "Performance tip"
#md #     Advanced users may want to consider the section [Some common gotchas](@ref)
#md #     to improve the code that is generated by `@taylorize` by re-writing
#md #     the right-hand side of `dx[2]` with the use of auxiliary variables.

#jl function vanderpol()
#jl    X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])
#jl    ivp = @ivp(x' = vanderpol!(x), dim: 2, x(0) ∈ X0)
#jl    tspan = (0.0, 5.0)
#jl    return ivp, tspan
#jl end

# ## Safety verification
#
# We set the initial condition ``x(0) ∈ [1.25, 1.55]``, ``y(0) ∈ [2.35,2.45]``.
# The *unsafe set* is given by ``y ≥ 2.75`` for a time span ``[0, 7]``.
# In other words, we would like to prove that there doesn't exist a solution of
# the model with a ``y`` value which is greater than 2.75, for any initial condition
# on the given domain. The time horizon of ``T = 7`` is chosen such that the oscillator
# can do at least one full cycle.

# We proceed by defining the initial conditions as a hyperrectangular set according
# to the problem's specifications. Then we build the initial-value problem and pass
# it to the `solve` function. We specify using `TMJets` algorithm with default options.

X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45]) #!jl
prob = @ivp(x' = vanderpol!(x), dim=2, x(0) ∈ X0) #!jl
sol = solve(prob, T=7.0, alg=TMJets(abstol=1e-12)); #!jl

# For further computations, it is convenient to work with a [zonotopic](https://juliareach.github.io/ReachabilityAnalysis.jl/dev/man/basics/#Zonotopes-1)
# overapproximation of the flowpipe.

solz = overapproximate(sol, Zonotope); #!jl

# The maximum value of variable ``y`` is obtained by computing the [support function](https://juliareach.github.io/ReachabilityAnalysis.jl/dev/man/basics/#Representing-sets-with-support-functions-1)
# of the flowpipe along direction ``[0, 1]``:

ρ([0.0, 1.0], solz) #!jl

# That shows that the property is satisfied. Below we plot the flowpipe in the
# x-y plane, together with the horizontal line ``y = 2.75``.

plot(solz, vars=(1, 2), lw=0.2, xlims=(-2.5, 2.5), xlab="x", ylab="y") #!jl
plot!(x -> 2.75, color=:red, lab="y = 2.75", style=:dash, legend=:bottomright) #!jl

# We can also plot the state variables ``x(t)`` and ``y(t)`` as a function of time
# (recall that `0` in `vars` is used to denote the time variable):

plot(solz, vars=(0, 1), lw=0.2, xlab="t", ylab="x") #!jl

#-

plot(solz, vars=(0, 2), lw=0.2,  xlab="t", ylab="y") #!jl

# ## Invariant Set
#
# We can use the reachability result to examine an invariant of the system. In
# other words, we can algorithmically prove that the flowpipe re-enters from where
# it started after giving one loop, using inclusion checks.

plot(solz, vars=(1, 2), lw=0.2, xlims=(0.0, 2.5), ylims=(1.6, 2.8), xlab="x", ylab="y") #!jl
plot!(X0, color=:orange, lab="X0") #!jl
plot!(solz[1:13], vars=(1, 2), color=:green, lw=1.0, alpha=0.5, lab="F[1:13]") #!jl
plot!(solz[388], vars=(1, 2), color=:red, lw=1.0, alpha=0.6, lab="F[388]") #!jl

# It is seen that the reach-set corresponding to the time-span

tspan(solz[388])  #!jl

# is included in the set union ``F[1] \cup \cdots \cup F[13]`` of previously
# computed reach-sets. Notice that all future trajectories starting from
# the 388-th reach-set are already covered by the flowpipe. Therefore, we can claim
# that an *invariant set* was found.

# ## Limit cycle

# To examine the limit cycle we can intersect a line segment perpendicular to the
# flowpipe, that will allow us to get a cross-section of the sets in order to prove
# that after one cycle the intersection segment actually shrinks. This approach
# is similar to the method of [Poincaré sections](https://en.wikipedia.org/wiki/Poincar%C3%A9_map).

line = LineSegment([1, 2.], [2., 2.5]) #!jl
plot(solz, vars=(1, 2), lw=0.2, xlims=(0.0, 2.5), ylims=(1.6, 2.8), xlab="x", ylab="y") #!jl
plot!(X0, color=:orange, lab="X0") #!jl
plot!(solz[1:13], vars=(1, 2), color=:green, lw=1.0, alpha=0.5, lab="F[1:13]") #!jl
plot!(solz[388], vars=(1, 2), color=:red, lw=1.0, alpha=0.6, lab="F[388]") #!jl
plot!(line, lw=2.0) #!jl

# Then we can define a function to get the cross section of the flowpipe. The
# function needs the flowpipe, a line segment that cuts the flowpipe and the
# indices of the subsets to cut.

function cross_section(line::LineSegment, fp, idx) #!jl
    p = VPolygon() #!jl
    for i in idx #!jl
        x = intersection(line, set(fp[i])) #!jl
        if !isempty(x)  #!jl
            p = convex_hull(p, x) #!jl
        end  #!jl
    end #!jl
    vl = vertices_list(p) #!jl
    @assert length(vl) == 2 #!jl
    return LineSegment(vl[1], vl[2]) #!jl
end #!jl

# Then we can get the cross section of the first five sets and the last set,
# calling them `i1` and `i2` respectively.

ifirst = cross_section(line, solz, 1:13) #!jl
ilast = cross_section(line, solz, [388]) #!jl

# We can also calculate the length of each cross section, remember that the
# system is 2D, so the cross section will be a line segment.

lfirst = norm(ifirst.q - ifirst.p) #!jl
#-
llast = norm(ilast.q - ilast.p); #!jl

#-

plot(ifirst, lw=3.0, alpha=1.0, label="First subsets", legend=:bottomright) #!jl
plot!(ilast, lw=5.0, alpha=1.0, label="Last subset") #!jl

#-

# The inclusion check succeeds:

ilast ⊆ ifirst #!jl

# We can see, the cross section of the last subset is a subset of the first few
# sets, thus, the cycle will continue, presumably getting smaller each revolution.

# ## References

# [^ARCHCOMP19LIN]: Althoff, M., Bak, S., Forets, M., Frehse, G., Kochdumper, N., Ray, R., ... & Schupp, S. (2019, May). ARCH-COMP19 Category Report: Continuous and Hybrid Systems with Linear Continuous Dynamics. In ARCH@ CPSIoTWeek (pp. 14-40).


text:
using NBodyIPs, JuLIP, Test, IPFitting, DataFrames
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB
using NBodyIPs: BondLengthDesc, BondAngleDesc
using NBodyIPs.Polys: NBPoly

const Lsq = IPFitting.Lsq
const Err = IPFitting.Errors

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.25*r0, 100, calc)

TRANSFORM = "exp( - 2 * (r/$r0 - 1.5) )"

rcut2 = cutoff(calc)*1.6
CUTOFF2 = "(:cos, $(rcut2-1), $rcut2)"
D2 = BondLengthDesc(TRANSFORM, CUTOFF2)

rcut3 = cutoff(calc)*1.6
CUTOFF3 = "(:cos, $(rcut3-1), $rcut3)"
D3 = BondLengthDesc(TRANSFORM, CUTOFF3)

cw = Dict("rand" => 1.0)
dw   = Dict("E" => 100.0, "F" => 1.0)

##

# create a monomial database
B = nbpolys(3, D3, 8)
@show length(B)
db = LsqDB("", B, data)

# orthogonalise the basis
Bo = Lsq.onb(db, E0 = 0.0, configweights = cw, dataweights = dw)
# extract the pure 3-body terms
Bo = [ b.components[1] for b in Bo ]

# checkout how "nice" the coefficients are
@show [ norm(b.c) for b in Bo ]

# construct a new database, based on the data-ONB
db_onb = LsqDB("", Bo, data)

# fit with ONB
IP_onb, errs = Lsq.lsqfit( db_onb,
      E0 = 0.0,
      configweights = Dict("rand" => 1.0),
      dataweights   = Dict("E" => 100.0, "F" => 1.0)
   )

# fit the monomial basis
IP, errs = Lsq.lsqfit( db,
      E0 = 0.0,
      configweights = Dict("rand" => 1.0),
      dataweights   = Dict("E" => 100.0, "F" => 1.0)
   )

# check for errors
@show maximum( abs(energy(IP_onb, d.at) - energy(IP, d.at))  for d in data  )


text:
"""
[Based on arXiv:math/9903025v2 by Kenyon, Propp, Wilson (2000)]
Generates samples from the dimer model on a 2-d rectangular lattice.
Also supports loop erased random walks and Wilson's algorithm on an
arbitrary graph.

Authors: Samuel S. Watson, Doruk Efe Gökmen
Date: 17/01/2020

TODO: implement the dataformat used in the MKJ-ZR paper.
"""

__precompile__(true)

module Dimers

export Wilson,
       LERW,
       gridgraph,
       rotate,
       flatten,
       midpoint,
       dimersample,
       dimerheight,
       dimerSitesBonds

import Graphs

function LERW(Γ::Graphs.AbstractGraph,
                   startingvertex::Int64,
                   roots::Array{Bool,1};
                   maxiter::Integer=10^7)

    X = [startingvertex]
    cntr = 1
    while ~(roots[X[end]])
        neighbors = Γ.adjlist[X[end]]
        push!(X,Graphs.vertex_index(neighbors[rand(1:length(neighbors))],Γ))
        i = 1
        while X[i] != X[end]
            i += 1
        end
        X = X[1:i]
        cntr += 1
        if cntr >= maxiter
            error("Maximum iterations hit in LERW")
        end
    end
    vert = Graphs.vertices(Γ)
    return vert[X]
end

function Wilson(Γ::Graphs.AbstractGraph,roots::Array{Bool,1})

    if length(Graphs.connected_components(Γ)) > 1
        error("Graph not connected")
    end

    maxiter = length(Graphs.vertices(Γ))
    cntr = 1
    UST = Graphs.adjlist(typeof(Graphs.vertices(Γ)[1]))
    for v in Graphs.vertices(Γ)
        Graphs.add_vertex!(UST,v)
    end
    discovered = roots
    while ~all(discovered)
        i = 1
        while discovered[i]
            i += 1
        end
        lerw = LERW(Γ,i,discovered)
        for k=1:length(lerw)-1
            Graphs.add_edge!(UST,lerw[k+1],lerw[k])
            discovered[Graphs.vertex_index(lerw[k],Γ)] = true
        end
        cntr += 1
        if cntr >= maxiter
            error("Something's gone wrong with Wilson's algorithm")
            break
        end
    end
    return UST
end

rotate(segment::Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}) =
(div(segment[1][1] + segment[2][1] - segment[1][2] + segment[2][2],2),
div(segment[1][1] - segment[2][1] + segment[1][2] + segment[2][2],2)),
(div(segment[1][1] + segment[2][1] + segment[1][2] - segment[2][2],2),
div(-segment[1][1] + segment[2][1] + segment[1][2] + segment[2][2],2))

flatten(e::Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}) = vcat(map(x->vcat(x...),vcat(e...))...)

midpoint(point1::Tuple{Int64,Int64},point2::Tuple{Int64,Int64}) = (div(point1[1] + point2[1],2),div(point1[2] + point2[2],2))

function gridgraph(n::Int64)

    Γ = Graphs.adjlist(Tuple{Int64,Int64},is_directed=false)

    for i=1:n
        for j=1:n
            Graphs.add_vertex!(Γ,(i,j))
        end
    end

    gridedges = Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}[]

    for i=1:n
        for j=1:n
            if i < n
                Graphs.add_edge!(Γ,(i,j),(i+1,j))
            end
            if j < n
                Graphs.add_edge!(Γ,(i,j),(i,j+1))
            end
        end
    end

    return Γ
end

function dimersample(m::Int64,n::Int64)

    function add_edge_and_continue(vertex::Tuple{Int64,Int64},prev_vertex::Tuple{Int64,Int64})
        Graphs.add_edge!(dualtree_ordered,prev_vertex,vertex)
        if (vertex[1]+2,vertex[2]) != prev_vertex && (vertex[1]+2,vertex[2]) in Graphs.out_neighbors(vertex,dualtree)
            add_edge_and_continue((vertex[1]+2,vertex[2]),vertex)
        end
        if (vertex[1]-2,vertex[2]) != prev_vertex && (vertex[1]-2,vertex[2]) in Graphs.out_neighbors(vertex,dualtree)
            add_edge_and_continue((vertex[1]-2,vertex[2]),vertex)
        end
        if (vertex[1],vertex[2]+2) != prev_vertex && (vertex[1],vertex[2]+2) in Graphs.out_neighbors(vertex,dualtree)
            add_edge_and_continue((vertex[1],vertex[2]+2),vertex)
        end
        if (vertex[1],vertex[2]-2) != prev_vertex && (vertex[1],vertex[2]-2) in Graphs.out_neighbors(vertex,dualtree)
            add_edge_and_continue((vertex[1],vertex[2]-2),vertex)
        end
    end

    Γ = Graphs.adjlist(Tuple{Int64,Int64},is_directed=false)

    for i=1:m+1
        for j=1:n+1
            if (i,j) != (m+1,n+1)
                Graphs.add_vertex!(Γ,(2i-1,2j-1))
            end
        end
    end

    primaledges = Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}[]

    for i=1:m
        for j=1:n
            push!(primaledges,((2i-1,2j-1),(2i-1+2,2j-1)))
            push!(primaledges,((2i-1,2j-1),(2i-1,2j-1+2)))
        end
    end

    for e in primaledges
        Graphs.add_edge!(Γ,e...)
    end

    dualtree = Graphs.adjlist(Tuple{Int64,Int64},is_directed=false)

    for i=0:m
        for j=0:n
            if (i,j) != (0,0)
                Graphs.add_vertex!(dualtree,(2i,2j))
            end
        end
    end

    roots = [false for i=1:length(Graphs.vertices(Γ))]
    dualroots = [false for i=1:length(Graphs.vertices(dualtree))]

    for i=1:m
        roots[Graphs.vertex_index((2i-1,2n+1),Γ)] = true
        dualroots[Graphs.vertex_index((2i,0),dualtree)] = true
    end

    for j=1:n
        roots[Graphs.vertex_index((2m+1,2j-1),Γ)] = true
        dualroots[Graphs.vertex_index((0,2j),dualtree)] = true
    end

    UST = Wilson(Γ,roots)

    for e in primaledges
        if ~(e[2] in Graphs.out_neighbors(e[1],UST) || e[1] in Graphs.out_neighbors(e[2],UST))
            newedge = rotate(e)
            Graphs.add_edge!(dualtree,newedge...)
        end
    end

    dualtree_ordered = Graphs.inclist(Tuple{Int64,Int64},is_directed=true)

    for v in Graphs.vertices(dualtree)
        Graphs.add_vertex!(dualtree_ordered,v)
    end

    for k=2:2:2n
        if (2,k) in Graphs.out_neighbors((0,k),dualtree) || (0,k) in Graphs.out_neighbors((2,k),dualtree)
            add_edge_and_continue((2,k),(0,k))
        end
        if (k,0) in Graphs.out_neighbors((k,2),dualtree) || (k,2) in Graphs.out_neighbors((k,0),dualtree)
            add_edge_and_continue((k,2),(k,0))
        end
    end

    dimergraph = Graphs.adjlist(typeof(Graphs.vertices(Γ)[1]))
    dimer_vertices = Tuple{Int64,Int64}[]
    dimer_edges = Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}[]

    for i=1:2m
        for j=1:2n
            push!(dimer_vertices,(i,j))
            if i < 2m
                push!(dimer_edges,((i,j),(i+1,j)))
            end
            if j < 2n
                push!(dimer_edges,((i,j),(i,j+1)))
            end
            Graphs.add_vertex!(dimergraph,(i,j))
        end
    end

    for v in Graphs.vertices(UST)
        for w in Graphs.out_neighbors(v,UST)
            Graphs.add_edge!(dimergraph,midpoint(v,w),w)
        end
    end

    for v in Graphs.vertices(dualtree_ordered)
        for w in Graphs.out_neighbors(v,dualtree_ordered)
            Graphs.add_edge!(dimergraph,midpoint(v,w),w)
        end
    end

    return dimergraph

end

dimersample(n::Integer) = dimersample(n,n)

"""
function drawgraph{V,E}(Γ::Graphs.AbstractGraph{V,E};
                         pointsize=0.002,
                         linewidth=1.0)

    all_points = Graphics2D.GraphicElement[]
    all_edges = Graphics2D.GraphicElement[]

    for v in Graphs.vertices(Γ)
        push!(all_points,Graphics2D.Point(v...;pointsize=pointsize))
        for w in Graphs.out_neighbors(v,Γ)
            push!(all_edges,Graphics2D.Line([v[1] v[2]; w[1] w[2]];linewidth=linewidth))
        end
    end

    return [all_points; all_edges]
end
"""

function dimerSitesBonds(N_white::Int64, corr_spin::Bool)
    g = dimersample(N_white)

    A = zeros(Int8, 4N_white, 4N_white)

    for i=1:N_white
        for j=1:N_white
            k = 2i-1
            l = 2j
            delta1 = [t for t in Graphs.out_neighbors((k, l), g)[1]] - [k, l]

            m = 2i
            n = 2j-1
            delta2 = [t for t in Graphs.out_neighbors((m, n), g)[1]] - [m, n]

            A[([2k-1, 2l-1] + delta1)[1], ([2k-1, 2l-1] + delta1)[2]] = 1
            A[([2m-1, 2n-1] + delta2)[1], ([2m-1, 2n-1] + delta2)[2]] = 1

            for p=0:1
                for r=0:1
                    A[2(i+p*N_white)-1, 2(j+r*N_white)-1] = rand(0:1)
                    if corr_spin
                        A[2(i+p*N_white), 2(j+r*N_white)] = A[2(i+p*N_white)-1, 2(j+r*N_white)-1]
                    else
                        A[2(i+p*N_white), 2(j+r*N_white)] = rand(0:1)
                    end
                end
            end
        end
    end

    return A
end

function dimerheight(dimergraph::Graphs.AbstractGraph)

    all_edges = Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}[]

    m = div(maximum(map(x->x[1],Graphs.vertices(dimergraph))),2)
    n = div(maximum(map(x->x[2],Graphs.vertices(dimergraph))),2)

    for i=1:2m
        for j=1:2n
            if i < 2m
                push!(all_edges,((i,j),(i+1,j)))
            end
            if j < 2n
                push!(all_edges,((i,j),(i,j+1)))
            end
        end
    end

    down_edges = zeros(Bool,2m+1,2n+1)

    for i=1:2m-1
        for j=1:2n
            if (i,j) in Graphs.out_neighbors((i+1,j),dimergraph) || (i+1,j) in Graphs.out_neighbors((i,j),dimergraph)
                down_edges[i,j] = true
            end
        end
    end

    h = zeros(Int64,2n+1,2n+1)

    for j=2:2:2n+1
        h[1,j] = 1
    end

    for i=2:2n+1
        h[i,1] = isodd(i) ? 0 : -1
        for j=2:2n+1
            if down_edges[i-1,j-1]
                h[i,j] = h[i-1,j] + (-1)^(i+j)
            else
                h[i,j] = h[i,j-1] + (-1)^(i+j+1)
            end
        end
    end

    return h
end

end # module


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ a2a09c0b-d3ec-4b46-9bca-7c69b794f8ce
md"# Solutions to exercises"

# ╔═╡ 40956165-26d0-4d61-b795-033f6f2a0674
md"## Setup"

# ╔═╡ 197fd00e-9068-46ea-935f-ddf6a6bfbbdd
md"Inspect Julia version:"

# ╔═╡ f6d4f8c4-e441-45c4-af2a-25235e544a31
VERSION

# ╔═╡ 45740c4d-b789-45dc-8af5-148d95ea2900
md"The following instantiates a package environment."

# ╔═╡ 42b0f1e1-16c9-4238-a6bf-47194d7e8e12
md"""
The package environment has been created using **Julia 1.6** and may not
instantiate properly for other Julia versions.
"""

# ╔═╡ 754d60ac-4cb0-443a-828a-4cc485149963
begin
  using Pkg
  Pkg.activate(@__DIR__)
  Pkg.instantiate()
end

# ╔═╡ 19dae74b-24bd-428f-9e54-3f7b3ca87ecb
md"## Exercise 1"

# ╔═╡ 3af8467f-4988-44dc-ba1f-6363f43ec697
md"Write a function named `total` that adds the elements of its vector input."

# ╔═╡ cd6e9fce-e54a-48c3-95ea-2955abd45275
md"### Solution"

# ╔═╡ 992f71c3-12d6-4ab6-b1b3-dd47e367ba54
function total(v)
    sum = 0.0 # better is zero(v)
    for i in 1:length(v)
        sum = sum + v[i]
    end
    return sum
end

# ╔═╡ 61304232-d701-4218-8d7e-8e449afd1c48
total(1:10)

# ╔═╡ 18960c40-199a-48fa-a949-7f3bd292fe31
md"The built-in function is called `sum`."

# ╔═╡ b04e2fdc-8f1a-4964-8514-99938a2932db
md"## Exercise 2"

# ╔═╡ 860bb453-c2f1-446d-a0de-1721c1bc2df2
md"""
Generate a 1000 random samples from the standard normal
distribution. Create a second such sample, and add the two samples
point-wise.  Compute the (sample) mean and variance of the combined
samples. In the same plot, show a frequency-normalized histogram of
the combined samples and a plot of the pdf for normal distribution
with zero mean and variance `2`.
"""

# ╔═╡ cd6e9fce-e54a-48c3-bca8-7eec7950fd82
md"### Solution"

# ╔═╡ fe0bd76b-d2bb-4f68-9873-1b1430e635cc
begin
  using Distributions, Plots, Statistics

  samples1 = randn(1000); # or rand(Normal(), 1000)
  samples2 = randn(1000);

  samples = samples1 .+ samples2;

  mu = mean(samples)
  var = std(samples)^2

  @show mu var
end

# ╔═╡ 37a8bacf-d991-42ac-b43e-d0ebe87da176
begin
  d = Normal(0, sqrt(2))
  f(x) = pdf(d, x)

  xs = -5:(0.1):5
  ys = f.(xs);

  using Plots

  plt = histogram(samples, normalize=true)
  plot!(xs, ys)
  plt
end

# ╔═╡ 4f939d4a-e802-4b7d-9009-4cfb2012998f
md"## Exercise 3"

# ╔═╡ 63f410ea-c37a-4761-abd3-cb77d7a79683
md"The following shows that named tuples share some behaviour with dictionaries:"

# ╔═╡ cb2c42bf-21b5-4e04-879e-dc128f3db7b3
begin
  t = (x = 1, y = "cat", z = 4.5)
  keys(t)
end

# ╔═╡ b8996fa3-c046-4fb1-a367-9aa3c6cf34d0
t[:y]

# ╔═╡ 6584c100-ffd6-4fa1-bf32-94657e65284e
md"""
Write a function called `dict` that converts a named tuple to an
actual dictionary. You can create an empty dictionary using `Dict()`.
"""

# ╔═╡ cd6e9fce-e54a-48c3-9afd-65a5b5facac9
md"### Solution"

# ╔═╡ 834576cc-c6dc-4d3e-b6c8-54f5ed90a964
function dict(t)
    d = Dict()
    for k in keys(t)
        d[k] = t[k]
    end
    return d
end

# ╔═╡ dfdc3d18-4274-491d-9292-fe822525fc77
dict(t)

# ╔═╡ a8888dc5-ef44-4ee7-ae5e-4b83dcbc9675
md"Or alternatively:"

# ╔═╡ 5b64ff6d-923b-4eab-8a29-32f11452654a
dict(t) = Dict(k => t[k] for k in keys(t))

# ╔═╡ 135dac9b-0bd9-4e1d-a5f1-151fcbe229a2
md"""
---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
"""

# ╔═╡ Cell order:
# ╟─a2a09c0b-d3ec-4b46-9bca-7c69b794f8ce
# ╟─40956165-26d0-4d61-b795-033f6f2a0674
# ╟─197fd00e-9068-46ea-935f-ddf6a6bfbbdd
# ╠═f6d4f8c4-e441-45c4-af2a-25235e544a31
# ╟─45740c4d-b789-45dc-8af5-148d95ea2900
# ╟─42b0f1e1-16c9-4238-a6bf-47194d7e8e12
# ╠═754d60ac-4cb0-443a-828a-4cc485149963
# ╟─19dae74b-24bd-428f-9e54-3f7b3ca87ecb
# ╟─3af8467f-4988-44dc-ba1f-6363f43ec697
# ╟─cd6e9fce-e54a-48c3-95ea-2955abd45275
# ╠═992f71c3-12d6-4ab6-b1b3-dd47e367ba54
# ╠═61304232-d701-4218-8d7e-8e449afd1c48
# ╟─18960c40-199a-48fa-a949-7f3bd292fe31
# ╟─b04e2fdc-8f1a-4964-8514-99938a2932db
# ╟─860bb453-c2f1-446d-a0de-1721c1bc2df2
# ╟─cd6e9fce-e54a-48c3-bca8-7eec7950fd82
# ╠═fe0bd76b-d2bb-4f68-9873-1b1430e635cc
# ╠═37a8bacf-d991-42ac-b43e-d0ebe87da176
# ╟─4f939d4a-e802-4b7d-9009-4cfb2012998f
# ╟─63f410ea-c37a-4761-abd3-cb77d7a79683
# ╠═cb2c42bf-21b5-4e04-879e-dc128f3db7b3
# ╠═b8996fa3-c046-4fb1-a367-9aa3c6cf34d0
# ╟─6584c100-ffd6-4fa1-bf32-94657e65284e
# ╟─cd6e9fce-e54a-48c3-9afd-65a5b5facac9
# ╠═834576cc-c6dc-4d3e-b6c8-54f5ed90a964
# ╠═dfdc3d18-4274-491d-9292-fe822525fc77
# ╟─a8888dc5-ef44-4ee7-ae5e-4b83dcbc9675
# ╠═5b64ff6d-923b-4eab-8a29-32f11452654a
# ╟─135dac9b-0bd9-4e1d-a5f1-151fcbe229a2


text:
using Test

using Distributions, LinearAlgebra, Statistics
using EnKF

include("inflation.jl")

include("state.jl")

include("initial.jl")

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
# Automatically generated using Clang.jl

const XC_UNPOLARIZED = 1
const XC_POLARIZED = 2
const XC_NON_RELATIVISTIC = 0
const XC_RELATIVISTIC = 1
const XC_EXCHANGE = 0
const XC_CORRELATION = 1
const XC_EXCHANGE_CORRELATION = 2
const XC_KINETIC = 3
const XC_FAMILY_UNKNOWN = -1
const XC_FAMILY_LDA = 1
const XC_FAMILY_GGA = 2
const XC_FAMILY_MGGA = 4
const XC_FAMILY_LCA = 8
const XC_FAMILY_OEP = 16
const XC_FAMILY_HYB_GGA = 32
const XC_FAMILY_HYB_MGGA = 64
const XC_FAMILY_HYB_LDA = 128
const XC_FLAGS_HAVE_EXC = 1 << 0
const XC_FLAGS_HAVE_VXC = 1 << 1
const XC_FLAGS_HAVE_FXC = 1 << 2
const XC_FLAGS_HAVE_KXC = 1 << 3
const XC_FLAGS_HAVE_LXC = 1 << 4
const XC_FLAGS_1D = 1 << 5
const XC_FLAGS_2D = 1 << 6
const XC_FLAGS_3D = 1 << 7
const XC_FLAGS_HYB_CAM = 1 << 8
const XC_FLAGS_HYB_CAMY = 1 << 9
const XC_FLAGS_VV10 = 1 << 10
const XC_FLAGS_HYB_LC = 1 << 11
const XC_FLAGS_HYB_LCY = 1 << 12
const XC_FLAGS_STABLE = 1 << 13
const XC_FLAGS_DEVELOPMENT = 1 << 14
const XC_FLAGS_NEEDS_LAPLACIAN = 1 << 15
const XC_FLAGS_HAVE_ALL = (((XC_FLAGS_HAVE_EXC | XC_FLAGS_HAVE_VXC) | XC_FLAGS_HAVE_FXC) | XC_FLAGS_HAVE_KXC) | XC_FLAGS_HAVE_LXC
const XC_EXT_PARAMS_DEFAULT = -999998888
const XC_TAU_EXPLICIT = 0
const XC_TAU_EXPANSION = 1
const XC_MAX_REFERENCES = 5

struct func_reference_type
    ref::Cstring
    doi::Cstring
    bibtex::Cstring
end

struct func_params_type
    n::Cint
    names::Ptr{Cstring}
    descriptions::Ptr{Cstring}
    values::Ptr{Cdouble}
    set::Ptr{Cvoid}
end

struct xc_func_info_type
    number::Cint
    kind::Cint
    name::Cstring
    family::Cint
    refs::NTuple{5, Ptr{func_reference_type}}
    flags::Cint
    dens_threshold::Cdouble
    ext_params::func_params_type
    init::Ptr{Cvoid}
    _end::Ptr{Cvoid}
    lda::Ptr{Cvoid}
    gga::Ptr{Cvoid}
    mgga::Ptr{Cvoid}
end

struct xc_dimensions
    rho::Cint
    sigma::Cint
    lapl::Cint
    tau::Cint
    zk::Cint
    vrho::Cint
    vsigma::Cint
    vlapl::Cint
    vtau::Cint
    v2rho2::Cint
    v2rhosigma::Cint
    v2rholapl::Cint
    v2rhotau::Cint
    v2sigma2::Cint
    v2sigmalapl::Cint
    v2sigmatau::Cint
    v2lapl2::Cint
    v2lapltau::Cint
    v2tau2::Cint
    v3rho3::Cint
    v3rho2sigma::Cint
    v3rho2lapl::Cint
    v3rho2tau::Cint
    v3rhosigma2::Cint
    v3rhosigmalapl::Cint
    v3rhosigmatau::Cint
    v3rholapl2::Cint
    v3rholapltau::Cint
    v3rhotau2::Cint
    v3sigma3::Cint
    v3sigma2lapl::Cint
    v3sigma2tau::Cint
    v3sigmalapl2::Cint
    v3sigmalapltau::Cint
    v3sigmatau2::Cint
    v3lapl3::Cint
    v3lapl2tau::Cint
    v3lapltau2::Cint
    v3tau3::Cint
    v4rho4::Cint
    v4rho3sigma::Cint
    v4rho3lapl::Cint
    v4rho3tau::Cint
    v4rho2sigma2::Cint
    v4rho2sigmalapl::Cint
    v4rho2sigmatau::Cint
    v4rho2lapl2::Cint
    v4rho2lapltau::Cint
    v4rho2tau2::Cint
    v4rhosigma3::Cint
    v4rhosigma2lapl::Cint
    v4rhosigma2tau::Cint
    v4rhosigmalapl2::Cint
    v4rhosigmalapltau::Cint
    v4rhosigmatau2::Cint
    v4rholapl3::Cint
    v4rholapl2tau::Cint
    v4rholapltau2::Cint
    v4rhotau3::Cint
    v4sigma4::Cint
    v4sigma3lapl::Cint
    v4sigma3tau::Cint
    v4sigma2lapl2::Cint
    v4sigma2lapltau::Cint
    v4sigma2tau2::Cint
    v4sigmalapl3::Cint
    v4sigmalapl2tau::Cint
    v4sigmalapltau2::Cint
    v4sigmatau3::Cint
    v4lapl4::Cint
    v4lapl3tau::Cint
    v4lapl2tau2::Cint
    v4lapltau3::Cint
    v4tau4::Cint
end

struct xc_func_type
    info::Ptr{xc_func_info_type}
    nspin::Cint
    n_func_aux::Cint
    func_aux::Ptr{Ptr{xc_func_type}}
    mix_coef::Ptr{Cdouble}
    cam_omega::Cdouble
    cam_alpha::Cdouble
    cam_beta::Cdouble
    nlc_b::Cdouble
    nlc_C::Cdouble
    dim::xc_dimensions
    params::Ptr{Cvoid}
    dens_threshold::Cdouble
    zeta_threshold::Cdouble
    sigma_threshold::Cdouble
    tau_threshold::Cdouble
end


text:
using Documenter
using SuiteSparseGraphBLAS
using LinearAlgebra
using SparseArrays
makedocs(
    modules = [SuiteSparseGraphBLAS],
    sitename="SuiteSparseGraphBLAS.jl",
    pages = [
        "Introduction" => "index.md",
        "Arrays" => "arrays.md",
        "Operations" => "operations.md",
        "Operators" => [
            "operators.md",
            "unaryops.md",
            "binaryops.md",
            "monoids.md",
            "semirings.md",
            "selectops.md",
            "udfs.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/JuliaSparse/SuiteSparseGraphBLAS.jl.git",
)


text:
#include "dmrg.jl"
function tdvptwosite(H::MPO, psi0::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    check_hascommoninds(siteinds, H, psi)
    check_hascommoninds(siteinds, H, psi')
    # Permute the indices to have a better memory layout
    # and minimize permutations
    H = permute(H, (linkind, siteinds, linkind))
    PH = ProjMPO(H)
    psi = copy(psi0)
    return tdvptwosite(PH, psi, dt, order; kwargs...)
  end

function tdvptwosite!(H::MPO, psi::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    check_hascommoninds(siteinds, H, psi)
    check_hascommoninds(siteinds, H, psi')
    # Permute the indices to have a better memory layout
    # and minimize permutations
    H = permute(H, (linkind, siteinds, linkind))
    PH = ProjMPO(H)
    return tdvptwosite(PH, psi, dt, order; kwargs...)
  end

function forward!(PH::ProjMPO, psi::MPS,dt::Number,kwargs)
    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")

    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)

    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )

    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)

    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)

    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)

    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end

    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end

    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end

    N = length(psi)
    for (b, ha) in sweepnext(N)
        if ha == 2
            return PH, psi
        end
        ##Taken care of in two steps below, assuming that before the sweep position was set
        #@timeit_debug timer "dmrg: position!" begin
        #  position!(PH, psi, b)
        #end

        @debug_check begin
            checkflux(psi)
            checkflux(PH)
        end

        @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
            phi = psi[b] * psi[b + 1]
        end

        @timeit_debug timer "tdvp2: expintegrator-2site" begin
            phi, propinfo = exponentiate(
            PH,
            dt,
            phi,
            ishermitian=ishermitian,
            tol=eigsolve_tol,
            krylovdim=eigsolve_krylovdim,
            #maxiter=eigsolve_maxiter,
            )
        end

        ortho = ha == 1 ? "left" : "right"

        @debug_check begin
            checkflux(phi)
        end

        @timeit_debug timer "dmrg: replacebond!" begin
            spec = replacebond!(
            psi,
            b,
            phi;
            maxdim=maxdim,
            cutoff=cutoff,
            #eigen_perturbation=drho,
            ortho=ortho,
            normalize=true,
            which_decomp=which_decomp,
            svd_alg=svd_alg,
            )
        end

        maxtruncerr = max(maxtruncerr, spec.truncerr)
        #@show inds(psi[b+1])

        if b+1<N
            position_leftonly!(PH, psi, b+1)
            phi=psi[b + 1]
            #@show inds(psi[b+1])
            #@show inds(phi)
            @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                psi[b+1], propinfo = exponentiate(
                PH,
                -dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #@show inds(vecs)
            position_rightonly!(PH, psi, b+1)
        else
            position!(PH,psi,b)
        end

        @debug_check begin
            checkflux(psi)
            checkflux(PH)
        end
    end
  end

function backward!(PH::ProjMPO, psi::MPS,dt::Number,kwargs)
    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")

    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)

    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )

    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)

    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)

    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)

    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end

    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end

    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end

    N = length(psi)
    for (b, ha) in sweepnext(N)
        if ha == 2
            ortho = ha == 1 ? "left" : "right"
            if b+1<N
                phi=psi[b + 1]

                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1] = vecs
                position_leftonly!(PH, psi, b)
            end
            #position_leftonly!(PH, psi, b)
            #position_rightonly!(PH, psi, b)

            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end

            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #energy::Number = vals[1]
            phi = vecs
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            maxtruncerr = max(maxtruncerr, spec.truncerr)

            if b > 1    #after completed backward-forward sweep we want to be left with 2-site PH
                position_rightonly!(PH, psi, b-1)
                #position_leftonly!(PH, psi, b-1)
            else
                return PH,psi
            end

        end
    end
  end

function tdvptwosite!(PH::ProjMPO, psi::MPS, dt::Number, order::Int; kwargs...)::Tuple{Number,MPS}
    if length(psi) == 1
      error(
        "`tdvp2` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.eigsolve`, etc.",
      )
    end

    @debug_check begin
      # Debug level checks
      # Enable with ITensors.enable_debug_checks()
      checkflux(psi)
      checkflux(PH)
    end

    which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
    #svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
    svd_alg::String = get(kwargs, :svd_alg, "recursive")

    obs = get(kwargs, :observer, NoObserver())
    outputlevel::Int = get(kwargs, :outputlevel, 1)

    write_when_maxdim_exceeds::Union{Int,Nothing} = get(
      kwargs, :write_when_maxdim_exceeds, nothing
    )

    # eigsolve kwargs
    eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
    eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
    eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)

    # SVD kwargseigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
    maxdim::Int = get(kwargs, :maxdim,100)
    cutoff::Float64 = get(kwargs, :cutoff,0.0)

    # TODO: add support for non-Hermitian DMRG
    ishermitian::Bool = get(kwargs, :ishermitian, true)

    # TODO: use this as preferred syntax for passing arguments
    # to eigsolve
    #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
    #                         verbosity = 0, ishermitian = true,
    #                         which_eigenvalue = :SR)
    #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)

    # Keyword argument deprecations
    if haskey(kwargs, :maxiter)
      error("""maxiter keyword has been replaced by eigsolve_krylovdim.
               Note: compared to the C++ version of ITensor,
               setting eigsolve_krylovdim 3 is the same as setting
               a maxiter of 2.""")
    end

    if haskey(kwargs, :errgoal)
      error("errgoal keyword has been replaced by eigsolve_tol.")
    end

    if haskey(kwargs, :quiet)
      error("quiet keyword has been replaced by outputlevel")
    end

    N = length(psi)

    if !isortho(psi) || orthocenter(psi) != 1
      orthogonalize!(psi, 1)
    end
    @assert isortho(psi) && orthocenter(psi) == 1

    position!(PH, psi, 1)
    energy = 0.0
    if order==2
        sweeps=1
    end
    for sw in 1:sweeps
      sw_time = @elapsed begin
        maxtruncerr = 0.0
        #=
        if !isnothing(write_when_maxdim_exceeds) &&
          maxdim(sweeps, sw) > write_when_maxdim_exceeds
          if outputlevel >= 2
            println(
              "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sw)), writing environment tensors to disk",
            )
          end
          PH = disk(PH)
        end
        =#
        print("forward\n")
        PH,psi=forward!(PH,psi,dt/2.,kwargs)

        print("backward\n")
        PH,psi=backward!(PH,psi,dt/2.,kwargs)
        #=
        for (b, ha) in sweepnext(N)
          @debug_check begin
            checkflux(psi)
            checkflux(PH)
          end

          if ha == 1
            ##Taken care of in two steps below, assuming that before the sweep position was set
            #@timeit_debug timer "dmrg: position!" begin
            #  position!(PH, psi, b)
            #end

            @debug_check begin
                checkflux(psi)
                checkflux(PH)
            end

            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end
            #@show inds(phi)

            #@show phi
            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,

                #maxiter=eigsolve_maxiter,
                )
            end
            #print(size(PH),size(vecs))
            #energy::Number = vals[1]
            #print(vecs)
            phi = vecs
            #@show inds(phi);
            ortho = ha == 1 ? "left" : "right"

            #drho = nothing
            #if noise(sweeps, sw) > 0.0
            #  @timeit_debug timer "dmrg: noiseterm" begin
            #    # Use noise term when determining new MPS basis
            #    drho = noise(sweeps, sw) * noiseterm(PH, phi, ortho)
            #  end
            #end

            @debug_check begin
                checkflux(phi)
            end

            #for j in 1:N
            #    @show inds(psi[j])
            #end
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            #for j in 1:N
            #    @show inds(psi[j])
            #end
            maxtruncerr = max(maxtruncerr, spec.truncerr)
            #@show inds(psi[b+1])

            if b+1<N
                position_leftonly!(PH, psi, b+1)
                phi=psi[b + 1]
                #@show inds(psi[b+1])
                #@show inds(phi)
                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1]=vecs        ##FIXME: necessary? was the prior assignment occurence of this assignment to new var?
                position_rightonly!(PH, psi, b+1)
            else
                position!(PH,psi,b)
            end

            @debug_check begin
                checkflux(psi)
                checkflux(PH)
            end

          elseif ha == 2
            #return (energy, psi)
            ortho = ha == 1 ? "left" : "right"
            if b+1<N
                phi=psi[b + 1]
                #@show inds(phi)
                #println(b,' ', ha)
                #println(indices(phi))
                #println(indices(PH))

                @timeit_debug timer "tdvp2: expintegrator-1site-backward" begin
                    vecs, propinfo = exponentiate(
                    PH,
                    -dt,
                    phi,
                    ishermitian=ishermitian,
                    tol=eigsolve_tol,
                    krylovdim=eigsolve_krylovdim,
                    #verbosity=1
                    #maxiter=eigsolve_maxiter,
                    )
                end
                #@show inds(vecs)
                psi[b+1] = vecs
                position_leftonly!(PH, psi, b)
            end
            #position_leftonly!(PH, psi, b)
            #position_rightonly!(PH, psi, b)

            @timeit_debug timer "tdvp2: psi[b]*psi[b+1]" begin
                phi = psi[b] * psi[b + 1]
            end

            @timeit_debug timer "tdvp2: expintegrator-2site" begin
                vecs, propinfo = exponentiate(
                PH,
                dt,
                phi,
                ishermitian=ishermitian,
                tol=eigsolve_tol,
                krylovdim=eigsolve_krylovdim,
                #verbosity=1
                #maxiter=eigsolve_maxiter,
                )
            end
            #energy::Number = vals[1]
            phi = vecs
            @timeit_debug timer "dmrg: replacebond!" begin
                spec = replacebond!(
                psi,
                b,
                phi;
                maxdim=maxdim,
                cutoff=cutoff,
                #eigen_perturbation=drho,
                ortho=ortho,
                normalize=true,
                which_decomp=which_decomp,
                svd_alg=svd_alg,
                )
            end
            maxtruncerr = max(maxtruncerr, spec.truncerr)

            if b > 1    #after completed backward-forward sweep we want to be left with 2-site PH
                position_rightonly!(PH, psi, b-1)
                #position_leftonly!(PH, psi, b-1)

            end

            end
          if outputlevel >= 2
            @printf(
              "Sweep %d, half %d, bond (%d,%d) energy=%.12f\n", sw, ha, b, b + 1, energy
            )
            @printf(
              "  Truncated using cutoff=%.1E maxdim=%d mindim=%d\n",
              cutoff(sweeps, sw),
              maxdim(sweeps, sw),
              mindim(sweeps, sw)
            )
            @printf(
              "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b))
            )
            flush(stdout)
          end

          sweep_is_done = (b == 1 && ha == 2)
          =#
          #=
          measure!(
            obs;
            energy=energy,
            psi=psi,
            bond=b,
            sweep=sw,
            half_sweep=ha,
            spec=spec,
            outputlevel=outputlevel,
            sweep_is_done=sweep_is_done,
          )
          =#
        end
      end
      if outputlevel >= 1
        #=
        @printf(
          "After sweep %d energy=%.12f maxlinkdim=%d maxerr=%.2E time=%.3f\n",
          sw,
          energy,
          maxlinkdim(psi),
          maxtruncerr,
          sw_time
        )
        =#
        flush(stdout)
      end

    return (energy, psi)
  end


text:
include("lidar_sensors.jl")

type MultiFeatureExtractor <: AbstractFeatureExtractor
    extract_core::Bool # [8] if true, extract lane offset, lane relative heading, speed, vehicle width, vehicle length, marker dist left+right, lane curvature
    extract_temporal::Bool # [6]  if true, extract accel, jerk, turnrate G, angular rate G, turnrate F, angular rate F
    extract_well_behaved::Bool # [3] if true, extract is_colliding, is_offroad, is_reversing
    extract_neighbor_features::Bool # [28] if true, extract everything that can be missing in tim2d features (ie, speed/dist fore, etc.)
    extract_carlidar_rangerate::Bool # [nbeams] if true, extract range rate as well as range

    carlidar::LidarSensor # [2*nbeams] NOTE: range rate can be turned off.
    roadlidar::RoadlineLidarSensor # [nbeams × depth]
    road_lidar_culling::RoadwayLidarCulling
end
function MultiFeatureExtractor(
    extract_core::Bool,
    extract_temporal::Bool,
    extract_well_behaved::Bool,
    extract_neighbor_features::Bool,
    extract_carlidar_rangerate::Bool,
    carlidar_nbeams::Int,
    roadlidar_nbeams::Int,
    roadlidar_nlanes::Int,
    ;
    carlidar_max_range::Float64 = 100.0, # [m]
    roadlidar_max_range::Float64 = 100.0, # [m]
    )

    carlidar = LidarSensor(carlidar_nbeams, max_range=carlidar_max_range, angle_offset=-π)
    roadlidar = RoadlineLidarSensor(roadlidar_nbeams, max_range=roadlidar_max_range, angle_offset=-π, max_depth=roadlidar_nlanes)

    MultiFeatureExtractor(
        extract_core,
        extract_temporal,
        extract_well_behaved,
        extract_neighbor_features,
        extract_carlidar_rangerate,
        carlidar,
        roadlidar,
        RoadwayLidarCulling(),
        )
end

AutomotiveDrivingModels.rec_length(ext::MultiFeatureExtractor) = ext.extract_temporal ? 3 : 1
function Base.length(ext::MultiFeatureExtractor)
    nbeams_carlidar = nbeams(ext.carlidar)
    nbeams_roadlidar = nbeams(ext.roadlidar)
    nlanes_roadlidar = nlanes(ext.roadlidar)

    8 * ext.extract_core +
    6 * ext.extract_temporal +
    3 * ext.extract_well_behaved +
    28 * ext.extract_neighbor_features +
    nbeams_carlidar * (1 + ext.extract_carlidar_rangerate) +
    nbeams_roadlidar * nlanes_roadlidar
end
function get_start_of_feature_section(ext::MultiFeatureExtractor)

    nbeams_carlidar = nbeams(ext.carlidar)
    nbeams_roadlidar = nbeams(ext.roadlidar)
    nlanes_roadlidar = nlanes(ext.roadlidar)

    i = 0
    core = temporal = well_behaved = neighbor = carlidar_range = carlidar_range_rate = roadlidar_range = -1

    if ext.extract_core
        core = 1
        i += 8
    end
    if ext.extract_temporal
        temporal = i
        i += 6
    end
    if ext.extract_well_behaved
        well_behaved = i
        i += 3
    end
    if ext.extract_neighbor_features
        neighbor = i
        i += 28
    end
    if nbeams_carlidar > 0
        carlidar_range = i
        i += nbeams_carlidar

        if ext.extract_carlidar_rangerate
            carlidar_range_rate = i
            i += nbeams_carlidar
        end
    end
    if nbeams_roadlidar > 0
        roadlidar_range = i
        i += nbeams_roadlidar * nlanes_roadlidar
    end

    (core, temporal, well_behaved, neighbor, carlidar_range, carlidar_range_rate, roadlidar_range, i)
end
function AutomotiveDrivingModels.pull_features!{F<:AbstractFloat}(ext::MultiFeatureExtractor,
    features::Vector{F},
    rec::SceneRecord,
    roadway::Roadway,
    vehicle_index::Int,
    pastframe::Int=0,
    )

    feature_index = 0
    scene = get_scene(rec, pastframe)
    veh_ego = scene[vehicle_index]

    d_ml = convert(Float64, get(MARKERDIST_LEFT, rec, roadway, vehicle_index, pastframe))
    d_mr = convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, vehicle_index, pastframe))

    if ext.extract_core
        #=
        1 - lane offset (positive to left) [m]
        2 - lane relative heading (positive to left) [rad]
        3 - speed [m/s]
        4 - vehicle length [m]
        5 - vehicle width [m]
        6 - lane curvature [1/m]
        7 - Marker Dist left
        8 - Marker Dist right
        =#

        features[feature_index+=1] = veh_ego.state.posF.t
        features[feature_index+=1] = veh_ego.state.posF.ϕ
        features[feature_index+=1] = veh_ego.state.v
        features[feature_index+=1] = veh_ego.def.length
        features[feature_index+=1] = veh_ego.def.width
        features[feature_index+=1] = convert(Float64, get(LANECURVATURE, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = d_ml
        features[feature_index+=1] = d_mr
    end

    if ext.extract_temporal
        #=
        1 - previous acceleration [m/s²]
        2 - previous jerk [m/s³]
        3 - previous turnrate G [rad/s]
        4 - previous angular rate G [rad/s²]
        5 - previous turnrate F [rad/s]
        6 - previous angular rate F [rad/s²]
        =#

        features[feature_index+=1] = convert(Float64, get(ACC, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, get(JERK, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, get(TURNRATEG, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, get(ANGULARRATEG, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, get(TURNRATEF, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, get(ANGULARRATEF, rec, roadway, vehicle_index, pastframe))
    end

    if ext.extract_well_behaved
        #=
        1 - is colliding
        2 - is offroad
        3 - is reversing
        =#

        features[feature_index+=1] = convert(Float64, get(IS_COLLIDING, rec, roadway, vehicle_index, pastframe))
        features[feature_index+=1] = convert(Float64, d_ml < -1.0 || d_mr < -1.0)
        features[feature_index+=1] = convert(Float64, veh_ego.state.v < 0.0)
    end

    if ext.extract_neighbor_features
        #=
        can be missing: - first is feature, 2nd is indicator whether it is missing (0 if exists, 1 if missing)
        (1, 2) - lane offset from left lane
        (3, 4) - lane offset from right lane
        (5, 6) - speed fore middle
        (7, 8) - dist fore middle (distance along the lane from the front of the ego car to the rear of the ego car in the same lane)
        (9, 10) - speed fore left
        (11,12) - dist fore left
        (13,14) - speed fore right
        (15,16) - dist fore right
        (17,18) - speed rear middle
        (19,20) - dist fore middle
        (21,22) - speed rear left
        (23,24) - dist rear left
        (25,26) - speed rear right
        (27,28) - dist rear right
        =#

        vtpf = VehicleTargetPointFront()
        vtpr = VehicleTargetPointRear()
        fore_M = get_neighbor_fore_along_lane(      scene, vehicle_index, roadway, vtpf, vtpr, vtpf)
        fore_L = get_neighbor_fore_along_left_lane( scene, vehicle_index, roadway, vtpf, vtpr, vtpf)
        fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, vtpf, vtpr, vtpf)
        rear_M = get_neighbor_rear_along_lane(      scene, vehicle_index, roadway, vtpr, vtpf, vtpr)
        rear_L = get_neighbor_rear_along_left_lane( scene, vehicle_index, roadway, vtpr, vtpf, vtpr)
        rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, vtpr, vtpf, vtpr)

        function set_feature_missing!(features::Vector, i::Int)
            features[i] = 0.0
            features[i+1] = 1.0
            features
        end
        function set_feature!(features::Vector, i::Int, v::Float64)
            features[i] = v
            features[i+1] = 0.0
            features
        end
        function set_dual_feature!(features::Vector, i::Int, f::FeatureValue)
            if f.i == FeatureState.MISSING
                set_feature_missing!(features, i)
            else
                set_feature!(features, i, f.v)
            end
            features
        end
        function set_speed_and_distance!(features::Vector, i::Int, neigh::NeighborLongitudinalResult)
            neigh.ind != 0 ? set_feature!(features, i, scene[neigh.ind].state.v) :
                              set_feature_missing!(features, i)
            neigh.ind != 0 ? set_feature!(features, i+2, neigh.Δs) :
                              set_feature_missing!(features, i+2)
            features
        end

        set_dual_feature!(features, feature_index+=1, get(LANEOFFSETLEFT, rec, roadway, vehicle_index, pastframe))
        feature_index+=1
        set_dual_feature!(features, feature_index+=1, get(LANEOFFSETRIGHT, rec, roadway, vehicle_index, pastframe))
        feature_index+=1

        set_speed_and_distance!(features, feature_index+=1, fore_M); feature_index+=3
        set_speed_and_distance!(features, feature_index+=1, fore_L); feature_index+=3
        set_speed_and_distance!(features, feature_index+=1, fore_R); feature_index+=3
        set_speed_and_distance!(features, feature_index+=1, rear_M); feature_index+=3
        set_speed_and_distance!(features, feature_index+=1, rear_L); feature_index+=3
        set_speed_and_distance!(features, feature_index+=1, rear_R); feature_index+=3
    end

    nbeams_carlidar = nbeams(ext.carlidar)
    if nbeams_carlidar > 0
        observe!(ext.carlidar, scene, roadway, vehicle_index)
        copy!(features, feature_index+=1, ext.carlidar.ranges, 1)
        feature_index += nbeams_carlidar - 1
        if ext.extract_carlidar_rangerate
            copy!(features, feature_index+=1, ext.carlidar.range_rates, 1)
            feature_index += nbeams_carlidar - 1
        end
    end

    nbeams_roadlidar = nbeams(ext.roadlidar)
    if nbeams_roadlidar > 0
        if ext.road_lidar_culling.is_leaf
            observe!(ext.roadlidar, scene, roadway, vehicle_index)
        else
            observe!(ext.roadlidar, scene, roadway, vehicle_index, ext.road_lidar_culling)
        end
        copy!(features, feature_index+=1, ext.roadlidar.ranges)
        feature_index += length(ext.roadlidar.ranges) - 1
    end

    @assert(feature_index == length(features))
    if findfirst(v->isnan(v), features) != 0
        error("feature $(findfirst(v->isnan(v), features)) is nan")
    end

    features
end

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
#=##############################################################################################
Filename: system.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this file defines freestream and interstream properties
=###############################################################################################

"""
    Freestream{F, VF <: AbstractArray}

Describes the freestream.

# Arguments:

* `Vinf::Float64` : magnitude of the freestream
* `alpha::Float64` : angle of attack of the aircraft
* `beta::Float64` : sideslip angle of the aircraft
* `Omega::Vector{Float64}` : angular velocity about roll, pitch, and yaw axes

"""
struct Freestream{TF1,TF2,TF3,TF4}
    vinf::TF1
    alpha::TF2
    beta::TF3
    Omega::StaticArrays.SArray{Tuple{3},TF4,1,3}
end

"""
    freestream2vector(freestream)

# Argument:

* `freestream::Freestream`: freestream properties

# Returns:

* `Vinf::Vector{Float64}`: freestream velocity in the body frame (aft-starboard-up)

"""
function freestream2vector(freestream)

    vinf = LA.norm(freestream.vinf)
    alpha = freestream.alpha
    beta = freestream.beta
    Vinf = vinf .*
        [   cos(alpha) * cos(beta),
            -sin(beta),
            sin(alpha) * cos(beta)
        ]

    return Vinf
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
module GeneralizedMonteCarlo

using Pkg, Pkg.TOML, Distributions, DocStringExtensions, Documenter, Symbolics
"""
$(SIGNATURES)

Either makes a random draw of a float from the distribution based on its specified parameters/pmf or returns the input as a float

#Example

```julia-repl
julia> issubset(GeneralizedMonteCarlo.inputPassthrough([.01 .02 .03 .04 .05 .08; .1 .2 .3 .25 .10 .05]), [.01 .02 .03 .04 .05 .08])
true
```

```julia-repl
julia> 0 ≤ GeneralizedMonteCarlo.inputPassthrough(Uniform(0, 1)) ≤ 1
true
```

```julia-repl
julia> GeneralizedMonteCarlo.inputPassthrough(3) == 3.0
true
```

"""
function inputPassthrough(x :: Union{Real, Sampleable, Matrix{Float64}})
    if typeof(x) <: Matrix
        return rand(Distributions.DiscreteNonParametric(x[1, :], vec(x[2:size(x)[1], :])))*1.0 #draws randomly from given pmf
    elseif typeof(x) <: Sampleable
        return rand(x)*1.0
    elseif typeof(x) <: Real
        return x*1.0
    end
end

"""
$(SIGNATURES)

Creates a new Monte Carlo simulation of calculating revenue according to a model of the users preferences, when given the number of simulations to calculate, as well as the distributions (can be user defined with a pmf or from https://github.com/JuliaStats/Distributions.jl). Or the user can just pass through a Real number. Returns array of random draws and the values from the equation.

Notice that:
julia> @variables var1 var2 ...
needs to be ran before the function, as shown in the example below.

#Example

```julia-repl
julia> @variables iws ms p; GeneralizedMonteCarlo.GMC(10_000, [iws, ms, 190 + ms + p], [iws, ms, p], [[10_000_000 10_500_000 11_000_000; 0.2 0.6 0.2], [0.01 0.02 0.03 0.04 0.05 0.08; 0.10 0.20 0.30 0.25 0.10 0.05], [-3.0, 3.0; 0.5 0.5]])
```

"""
function GMC(NSimulations :: Int, equations :: Array, symbols :: Array, probs :: Array)
    draws = zeros(Float64, NSimulations, length(equations))
    values = zeros(Num, NSimulations, length(equations))
    floatvalues = zeros(Float64, NSimulations, length(equations))
    for n in 1:NSimulations
        for l in 1:length(equations)
            draws[n, l] = inputPassthrough(probs[l])
        end
        for l in 1:length(equations)
            values[n,l] = substitute(Num(equations[l]), Dict(zip(symbols, draws[n, :])))[1]
            floatvalues[n, l] = parse(Float64, string(values[n, l]))
        end
    end
    return floatvalues, draws
end
end


text:
# this script does not run as a whole
# (the code snippets are for demo-purposes only,
# and should be used on their own, and/or in the REPL
# on Julia prompt:
# julia> stdin
# Base.TTY(Base.Libc.WindowsRawSocket(0x000000000000027c) open, 0 bytes waiting)
# julia> stdout
# Base.TTY(Base.Libc.WindowsRawSocket(0x0000000000000280) open, 0 bytes waiting)

# read(stdin, Char)
write(stdout, "Julia") #> Julia5
# read(stdin,3)
# readline(stdin)
# input:  # Julia
# output: # "Julia\r\n" (Windows) - "Julia\n" (Linux)

# stream can be stdin or any other input stream:
stream = stdin
for line in eachline(stream)
    println("Found $line")
end

while !eof(stream)
    x = read(stream, Char)
    println("Found: $x")
end

function process(item)
    # do the processing
end

# files:
fieldnames(IOStream)
#> (:handle, :ios, :name, :mark)

IOStream.types
#> svec(Ptr{Nothing}, Array{UInt8,1}, AbstractString, Int64)

# reading:
fname = "example.dat"
f1 = open(fname)
# IOStream(<file example.dat>)

data = readlines(f1)
# 3-element Array{String,1}:
#  "this is line 1."
#  "42; 3.14"
#  "this is line 3."

for line in data
	println(line)
	process(line)
end
close(f1)

open(fname) do file
    process(file)
end

open(fname) do file
    for line in eachline(file)
        print(line)
        process(line)
    end
end

# writing:
fname = "example2.dat"
f2 = open(fname, "w")
write(f2, "I write myself to a file\n") #> 24
println(f2, "even with println!")
close(f2)

# looping over all files in a folder:
for file in readdir()
    process(file)
end

# example:
mydir = pwd();
cd("..")

for fn in readdir()
    println(fn)
end
# possible output:
#  .DS_Store
# .git
# Chapter 1
# Chapter 2
# Chapter 3
# Chapter 4
# Chapter 5
# Chapter 6
# Chapter 7
# Chapter 8
# Chapter 9
# LICENSE
# README.md

cd(mydir)


text:
one_degree_interface_heights() =
    [-5244.5
     -4834.0
     -4446.5
     -4082.0
     -3740.5
     -3422.0
     -3126.5
     -2854.0
     -2604.5
     -2378.0
     -2174.4
     -1993.6
     -1834.7
     -1695.6
     -1572.8
     -1461.4
     -1356.9
     -1255.5
     -1155.5
     -1056.3
      -958.0
      -861.4
      -767.5
      -677.3
      -592.2
      -513.3
      -441.7
      -378.2
      -323.2
      -276.7
      -238.3
      -207.2
      -182.3
      -162.5
      -146.5
      -133.0
      -121.3
      -110.5
      -100.1
       -90.0
       -80.0
       -70.0
       -60.0
       -50.0
       -40.0
       -30.0
       -20.0
       -10.0
         0.0]