text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra

import ClimaCore: slab, Domains, Meshes, Topologies, Spaces
import ClimaCore.Geometry: Geometry

@testset "1d domain space" begin
    FT = Float64
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-3)..Geometry.XPoint{FT}(5),
        Geometry.YPoint{FT}(0)..Geometry.YPoint{FT}(0),
        x1periodic = true,
        x2periodic = false,
        x2boundary = (:south, :north),
    )
    mesh = Meshes.EquispacedRectangleMesh(domain, 1, 1)
    grid_topology = Topologies.GridTopology(mesh)

    quad = Spaces.Quadratures.GLL{4}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace1D(grid_topology, quad)

    coord_data = Spaces.coordinates_data(space)
    @test eltype(coord_data) == Geometry.XPoint{Float64}

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (4, 1, 1)
    coord_slab = slab(Spaces.coordinates_data(space), 1)
    @test coord_slab[1] == Geometry.XPoint{FT}(-3)
    @test coord_slab[4] == Geometry.XPoint{FT}(5)

    local_geometry_slab = slab(space.local_geometry, 1)
    dss_weights_slab = slab(space.dss_weights, 1)

    for i in 1:4
        @test Geometry.components(local_geometry_slab[i].∂x∂ξ) ≈
              @SMatrix [8 / 2]
        @test Geometry.components(local_geometry_slab[i].∂ξ∂x) ≈
              @SMatrix [2 / 8]
        @test local_geometry_slab[i].J ≈ (8 / 2)
        @test local_geometry_slab[i].WJ ≈ (8 / 2) * weights[i]
        if i in (1, 4)
            @test dss_weights_slab[i] ≈ 1 / 2
        else
            @test dss_weights_slab[i] ≈ 1
        end
    end

end

@testset "1×1 domain space" begin
    FT = Float32
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-3)..Geometry.XPoint{FT}(5),
        Geometry.YPoint{FT}(-2)..Geometry.YPoint{FT}(8),
        x1periodic = true,
        x2periodic = false,
        x2boundary = (:south, :north),
    )
    mesh = Meshes.EquispacedRectangleMesh(domain, 1, 1)
    grid_topology = Topologies.GridTopology(mesh)

    quad = Spaces.Quadratures.GLL{4}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (4, 4, 2, 1)
    coord_slab = slab(Spaces.coordinates_data(space), 1)
    @test coord_slab[1, 1] ≈ Geometry.XYPoint{FT}(-3.0, -2.0)
    @test coord_slab[4, 1] ≈ Geometry.XYPoint{FT}(5.0, -2.0)
    @test coord_slab[1, 4] ≈ Geometry.XYPoint{FT}(-3.0, 8.0)
    @test coord_slab[4, 4] ≈ Geometry.XYPoint{FT}(5.0, 8.0)

    local_geometry_slab = slab(space.local_geometry, 1)
    dss_weights_slab = slab(space.dss_weights, 1)

    for i in 1:4, j in 1:4
        @test Geometry.components(local_geometry_slab[i, j].∂x∂ξ) ≈
              @SMatrix [8/2 0; 0 10/2]
        @test Geometry.components(local_geometry_slab[i, j].∂ξ∂x) ≈
              @SMatrix [2/8 0; 0 2/10]
        @test local_geometry_slab[i, j].J ≈ (10 / 2) * (8 / 2)
        @test local_geometry_slab[i, j].WJ ≈
              (10 / 2) * (8 / 2) * weights[i] * weights[j]
        if i in (1, 4)
            @test dss_weights_slab[i, j] ≈ 1 / 2
        else
            @test dss_weights_slab[i, j] ≈ 1
        end
    end

    @test length(space.boundary_surface_geometries) == 2
    @test keys(space.boundary_surface_geometries) == (:south, :north)
    @test sum(parent(space.boundary_surface_geometries.north.sWJ)) ≈ 8
    @test parent(space.boundary_surface_geometries.north.normal)[1, :, 1] ≈
          [0.0, 1.0]
end


text:
"""
    classnames() -> Vector{String}

Return the 10 names for the Fashion-MNIST classes as a vector of
strings.
"""
classnames() = CLASSES

"""
    traintensor([T = N0f8], [indices]; [dir]) -> Array{T}

Same as [`MNIST.traintensor`](@ref) but for the `FashionMNIST` dataset.

$(download_docstring("FashionMNIST", DEPNAME))
"""
function traintensor(::Type{T}, args...; dir = nothing) where T
    path = datafile(DEPNAME, TRAINIMAGES, dir)
    images = Reader.readimages(path, args...)
    bytes_to_type(T, images)
end

function traintensor(args...; dir = nothing)
    traintensor(N0f8, args...; dir = dir)
end

"""
    testtensor([T = N0f8], [indices]; [dir]) -> Array{T}

Same as [`MNIST.testtensor`](@ref) but for the `FashionMNIST` dataset.
```

$(download_docstring("FashionMNIST", DEPNAME))
"""
function testtensor(::Type{T}, args...; dir = nothing) where T
    path = datafile(DEPNAME, TESTIMAGES, dir)
    images = Reader.readimages(path, args...)
    bytes_to_type(T, images)
end

function testtensor(args...; dir = nothing)
    testtensor(N0f8, args...; dir = dir)
end

"""
    trainlabels([indices]; [dir])

Returns the Fashion-MNIST **trainset** labels corresponding to
the given `indices` as an `Int` or `Vector{Int}`. The values of
the labels denote the zero-based class-index that they represent
(see [`FashionMNIST.classnames`](@ref) for the corresponding
names). If `indices` is omitted, all labels are returned.

```julia-repl
julia> FashionMNIST.trainlabels() # full training set
60000-element Array{Int64,1}:
 9
 0
 ⋮
 0
 5

julia> FashionMNIST.trainlabels(1:3) # first three labels
3-element Array{Int64,1}:
 9
 0
 0

julia> y = FashionMNIST.trainlabels(1) # first label
9

julia> FashionMNIST.classnames()[y + 1] # corresponding name
"Ankle boot"
```

$(download_docstring("FashionMNIST", DEPNAME))
"""
function trainlabels(args...; dir = nothing)
    path = datafile(DEPNAME, TRAINLABELS, dir)
    Vector{Int}(Reader.readlabels(path, args...))
end

function trainlabels(index::Integer; dir = nothing)
    path = datafile(DEPNAME, TRAINLABELS, dir)
    Int(Reader.readlabels(path, index))
end

"""
    testlabels([indices]; [dir])

Returns the Fashion-MNIST **testset** labels corresponding to the
given `indices` as an `Int` or `Vector{Int}`. The values of the
labels denote the class-index that they represent (see
[`FashionMNIST.classnames`](@ref) for the corresponding names).
If `indices` is omitted, all labels are returned.

```julia-repl
julia> FashionMNIST.testlabels() # full test set
10000-element Array{Int64,1}:
 9
 2
 ⋮
 1
 5

julia> FashionMNIST.testlabels(1:3) # first three labels
3-element Array{Int64,1}:
 9
 2
 1

julia> y = FashionMNIST.testlabels(1) # first label
9

julia> FashionMNIST.classnames()[y + 1] # corresponding name
"Ankle boot"
```

$(download_docstring("FashionMNIST", DEPNAME))
"""
function testlabels(args...; dir = nothing)
    path = datafile(DEPNAME, TESTLABELS, dir)
    Vector{Int}(Reader.readlabels(path, args...))
end

function testlabels(index::Integer; dir = nothing)
    path = datafile(DEPNAME, TESTLABELS, dir)
    Int(Reader.readlabels(path, index))
end

"""
    traindata([T = N0f8], [indices]; [dir]) -> images, labels

Same as [`MNIST.traindata`](@ref) but for the `FashionMNIST` dataset.

$(download_docstring("FashionMNIST", DEPNAME))

Take a look at [`FashionMNIST.traintensor`](@ref) and
[`FashionMNIST.trainlabels`](@ref) for more information.
"""
function traindata(::Type{T}, args...; dir = nothing) where T
    (traintensor(T, args...; dir = dir),
     trainlabels(args...; dir = dir))
end

traindata(args...; dir = nothing) = traindata(N0f8, args...; dir = dir)

"""
    testdata([T = N0f8], [indices]; [dir]) -> images, labels

Same as [`MNIST.testdata`](@ref) but for the `FashionMNIST` dataset.

$(download_docstring("FashionMNIST", DEPNAME))

Take a look at [`FashionMNIST.testtensor`](@ref) and
[`FashionMNIST.testlabels`](@ref) for more information.
"""
function testdata(::Type{T}, args...; dir = nothing) where T
    (testtensor(T, args...; dir = dir),
     testlabels(args...; dir = dir))
end

testdata(args...; dir = nothing) = testdata(N0f8, args...; dir = dir)

ext: jl


text:
# apply is operations over substitution

#include("primitives.jl")

"""
apply substitution to a term, substitution
"""
function apply(vars::Vlist, sym::Number, subst::Tlist)
 return sym
end

function apply(vars::Vlist, sym::Symbol, subst::Tlist)
 for i in 1:length(vars)
  if sym == vars[i]; return subst[i] end
 end
 return sym
end

"""
apply: ΣxΣ → Σ
"""
function apply(vars::Vlist, subst1::Tlist, subst2::Tlist)
 nterm = []
 for arg in subst1
  narg = apply(vars, arg, subst2)
  if isempty(nterm)
   nterm = [narg]
  elseif ndims(subst1) == 1
     nterm = vcat(nterm, [narg])
  else
     nterm = hcat(nterm, [narg])
  end
 end
 nterm
end

function apply(vars::Vlist, term::Expr, subst::Tlist)
 nterm = deepcopy(term)
 for i in 1:length(term.args)
  arg = term.args[i]
  narg = apply(vars, arg, subst)
  nterm.args[i] = narg
 end
 nterm
end

"""
apply sigma2 to sigma1
simply make v<-t in σ1 to v<-t*σ2
"""
function apply2(vars::Vlist, σ1::Tlist, σ2::Tlist)
  map(s->apply(vars, s, σ2), σ1)
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using ExperienceAnalysis
using Test
using Dates

@testset "Anniversary" begin
    issue = Date(2016, 7, 4)
    termination = Date(2020, 1, 17)

    @testset "Year" begin
        exp = exposure(ExperienceAnalysis.Anniversary(Year(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 7, 4))
        @test exp[end] == (from = Date(2019, 7, 4), to = Date(2020, 1, 17))
    end

    @testset "Year with full exp" begin
    exp = exposure(ExperienceAnalysis.Anniversary(Year(1)), issue, termination,true)

    @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 7, 4))
    @test exp[end] == (from = Date(2019, 7, 4), to = Date(2020, 7, 4))
end

    @testset "Month" begin
        exp = exposure(ExperienceAnalysis.Anniversary(Month(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2016, 8, 4))
        @test exp[end] == (from = Date(2020, 1, 4), to = Date(2020, 1, 17))
    end
end

@testset "Calendar" begin
    issue = Date(2016, 7, 4)
    termination = Date(2020, 1, 17)

    @testset "Year" begin
        exp = exposure(ExperienceAnalysis.Calendar(Year(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 1, 1))
        @test exp[2]   == (from = Date(2017, 1, 1), to = Date(2018, 1, 1))
        @test exp[end] == (from = Date(2020, 1, 1), to = Date(2020, 1, 17))
    end

    @testset "Year with full exp" begin
    exp = exposure(ExperienceAnalysis.Calendar(Year(1)), issue, termination,true)

    @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 1, 1))
    @test exp[2]   == (from = Date(2017, 1, 1), to = Date(2018, 1, 1))
    @test exp[end] == (from = Date(2020, 1, 1), to = Date(2021, 1, 1))
end

    @testset "Month" begin
        exp = exposure(ExperienceAnalysis.Calendar(Month(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2016, 8, 1))
        @test exp[2]   == (from = Date(2016, 8, 1), to = Date(2016, 9, 1))
        @test exp[end] == (from = Date(2020, 1, 1), to = Date(2020, 1, 17))
    end
end

@testset "AnniversaryCalendar" begin
    issue = Date(2016, 7, 4)
    termination = Date(2020, 1, 17)

    @testset "Year/Year" begin
        exp = exposure(ExperienceAnalysis.AnniversaryCalendar(Year(1), Year(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 1, 1))
        @test exp[2]   == (from = Date(2017, 1, 1), to = Date(2017, 7, 4))
        @test exp[end] == (from = Date(2020, 1, 1), to = Date(2020, 1, 17))
    end

    @testset "Year/Year with full exp" begin
    exp = exposure(ExperienceAnalysis.AnniversaryCalendar(Year(1), Year(1)), issue, termination,true)

    @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2017, 1, 1))
    @test exp[2]   == (from = Date(2017, 1, 1), to = Date(2017, 7, 4))
    @test exp[end] == (from = Date(2020, 1, 1), to = Date(2020, 7, 4))
end

    @testset "Month/Year" begin
        exp = exposure(ExperienceAnalysis.AnniversaryCalendar(Month(1), Year(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2016, 8, 4))
        @test exp[6]   == (from = Date(2016, 12, 4), to = Date(2017, 1, 1))
        @test exp[7]   == (from = Date(2017, 1, 1), to = Date(2017, 1, 4))
        @test exp[end] == (from = Date(2020, 1, 4), to = Date(2020, 1, 17))
    end

    @testset "Month/Month" begin
        exp = exposure(ExperienceAnalysis.AnniversaryCalendar(Month(1), Month(1)), issue, termination)

        @test exp[1]   == (from = Date(2016, 7, 4), to = Date(2016, 8, 1))
        @test exp[2]   == (from = Date(2016, 8, 1), to = Date(2016, 8, 4))
        @test exp[end] == (from = Date(2020, 1, 4), to = Date(2020, 1, 17))
    end
end

@testset "Experience Study Calculations" begin
    # https://www.soa.org/globalassets/assets/files/research/experience-study-calculations.pdf

    @testset "Section 4.3.1" begin
        lives = [(
            start = Date(2010,5,10),
            termination = nothing,
            status = "Inforce",
        ),
        (
            start = Date(2010,9,27),
            termination = Date(2012,2,16),
            status = "Claim",
        ),
        (
            start = Date(2010,7,3),
            termination = Date(2012,10,21),
            status = "Lapse",
        ),
        (
            start = Date(2009,2,12),
            termination = nothing,
            status = "Inforce",
        ),
        (
            start = Date(2009,10,30),
            termination = Date(2013,12,27),
            status = "Claim",
        ),
        (
            start = Date(2009,7,5),
            termination = Date(2010,3,17),
            status = "Claim",
        ),
        ]

        study_start = Date(2010,1,1)
        study_end = Date(2014,1,1)

        exp = Dict(

            i => exposure(
                ExperienceAnalysis.Anniversary(Year(1)),
                l.start,
                isnothing(l.termination) ? study_end : min(study_end,l.termination),
                l.status == "Claim"
            )
            for (i,l) in enumerate(lives))

        days = [e.to - e.from for e in exp[1]]
        @test days == Day.([365,366,365,236])

        days = [e.to - e.from for e in exp[2]]
        @test days == Day.([365,366])

        days = [e.to - e.from for e in exp[3]]
        @test days == Day.([365,366,110])

        # these fail because the study start doesn't truncate the starting date
        # and the `from` argument needs to be the anniv for the right date iteration
        days = [e.to - e.from for e in exp[4]]
        @test_broken days == Day.([42,365,365,366,323])

        days = [e.to - e.from for e in exp[5]]
        @test_broken days == Day.([302,365,365,366,365])

        days = [e.to - e.from for e in exp[6]]
        @test_broken days == Day.([185])

    end

end

max_stars_repo_name: JuliaActuary/ExperienceAnalysis

max_issues_repo_name: JuliaActuary/ExperienceAnalysis

max_forks_repo_name: JuliaActuary/ExperienceAnalysis


text:
module GroundEffects

using .Meta: isexpr

@nospecialize

"""
    lower(ex::Expr) :: Expr

Desugar subexpressions of `ex` as long as they are representable by Julia's
surface syntax.  The lowered output can be used as an output of macro.
"""
lower(ex::Expr) = lower(defaulthandlers(), ex)
lower(handlers::AbstractVector, ex::Expr) = Dispatcher(handlers)(ex)

defaulthandlers() = Any[
    handle_macrocall,
    handle_vect,
    handle_vcat,
    handle_hcat,
    handle_typed_vcat,
    handle_typed_hcat,
    handle_ref,
    handle_assignment,
    handle_inplace_materialize,
    handle_dotcall,
    handle_dotupdate,
    handle_getproperty,
    handle_do,
    handle_recursion,
]

struct Dispatcher
    handlers::Vector{Any}
end

function (lower::Dispatcher)(ex)
    for h in lower.handlers
        y = h(lower, ex)
        y === nothing || return something(y)
    end
    return ex
end

handle_recursion(lower, ::Any) = nothing
handle_recursion(lower, ex::Expr) = Expr(ex.head, map(lower, ex.args)...)

handle_macrocall(_, ex) = isexpr(ex, :macrocall) ? ex : nothing

function handle_vect(lower, ex)
    isexpr(ex, :vect) || return nothing
    return Expr(:call, Base.vect, map(lower, ex.args)...)
end

function handle_vcat(lower, ex)
    isexpr(ex, :vcat) || return nothing
    if all(isexpr.(ex.args, :row))
        rows = Tuple(length(a.args) for a in ex.args)
        return :($(Base.hvcat)($rows, $((
            lower(r) for a in ex.args for r in a.args
        )...)))
    else
        return :($(Base.vcat)($(map(lower, ex.args)...)))
    end
end

handle_hcat(lower, ex) =
    isexpr(ex, :hcat) ? :($(Base.hcat)($(map(lower, ex.args)...))) : nothing

function handle_typed_vcat(lower, ex)
    isexpr(ex, :typed_vcat) || return nothing
    elements = ex.args[2:end]
    if all(isexpr.(elements, :row))
        rows = Tuple(length(a.args) for a in elements)
        return :($(Base.typed_hvcat)($(lower(ex.args[1])), $rows, $((
            lower(r) for a in elements for r in a.args
        )...)))
    else
        return :($(Base.typed_vcat)($(map(lower, ex.args)...)))
    end
end

function handle_typed_hcat(lower, ex)
    isexpr(ex, :typed_hcat) || return nothing
    return :($(Base.typed_hcat)($(map(lower, ex.args)...)))
end

function handle_ref(lower, ex)
    isexpr(ex, :ref) || return nothing
    statements, collection, indices = _handle_ref(lower, ex)
    push!(statements, Expr(:call, Base.getindex, collection, indices...))
    if length(statements) == 1
        return statements[1]
    else
        return Expr(:block, statements...)
    end
end

function _handle_ref(lower, ex)
    statements = []
    if length(ex.args) == 1
        collection = ex.args[1]
        indices = []
    else
        if ex.args[1] isa Symbol
            collection = ex.args[1]
        else
            @gensym collection
            push!(statements, :($collection = $(ex.args[1])))
        end
        indices = lower_indices(lower, collection, ex.args[2:end])
    end
    return statements, collection, indices
end

function handle_assignment(lower, ex)
    isexpr(ex, :(=)) || return nothing
    lhs = Any[ex.args[1]]
    rhs = ex.args[2]
    while isexpr(rhs, :(=))
        push!(lhs, rhs.args[1])
        rhs = rhs.args[2]
    end
    @gensym rhsvalue
    args = mapfoldl(append!, lhs, init=Any[:($rhsvalue = $(lower(rhs)))]) do l
        _handle_assignment(lower, l, rhsvalue)
    end
    return Expr(:block, args..., rhsvalue)
end

function _handle_assignment(lower, lhs, rhs)
    if isexpr(lhs, :ref)
        statements, collection, indices = _handle_ref(lower, lhs)
        push!(statements, Expr(:call, Base.setindex!, collection, rhs, indices...))
        return statements
    elseif isexpr(lhs, :.)
        @assert length(lhs.args) == 2
        return [Expr(
            :call,
            Base.setproperty!,
            lower(lhs.args[1]),
            lhs.args[2],
            rhs,
        )]
    end
    return [:($lhs = $rhs)]
end

lower_indices(lower, collection, indices) =
    map(index -> lower_index(lower, collection, index), indices)

function lower_index(lower, collection, index)
    ex = handle_dotcall(index) do ex
        lower_index(lower, collection, ex)
    end
    ex === nothing || return something(ex)

    if isexpr(index, :call)
        return Expr(:call, lower_indices(lower, collection, index.args)...)
    elseif index === :end
        return :($(Base.lastindex)($collection))
    end
    return lower(index)
end

function handle_inplace_materialize(lower, ex)
    isexpr(ex, :.=) || return nothing
    @assert length(ex.args) == 2
    a1, a2 = map(lower, ex.args)
    return :($(Base.materialize!)($a1, $(Base.broadcasted)(identity, $a2)))
end

function isdotopcall(ex)
    ex isa Expr && !isempty(ex.args) || return false
    op = ex.args[1]
    return op isa Symbol && Base.isoperator(op) && startswith(string(op), ".")
end

isdotcall(ex) =
    isexpr(ex, :.) && length(ex.args) == 2 && isexpr(ex.args[2], :tuple)

function handle_dotcall(lower, ex)
    isdotcall(ex) || isdotopcall(ex) || return nothing
    return Expr(:call, Base.materialize, handle_lazy_dotcall(lower, ex))
end

function handle_lazy_dotcall(lower, ex)
    if isdotcall(ex)
        args = [
            lower(ex.args[1])
            map(x -> handle_lazy_dotcall(lower, x), ex.args[2].args)
        ]
    elseif isdotopcall(ex)
        args = [
            Symbol(String(ex.args[1])[2:end])
            map(x -> handle_lazy_dotcall(lower, x), ex.args[2:end])
        ]
    else
        return lower(ex)
    end
    return Expr(:call, Base.broadcasted, args...)
end

function handle_dotupdate(lower, ex)
    ex isa Expr || return nothing
    m = match(r"^(\.)?(([^.]+)=)$", string(ex.head))
    m === nothing && return nothing
    # e.g., `ex.head == :.+=`
    op = Symbol(m.captures[3])  # e.g., `:+`
    @assert length(ex.args) == 2
    a1, a2 = map(lower, ex.args)
    if m.captures[1] === nothing
        # e.g., `ex.head == :+=`
        if !(a1 isa Symbol)
            error("Assignment destination must be a symbol. Got:\n", a1)
        end
        return :($a1 = $op($a1, $a2))
    elseif a1 isa Symbol
        return :($(Base.materialize!)($a1, $(Base.broadcasted)($op, $a1, $a2)))
    else
        @gensym lhs
        return Expr(
            :block,
            :($lhs = $a1),
            :($(Base.materialize!)($lhs, $(Base.broadcasted)($op, $lhs, $a2))),
        )
    end
end

#=
function handle_comparison(lower, ex)
end
=#

function handle_getproperty(lower, ex)
    isexpr(ex, :.) || return nothing
    return Expr(:call, Base.getproperty, lower(ex.args[1]), ex.args[2])
end

function handle_do(lower, ex)
    isexpr(ex, :do) || return nothing
    call, lambda = map(lower, ex.args)
    return Expr(call.head, call.args[1], lambda, call.args[2:end]...)
end

end # module


text:
"""
    Smooth_Cut(x_mc::SMCg{N,T},x_mc_int::SMCg{N,T})

An operator that cuts the `x_mc` object using the `x_mc_int bounds` in a
differentiable fashion.
"""
function Smooth_Cut(x_mc::SMCg{N,T},x_mc_int::SMCg{N,T}) where {N,T<:AbstractFloat}
  t_cv::SMCg{N,T} = max(x_mc,x_mc_int.Intv.lo)
  t_cc::SMCg{N,T} = min(x_mc,x_mc_int.Intv.hi)
  return SMCg{N,T}(t_cc.cc,t_cv.cv,t_cc.cc_grad,t_cv.cv_grad,
                   (x_mc.Intv ∩ x_mc_int.Intv),(t_cv.cnst && t_cc.cnst),
                   x_mc.IntvBox,x_mc.xref)
end

"""
    Final_Cut(x_mc::SMCg{N,T},x_mc_int::SMCg{N,T})

An operator that cuts the `x_mc` object using the `x_mc_int bounds` in a
differentiable or nonsmooth fashion as specified by the `MC_param.mu flag`.
"""
function Final_Cut(x_mc::SMCg{N,T},x_mc_int::SMCg{N,T}) where {N,T<:AbstractFloat}
  if (MC_param.mu < 1)
    Intv::Interval = x_mc.Intv ∩ x_mc_int.Intv
    if (x_mc.cc <= x_mc_int.cc)
      cc::T = x_mc.cc
      cc_grad::SVector{N,T} = x_mc.cc_grad
    else
      cc = x_mc_int.cc
      cc_grad = x_mc_int.cc_grad
    end
    if (x_mc.cv >= x_mc_int.cv)
      cv::T = x_mc.cv
      cv_grad::SVector{N,T} = x_mc.cv_grad
    else
      cv = x_mc_int.cv
      cv_grad = x_mc_int.cv_grad
    end
    x_mc::SMCg{N,T} = SMCg{N,T}(cc,cv,cc_grad,cv_grad,(x_mc.Intv ∩ x_mc_int.Intv),x_mc.cnst,x_mc.IntvBox,x_mc.xref)
  else
    x_mc = Smooth_Cut(x_mc,x_mc_int)
  end
  return x_mc
end

"""
    Rnd_Out_Z_Intv(z_mct::SMCg{N,T},epsvi::Float64)

Rounds the interval of the `z_mct` vector elements out by `epsvi`.
"""
function Rnd_Out_Z_Intv(z_mct::Vector{SMCg{N,T}},epsvi::S) where {N,S<:AbstractFloat,T<:AbstractFloat}
  epsv::T = convert(T,epsvi)
  return [SMCg{N,T}(z_mct[i].cc,z_mct[i].cv,
             z_mct[i].cc_grad, z_mct[i].cv_grad,
             Interval{T}(z_mct[i].Intv.lo-epsv, z_mct[i].Intv.hi+epsv),
             z_mct[i].cnst, z_mct[i].IntvBox,z_mct[i].xref) for i=1:length(z_mct)]
end

"""
    Rnd_Out_Z_All(z_mct::Vector{SMCg{N,T}},epsvi::S)

Rounds the interval and relaxation bounds of the `z_mct` vector elements out by `epsvi`.
"""
function Rnd_Out_Z_All(z_mct::Vector{SMCg{N,T}},epsvi::S) where {N,S<:AbstractFloat,T<:AbstractFloat}
  epsv::T = convert(T,epsvi)
  return [SMCg{N,T}(z_mct[i].cc+epsv,z_mct[i].cv-epsv,
             z_mct[i].cc_grad, z_mct[i].cv_grad,
             Interval{T}(z_mct[i].Intv.lo-epsv, z_mct[i].Intv.hi+epsv),
             z_mct[i].cnst, z_mct[i].IntvBox,z_mct[i].xref) for i=1:length(z_mct)]
end

"""
    Rnd_Out_H_Intv(z_mct::Vector{SMCg{N,T}},Y_mct::Array{SMCg{N,T},2},epsvi::S)

Rounds the interval bounds of the `z_mct` and `Y_mct` elements out by `epsvi`.
"""
function Rnd_Out_H_All(z_mct::Vector{SMCg{N,T}},Y_mct::Array{SMCg{N,T},2},epsvi::S) where {N,S<:AbstractFloat,T<:AbstractFloat}
  epsv::T = convert(T,epsvi)
  temp1::Vector{SMCg{N,T}} = [SMCg{N,T}(z_mct[i].cc+epsv,z_mct[i].cv-epsv,
                                        z_mct[i].cc_grad, z_mct[i].cv_grad,
                                        Interval{T}(z_mct[i].Intv.lo-epsv, z_mct[i].Intv.hi+epsv),
                                        z_mct[i].cnst, z_mct[i].IntvBox,z_mct[i].xref) for i=1:length(z_mct)]
  temp2::Array{SMCg{N,T},2} = [SMCg{N,T}(Y_mct[i,j].cc+epsv,Y_mct[i,j].cv-epsv,
                                        Y_mct[i,j].cc_grad, Y_mct[i,j].cv_grad,
                                        Interval{T}(Y_mct[i,j].Intv.lo-epsv, Y_mct[i,j].Intv.hi+epsv),
                                        Y_mct[i,j].cnst, Y_mct[i,j].IntvBox,Y_mct[i,j].xref) for i=1:length(z_mct), j=1:length(z_mct)]
  return temp1,temp2
end

"""
    Rnd_Out_H_All(z_mct::Vector{SMCg{N,T}},Y_mct::Array{SMCg{N,T},2},epsvi::S)

Rounds the interval and relaxation bounds of the `z_mct` and `Y_mct` elements out by `epsvi`.
"""
function Rnd_Out_H_Intv(z_mct::Vector{SMCg{N,T}},Y_mct::Array{SMCg{N,T},2},epsvi::S) where {N,S<:AbstractFloat,T<:AbstractFloat}
  epsv::T = convert(T,epsvi)
  temp1::Vector{SMCg{N,T}} = [SMCg(z_mct[i].cc,z_mct[i].cv,
             z_mct[i].cc_grad, z_mct[i].cv_grad,
             @interval(z_mct[i].Intv.lo-epsv, z_mct[i].Intv.hi+epsv),
             z_mct[i].cnst, z_mct[i].IntvBox,z_mct[i].xref) for i=1:length(z_mct)]
  temp2::Array{SMCg{N,T},2} = [SMCg(Y_mct[i,j].cc,Y_mct[i,j].cv,
             Y_mct[i,j].cc_grad, Y_mct[i,j].cv_grad,
             @interval(Y_mct[i,j].Intv.lo-epsv, Y_mct[i,j].Intv.hi+epsv),
             Y_mct[i,j].cnst, Y_mct[i,j].IntvBox,Y_mct[i,j].xref) for i=1:length(z_mct), j=1:length(z_mct)]
  return temp1,temp2
end

"""
    Precondition(hm::Vector{SMCg{N,T}},hJm::Union{Vector{SMCg{N,T}},Array{SMCg{N,T},2}},
                 Y::Union{Vector{T},Array{T,2}},nx::Int64)

Preconditions `hm` and `hJm` by `Y` in place where all dimensions are `nx`.
"""
function Precondition!(hm::Vector{SMCg{N,T}},hJm::Union{Vector{SMCg{N,T}},Array{SMCg{N,T},2}},
                      Y::Union{Vector{T},Array{T,2}},nx::Int64) where {N,T<:AbstractFloat}
  S1::SMCg{N,T},S2::SMCg{N,T} = zero(SMCg{N,T}),zero(SMCg{N,T})
  for i=1:nx
    S2 = zero(SMCg{N,T})
    for j=1:nx
      S1 = zero(SMCg{N,T})
      for k=1:nx
        S1 = S1 + Y[i,k]*hJm[k,j]
      end
      hJm[i,j] = S1
      S2 += Y[i,j]*hm[j]
    end
    hm[i] = S2
  end
  return hm,hJm
end


text:
"Generic PDDL planning problem."
@kwdef mutable struct GenericProblem <: Problem
    name::Symbol # Name of problem
    domain::Symbol # Name of associated domain
    objects::Vector{Const} # List of objects
    objtypes::Dict{Const,Symbol} # Types of objects
    init::Vector{Term} # Predicates that hold in initial state
    goal::Term # Goal formula
    metric::Union{Term,Nothing} # Metric formula
    constraints::Union{Term,Nothing} # Constraints formula
end

GenericProblem(problem::GenericProblem) = copy(problem)
function GenericProblem(problem::Problem)
    return GenericProblem(
        get_name(problem),
        get_domain_name(problem),
        collect(Const, get_objects(problem)),
        Dict{Const,Symbol}(pairs(get_objtypes(state))...),
        collect(Term, get_init_terms(problem)),
        get_goal(problem),
        get_metric(problem),
        get_constraints(problem)
    )
end

function GenericProblem(
    name::Symbol, header::Dict{Symbol,Any}, body::Dict{Symbol,Any}
)
    header = filter(item -> first(item) in fieldnames(GenericProblem), header)
    body = filter(item -> first(item) in fieldnames(GenericProblem), body)
    return GenericProblem(;name=name, header..., body...)
end

function GenericProblem(
    state::State;
    domain=:domain,
    name=Symbol(domain, "-problem"),
    goal::Term=Compound(:and, []),
    metric=nothing,
    constraints=nothing,
)
    objtypes = Dict{Const,Symbol}(pairs(get_objtypes(state))...)
    objects = collect(get_objects(state))
    init = Term[]
    for (name, val) in get_fluents(state)
        if val isa Bool && val == true # Handle Boolean predicates
            term = name
        else # Handle non-Boolean fluents
            val = valterm(val) # Express value as term
            term = Compound(:(==), Term[name, val]) # Assignment expression
        end
        push!(init, term)
    end
    return GenericProblem(Symbol(name), Symbol(domain),
                          objects, objtypes, init,
                          goal, metric, constraints)
end

Base.copy(problem::GenericProblem) = deepcopy(problem)

get_name(problem::GenericProblem) = problem.name

get_domain_name(problem::GenericProblem) = problem.domain

get_objects(problem::GenericProblem) = problem.objects

get_objtypes(problem::GenericProblem) = problem.objtypes

get_init_terms(problem::GenericProblem) = problem.init

get_goal(problem::GenericProblem) = problem.goal

get_metric(problem::GenericProblem) = problem.metric

get_constraints(problem::GenericProblem) = problem.constraints


text:
using ITensors

# This is a short code showing how a combiner
# can be used to "flip" the direction of an Index
i = Index([QN(0) => 2, QN(1) => 3], "i")
j = Index([QN(0) => 2, QN(1) => 3], "j")
A = randomITensor(i, dag(j))
C = combiner(dag(j); tags="jflip", dir=-dir(dag(j)))
@show inds(A)
@show inds(A * C)


text:
@testset "$TEST $G" begin

g = Network(10)
@test graphtype(g) == Network
@test digraphtype(g) == DiNetwork
@test edgetype(g) == IndexedEdge
@test vertextype(g) == Int

g = DiNetwork(10)
@test graphtype(g) == Network
@test digraphtype(g) == DiNetwork
@test edgetype(g) == IndexedEdge
@test vertextype(g) == Int

end # testset


text:
"""
    auth(password[, username])

Authenticate to the server.
"""
auth(password, username=""; client=get_global_client()) = execute(["AUTH", username, password], client)

"""
    select(database)

Change the selected database for the current connection.
"""
select(database; client=get_global_client()) = execute(["SELECT", database], client)

"""
    ping()

Ping the server.
"""
ping(; client=get_global_client()) = execute(["PING"], client)

"""
    flushdb([; async=false])

Remove all keys from the current database.
"""
flushdb(; async=false, client=get_global_client()) = execute(["FLUSHDB", async ? "ASYNC" : ""], client)

"""
    flushall([; async=false])

Remove all keys from all databases.
"""
flushall(; async=false, client=get_global_client()) = execute(["FLUSHALL", async ? "ASYNC" : ""], client)

"""
    quit()

Close the connection.
"""
quit(; client=get_global_client()) = execute(["QUIT"], client)

"""
    set(key, value)

Set the string value of a key.
"""
set(key, value; client=get_global_client()) = execute(["SET", key, value], client)

"""
    setnx(key, value)

Set the value of a key, only if the key does not exist.
"""
setnx(key, value; client=get_global_client()) = execute(["SETNX", key, value], client)

"""
    get(key)

Get the value of a key.
"""
Base.get(key; client=get_global_client()) = execute(["GET", key], client)

"""
    del(key[, keys...])

Delete a key.
"""
del(key, keys...; client=get_global_client()) = execute(["DEL", key, keys...], client)

"""
    exists(key[, keys...])

Determine if a key exists.
"""
exists(key, keys...; client=get_global_client()) = execute(["EXISTS", key, keys...], client)

"""
    hexists(key, field)

Determine if a hash field exists.
"""
hexists(key, field; client=get_global_client()) = execute(["HEXISTS", key, field], client)

"""
    keys(pattern)

Find all keys matching the pattern.

!!! compat "Julia 1.6"
    `Jedis.keys` is not exported as the interface conflicts with `Base.keys`.
"""
keys(pattern; client=get_global_client()) = execute(["KEYS", pattern], client)

"""
    hkeys(key)

Get all fields in a hash.
"""
hkeys(key; client=get_global_client()) = execute(["HKEYS", key], client)

"""
    setex(key, seconds, value)

Set the value and expiration of a key.
"""
setex(key, seconds, value; client=get_global_client()) = execute(["SETEX", key, seconds, value], client)

"""
    expire(key, seconds)

Set a key's tiem to live in seconds.
"""
expire(key, seconds; client=get_global_client()) = execute(["EXPIRE", key, seconds], client)

"""
    ttl(key)

Get the time to live for a key.
"""
ttl(key; client=get_global_client()) = execute(["TTL", key], client)

"""
    multi()

Mark the start of a transaction block.

# Examples
```julia-repl
julia> multi()
"OK"

julia> set("key", "value")
"QUEUED"

julia> get("key")
"QUEUED"

julia> exec()
2-element Array{String,1}:
 "OK"
 "value"
```
"""
multi(; client=get_global_client()) = execute(["MULTI"], client)

"""
    exec()

Execute all commands issued after MULTI.

# Examples
```julia-repl
julia> multi()
"OK"

julia> set("key", "value")
"QUEUED"

julia> get("key")
"QUEUED"

julia> exec()
2-element Array{String,1}:
 "OK"
 "value"
```
"""
exec(; client=get_global_client()) = execute(["EXEC"], client)

"""
    multi_exec(fn::Function)

Execute a MULTI/EXEC transction in a do block.

# Examples
```julia-repl
julia> multi_exec() do
           set("key", "value")
           get("key")
           get("key")
       end
3-element Array{String,1}:
 "OK"
 "value"
 "value"
```
"""
multi_exec(fn::Function; client=get_global_client()) = (multi(; client=client); fn(); exec(; client=client))

"""
    hset(key, field, value)

Set the string value of a hash field.
"""
hset(key, field, value, fields_and_values...; client=get_global_client()) = execute(["HSET", key, field, value, fields_and_values...], client)

"""
    hget(key, field)

Get the value of a hash field.
"""
hget(key, field; client=get_global_client()) = execute(["HGET", key, field], client)

"""
    hgetall(key)

Get all the fields and values in a hash.
"""
hgetall(key; client=get_global_client()) = execute(["HGETALL", key], client)

"""
    hmget(key, field[, fields...])

Get the values of all the given hash fields.
"""
hmget(key, field, fields...; client=get_global_client()) = execute(["HMGET", key, field, fields...], client)

"""
    hdel(key, field[, fields...])

Delete one or more hash fields.
"""
hdel(key, field, fields...; client=get_global_client()) = execute(["HDEL", key, field, fields...], client)

"""
    lpush(key, element[, elements...])

Prepend one or multiple elements to a list.
"""
lpush(key, element, elements...; client=get_global_client()) = execute(["LPUSH", key, element, elements...], client)

"""
    rpush(key, element[, elements...])

Append one or multiple elements to a list.
"""
rpush(key, element, elements...; client=get_global_client()) = execute(["RPUSH", key, element, elements...], client)

"""
    lpos(key, element[, rank, num_matches, len])

Return the index of matching element on a list.
"""
lpos(key, element, rank="", num_matches="", len=""; client=get_global_client()) = execute(["LPOS", key, element, [isempty(rank) ? "" : "RANK", rank]..., [isempty(num_matches) ? "" : "COUNT", num_matches]..., [isempty(len) ? "" : "MAXLEN", len]...], client)

"""
    lrem(key, count, element)

Remove elements from a list.
"""
lrem(key, count, element; client=get_global_client()) = execute(["LREM", key, count, element], client)

"""
    lpop(key)

Remove and get the first element in a list.
"""
lpop(key; client=get_global_client()) = execute(["LPOP", key], client)

"""
    rpop(key)

Remove and get the last element in a list.
"""
rpop(key; client=get_global_client()) = execute(["RPOP", key], client)

"""
    blpop(key[, key...; timeout=0])

Remove and get the first element in a list, or block until one is available.
"""
blpop(key, keys...; timeout=0, client=get_global_client()) = execute(["BLPOP", key, keys..., timeout], client)

"""
    brpop(key[, key...; timeout=0])

Remove and get the last element in a list, or block until one is available.
"""
brpop(key, keys...; timeout=0, client=get_global_client()) = execute(["BRPOP", key, keys..., timeout], client)

"""
    llen(key)

Get the length of a list.
"""
llen(key; client=get_global_client()) = execute(["LLEN", key], client)

"""
    lrange(key, start, stop)

Get a range of elements from a list.
"""
lrange(key, start, stop; client=get_global_client()) = execute(["LRANGE", key, start, stop], client)

"""
    incr(key)

Increment the integer value of a key by one.
"""
incr(key; client=get_global_client()) = execute(["INCR", key], client)

"""
    incrby(key, increment)

Increment the integer value of a key by the given amount.
"""
incrby(key, increment; client=get_global_client()) = execute(["INCRBY", key, increment], client)

"""
    incrbyfloat(key, increment)

Increment the float value of a key by the given amount.
"""
incrbyfloat(key, increment; client=get_global_client()) = execute(["INCRBYFLOAT", key, increment], client)

"""
    hincrby(key, field, increment)

Increment the integer value of a hash field by the given number.
"""
hincrby(key, increment; client=get_global_client()) = execute(["HINCRBY", key, field, increment], client)

"""
    hincrbyfloat(key, field, increment)

Increment the float value of a hash field by the given number.
"""
hincrbyfloat(key, field, increment; client=get_global_client()) = execute(["HINCRBYFLOAT", key, field, increment], client)

"""
    zincrby(key, field, member)

Increment the score of a member in a sorted set.
"""
zincrby(key, field, increment; client=get_global_client()) = execute(["ZINCRBY", key, field, increment], client)


text:
abstract type AbstractFixPoint end

"""
  Solvers.fixpoint(f,x0::Real,method=SSFixPoint())
does a fixpoint iteration convergence on a series of real numbers.
f is a function that should
the following strategies:
 - `SSFixPoint(dampingfactor = 1.0)`: performs succesive substitutions until convergence is met.
α = `dampingfactor` is determines a buffer for each iteration, defined as `x1 = α*f(x1) + (1-α)*x1`

- `Aitken()`: uses Aitken's delta-squared process to accelerate convergence of the series. recommended for harmonic iterates.

"""
function fixpoint end
struct SSFixPoint{T<:Real} <: AbstractFixPoint
    dampingfactor::T
end

SSFixPoint(;dampingfactor=1.0) = SSFixPoint(dampingfactor)

function promote_method(method::SSFixPoint,T)
    return SSFixPoint(T(method.dampingfactor))
end

struct AitkenFixPoint <: AbstractFixPoint end

function promote_method(method::AitkenFixPoint,T)
    return method
end

function convergence(xold,xi,atol,rtol)
    not_finite = false
    for xii in xi
        if !isfinite(xii)
            not_finite = false
            break
        end
    end
    not_finite && return (true,false) #terminate, with nan
    xi == xold && return (true,true) #terminate, with current number
    if xi isa Number
        Δx = abs(xi-xold)
    else
        Δx = norm((xi[i] - xold[i] for i in eachindex(xold,xi)))
    end
    normxi = norm(xi)
    if abs(Δx) < max(atol,normxi*rtol)
        return (true,true) #terminate, with current number
    end
    return (false,false) #keep iterating
end

function fixpoint(f,x0,
    method::AbstractFixPoint = SSFixPoint();
    atol=zero(eltype(x0)),
    rtol=8eps(one(eltype(x0))),
    max_iters=100)
    _,atol,rtol = promote(one(eltype(x0)),atol,rtol)
    method = promote_method(method,eltype(x0))
    return _fixpoint(f,x0,method,atol,rtol,max_iters)
end

function _fixpoint(f::F,
    x0::T,
    method::SSFixPoint,
    atol::T = zero(T),
    rtol::T =8*eps(T),
    max_iters=100) where {F,T<:Real}
    nan = (0*atol)/(0*atol)
    xi = f(x0)
    converged,finite = convergence(x0,xi,atol,rtol)
    converged && return ifelse(finite,xi,nan)
    itercount = 1
    xold = x0
    α = method.dampingfactor
    while itercount < max_iters
        xi = α*f(xi) + (1-α)*xi
        converged,finite = convergence(xold,xi,atol,rtol)
        converged && return ifelse(finite,xi,nan)
        itercount +=1
        xold = xi
    end
    return nan
end

function _fixpoint(f::F,
    x00::T,
    method::AitkenFixPoint,
    atol::T = zero(T),
    rtol::T =8*eps(T),
    max_iters=100) where {F,T}

    nan = (0*atol)/(0*atol)
    itercount = 2
    x3 = x00
    x1 = f(x3)
    converged,finite = convergence(x3,x1,atol,rtol)
    converged && return ifelse(finite,x1,nan)
    x2 = f(x1)
    converged,finite = convergence(x1,x2,atol,rtol)
    converged && return ifelse(finite,x2,nan)

    while itercount < max_iters
        itercount += 1
        λ2 = (x2 - x1)/(x1-x3)
        dx = -(λ2/(1-λ2))*(x2-x1)
        x3 = x2 + dx
        converged,finite = convergence(x2,x3,atol,rtol)
        converged && return ifelse(finite,x3,nan)

        itercount += 1
        x1 = f(x3)
        converged,finite = convergence(x3,x1,atol,rtol)
        converged && return ifelse(finite,x1,nan)

        itercount += 1
        x2 = f(x1)
        converged,finite = convergence(x1,x2,atol,rtol)
        converged && return ifelse(finite,x2,nan)
    end
    return nan
end

function _fixpoint(f!::F,
    x0::X where {X <:AbstractVector{T}},
    method::SSFixPoint,
    atol::T = zero(T),
    rtol::T =8*eps(T),
    max_iters=100) where {F,T<:Real}
    nan = (0*atol)/(0*atol)
    xi = copy(x0)
    xi = f!(xi,x0)
    converged,finite = convergence(x0,xi,atol,rtol)
    if converged
        if finite
            return xi
        else
            xi .= nan
            return xi
        end
    end
    itercount = 1
    α = method.dampingfactor
    xold = copy(x0)
    while itercount < max_iters
        xi = f!(xi,xold)
        xi .*= α
        xi .+= (1 .- α) .* xold
        converged,finite = convergence(xold,xi,atol,rtol)
        if converged
            if finite
                return xi
            else
                xi .= nan
                return xi
            end
        end
        itercount +=1
        xold .= xi
    end
    xi .= nan
    return xi
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    @assert precompile(Tuple{typeof(watch_manifest), String})
    @assert precompile(Tuple{typeof(revise_dir_queued), String})
    @assert precompile(Tuple{TaskThunk})
    @assert precompile(Tuple{typeof(revise)})
    @assert precompile(Tuple{typeof(includet), String})
    # setindex! doesn't fully precompile, but it's still beneficial to do it
    # (it shaves off a bit of the time)
    # See https://github.com/JuliaLang/julia/pull/31466
    @assert precompile(Tuple{typeof(setindex!), ExprsSigs, Nothing, RelocatableExpr})
    @assert precompile(Tuple{typeof(setindex!), ExprsSigs, Vector{Any}, RelocatableExpr})
    @assert precompile(Tuple{typeof(setindex!), ModuleExprsSigs, ExprsSigs, Module})
    @assert precompile(Tuple{typeof(setindex!), Dict{PkgId,PkgData}, PkgData, PkgId})
    @assert precompile(Tuple{typeof(setindex!), Dict{String,WatchList}, WatchList, String})

    MI = CodeTrackingMethodInfo
    @assert precompile(Tuple{typeof(minimal_evaluation!), MI, Core.CodeInfo})
    @assert precompile(Tuple{typeof(methods_by_execution!), Any, MI, DocExprs, Module, Expr})
    @assert precompile(Tuple{typeof(methods_by_execution!), Any, MI, DocExprs, JuliaInterpreter.Frame, Vector{Bool}})
    @assert precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                            NamedTuple{(:skip_include,),Tuple{Bool}},
                            typeof(methods_by_execution!), Function, MI, DocExprs, Module, Expr})
    @assert precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                            NamedTuple{(:define, :skip_include),Tuple{Bool,Bool}},
                            typeof(methods_by_execution!), Function, MI, DocExprs, Module, Expr})
    @assert precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                            NamedTuple{(:define, :skip_include),Tuple{Bool,Bool}},
                            typeof(methods_by_execution!), Function, MI, DocExprs, JuliaInterpreter.Frame, Vector{Bool}})

    @assert precompile(Tuple{typeof(get_def), Method})
    @assert precompile(Tuple{typeof(parse_pkg_files), PkgId})
    @assert precompile(Tuple{typeof(Base.stale_cachefile), String, String})
    @assert precompile(Tuple{typeof(filter_valid_cachefiles), String, Vector{String}})
    @assert precompile(Tuple{typeof(pkg_fileinfo), PkgId})
    @assert precompile(Tuple{typeof(push!), WatchList, Pair{String,PkgId}})
    @assert precompile(Tuple{typeof(init_watching), PkgData, Vector{String}})
    @assert precompile(Tuple{Type{ModuleExprsSigs}, Module})
    @assert precompile(Tuple{Type{FileInfo}, Module, String})
    @assert precompile(Tuple{Type{PkgData}, PkgId})
    @assert precompile(Tuple{typeof(Base._deleteat!), Vector{Tuple{Module,String,Float64}}, Vector{Int}})

    @assert precompile(Tuple{typeof(track), Module, String})
    @assert precompile(Tuple{typeof(add_require), String, Module, String, String, Expr})
    @assert precompile(Tuple{typeof(_add_require), String, Module, String, String, Expr})
    return nothing
end


text:
using Test
import GeoInterface
import ArchGDAL as AG
import GeoFormatTypes as GFT

@testset "test_cookbook_projection.jl" begin
    @testset "Reproject a Geometry" begin
        @testset "Method 1" begin
            AG.importEPSG(2927) do source
                AG.importEPSG(4326) do target
                    AG.createcoordtrans(source, target) do transform
                        AG.fromWKT("POINT (1120351.57 741921.42)") do point
                            @test AG.toWKT(point) ==
                                  "POINT (1120351.57 741921.42)"
                            AG.transform!(point, transform)
                            @test GeoInterface.coordinates(point) ≈
                                  [47.3488070138318, -122.5981499431438]
                        end
                    end
                end
            end
        end

        @testset "Method 2" begin
            AG.importEPSG(2927) do source
                AG.importEPSG(4326) do target
                    AG.createcoordtrans(source, target) do transform
                        xs = [47.348801, 47.348801]
                        ys = [-122.598135, -122.598135]
                        zs = [0.0, 0.0]
                        @test AG.transform!(xs, ys, zs, transform) == true
                        @test xs ≈ [45.151458, 45.151458]
                        @test ys ≈ [-126.863475, -126.863475]
                        @test zs ≈ [0.0, 0.0]
                    end
                end
            end
        end

        @testset "Use reproject" begin
            @testset "reciprocal reprojection of wkt" begin
                wktpoint = GFT.WellKnownText(
                    GFT.Geom(),
                    "POINT (1120351.57 741921.42)",
                )
                result = GFT.WellKnownText(
                    GFT.Geom(),
                    "POINT (47.3488070138318 -122.598149943144)",
                )
                @test AG.reproject(wktpoint, GFT.EPSG(2927), GFT.EPSG(4326)) ==
                      result
                @test convert(
                    AG.Geometry,
                    AG.reproject(result, GFT.EPSG(4326), GFT.EPSG(2927)),
                ) |> GeoInterface.coordinates ≈
                      [1.12035156999967e6, 741921.420000271]
            end
            @testset "reproject vector, vector of vector, or tuple" begin
                coord = [1120351.57, 741921.42]
                @test AG.reproject(coord, GFT.EPSG(2927), nothing) ≈
                      [1120351.57, 741921.42]
                @test AG.reproject(coord, GFT.EPSG(2927), GFT.EPSG(4326)) ≈
                      [47.3488070138318, -122.5981499431438]
                @test AG.reproject([coord], GFT.EPSG(2927), GFT.EPSG(4326)) ≈
                      [[47.3488070138318, -122.5981499431438]]
                coord = (1120351.57, 741921.42)
                @test AG.reproject(
                    coord,
                    GFT.EPSG(2927),
                    GFT.EPSG(4326);
                    order = :compliant,
                ) ≈ [47.3488070138318, -122.5981499431438]
                @test AG.reproject(
                    coord,
                    GFT.EPSG(2927),
                    GFT.EPSG(4326);
                    order = :trad,
                ) ≈ [-122.5981499431438, 47.3488070138318]
                @test AG.reproject(
                    [coord],
                    GFT.EPSG(2927),
                    convert(GFT.WellKnownText, GFT.EPSG(4326));
                    order = :compliant,
                ) ≈ [[47.3488070138318, -122.5981499431438]]
                @test AG.reproject(
                    [coord],
                    GFT.EPSG(2927),
                    convert(GFT.WellKnownText, GFT.EPSG(4326));
                    order = :trad,
                ) ≈ [[-122.5981499431438, 47.3488070138318]]
                # :compliant doesn't work on PROJ axis order, it loses authority information
                @test AG.reproject(
                    [coord],
                    GFT.EPSG(2927),
                    convert(GFT.ProjString, GFT.EPSG(4326));
                    order = :compliant,
                ) ≈ [[-122.5981499431438, 47.3488070138318]]
                @test AG.reproject(
                    [coord],
                    GFT.EPSG(2927),
                    convert(GFT.ProjString, GFT.EPSG(4326));
                    order = :trad,
                ) ≈ [[-122.5981499431438, 47.3488070138318]]
            end
        end
    end

    @testset "Get Projection" begin
        AG.read("ospy/data1/sites.shp") do dataset
            layer = AG.getlayer(dataset, 0)
            spatialref = AG.getspatialref(layer)
            @test AG.toWKT(spatialref)[1:6] == "PROJCS"
            @test AG.toWKT(spatialref, false)[1:6] == "PROJCS"
            @test AG.toPROJ4(spatialref) ==
                  "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"
            @test AG.toXML(spatialref)[1:17] == "<gml:ProjectedCRS"
            @test AG.toMICoordSys(spatialref) ==
                  "Earth Projection 8, 104, \"m\", -111, 0, 0.9996, 500000, 0"
            AG.nextfeature(layer) do feature
                @test AG.toPROJ4(AG.getspatialref(AG.getgeom(feature))) ==
                      "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"
            end
        end

        AG.importEPSG(26912) do spatialref
            proj4str = "+proj=utm +zone=12 +datum=NAD83 +units=m +no_defs"
            @test AG.toPROJ4(spatialref) == proj4str
            @test AG.toWKT(spatialref)[1:6] == "PROJCS"
            AG.morphtoESRI!(spatialref)
            @test AG.toPROJ4(spatialref) == proj4str
            AG.morphfromESRI!(spatialref)
            @test AG.toPROJ4(spatialref) == proj4str
            AG.importEPSGA!(spatialref, 4326)
            @test AG.toPROJ4(spatialref) ==
                  "+proj=longlat +datum=WGS84 +no_defs"
        end

        AG.importEPSGA(4326) do spatialref
            cloneref = AG.clone(spatialref)
            AG.clone(spatialref) do cloneref2
                @test AG.toWKT(cloneref) == AG.toWKT(cloneref2)
            end
            @test AG.toWKT(cloneref) == AG.toWKT(AG.importEPSGA(4326))
            @test AG.toPROJ4(spatialref) ==
                  "+proj=longlat +datum=WGS84 +no_defs"
        end
    end
end


text:
# VSMBasedUnwarping.jl
# Methods for unwarping geometric distortions in MRI images

function unwarp(VSM, distorted, dim)
    if dim == 2
        distorted = switchdim(distorted)
        VSM = switchdim(VSM)
    end

    unwarped = unwarp(VSM, distorted)

    if dim == 2
        unwarped = switchdim(unwarped)
    end
    unwarped
end

switchdim(v) = permutedims(v, [2, 1, (3:ndims(v))...])

unwarp(VSM, distorted) = unwarp!(similar(distorted), VSM, distorted)

function unwarp!(unwarped, VSM, distorted)
    xi = 1:size(distorted, 1)
    for J in CartesianIndices(size(distorted)[4:end])
        for I in CartesianIndices(size(distorted)[2:3])
            xtrue = xi .+ VSM[:,I]
            xregrid = (xtrue[1] .<= xi .<= xtrue[end]) # only use x values inside (no extrapolation)
            unwarped[.!xregrid,I,J] .= 0
            unwarped[xregrid,I,J] .= unwarpline(xtrue, distorted[:,I,J], xi[xregrid])
        end
    end
    unwarped
end

function unwarpline(xtrue, distorted, xnew)
    #TODO try better interpolation than linear
    interpolate((xtrue,), distorted, Gridded(Linear()))(xnew)
end

function getVSM(B0, rbw, dim, threshold = 5.0)
    VSM = B0 ./ (2π * rbw)
    thresholdforward(VSM, -0.9, threshold, dim)
end

function thresholdforward(VSM, tmin, tmax, dim)
    if dim == 2
        VSM = switchdim(VSM)
    end

    deltaVSM = VSM[2:end,:,:] .- VSM[1:(end-1),:,:]
    VSMret = copy(VSM)
    nx = size(VSM, 1)
    for I in CartesianIndices(size(VSM)[2:3])
        for x in 1:(nx-1)
            if deltaVSM[x,I] < tmin || deltaVSM[x,I] > tmax
                if deltaVSM[x,I] < tmin
                    diff = tmin - deltaVSM[x,I]
                else
                    diff = tmax - deltaVSM[x,I]
                end
                VSMret[x+1,I] += diff
                if x + 1 < nx
                    deltaVSM[x+1,I] -= diff
                end
            end
        end
    end

    if dim == 2
        VSMret = switchdim(VSMret)
    end
    VSMret
end

# TODO not properly working
threshold(VSM, tmin, tmax, dim) = (thresholdforward(VSM, tmin, tmax, dim) .+ thresholdforward(VSM[end:-1:1,:,:], -tmax, -tmin, dim))[end:-1:1,:,:] ./ 2.0


text:
export normalize!

"""
    normalize!(S)

Locally normalize the weights of a SPN using Algorithm 1 from Peharz et al.
"""
function normalize!(spn::SumProductNetwork; ϵ = 1e-10)
    αp = ones(length(spn))

    for (nid, node) in enumerate(values(spn))

        if isa(node, Leaf)
            continue
        end

        α = 0.0

        if isa(node, SumNode)
            α = Float64(sum(weights(node)))
            if α < ϵ
                α = ϵ
            end
            T = eltype(node)
            logweights(node)[:] .-= T(log(α))
            logweights(node)[weights(node) .< ϵ] .= T(ϵ)
        elseif isa(node, ProductNode)
            α = αp[nid]
            αp[nid] = 1.
        end

        for fnode in parents(node)
            if isa(fnode, SumNode)
                id = findfirst(children(fnode) .== node)
                @assert id > 0
                T = eltype(fnode)
                logweights(fnode)[id] = logweights(fnode)[id] + T(log(α))
            elseif isa(fnode, ProductNode)
                id = findfirst(nodes .== fnode)
                if id == 0
                    @error("Parent of the following node not found! ", nid)
                end
                @assert id > 0
                αp[id] = α * αp[id]
            end
        end
    end
end


text:
module mrs1
using Elfel.RefShapes: RefShapeInterval, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeInterval) == 1
    qr = quadrature(RefShapeInterval)
    @test npts(qr) == 1
    qr = quadrature(RefShapeInterval, (order = 3,))
    @test npts(qr) == 3
    @test isapprox(param_coords(qr), [-0.774596669241483; 0.0; 0.774596669241483])
    @test isapprox(weights(qr), [0.5555555555555556; 0.888888888888889; 0.5555555555555556])
end
end
using .mrs1
mrs1.test()

module mrs2
using Elfel.RefShapes: RefShapeTriangle, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeTriangle) == 2
    qr = quadrature(RefShapeTriangle)
    @test npts(qr) == 1
    qr = quadrature(RefShapeTriangle, (npts = 3,))
    @test npts(qr) == 3
    @test isapprox(param_coords(qr), [0.6666666666666666 0.16666666666666666; 0.16666666666666666 0.6666666666666666; 0.16666666666666666 0.16666666666666666])
    @test isapprox(weights(qr), [0.16666666666666666; 0.16666666666666666; 0.16666666666666666])
end
end
using .mrs2
mrs2.test()

module mrs3
using Elfel
using Elfel.RefShapes: RefShapeSquare, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeSquare) == 2
    qr = quadrature(RefShapeSquare)
    # @show qr
    @test npts(qr) == 1
    qr = quadrature(RefShapeSquare, (order = 2,))
    @test npts(qr) == 4
    # @show param_coords(qr)[3, :]
    @test isapprox(vec(param_coords(qr)[3, :]), [0.577350269189626, -0.577350269189626])
    # @test isapprox(weights(qr), [0.16666666666666666; 0.16666666666666666; 0.16666666666666666])
end
end
using .mrs3
mrs3.test()

module mrs4
using Elfel
using Elfel.RefShapes: RefShapeSquare, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeSquare) == 2
    qr = quadrature(RefShapeSquare)
    # @show qr
    for o in 1:10
        qr = quadrature(RefShapeSquare, (order = o,))
        @test npts(qr) == o^2
    end
    true
end
end
using .mrs4
mrs4.test()

module mrs5
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeTriangle) == 2
   for n in [1, 3, 4, 6, 7, 9, 12, 13]
        qr = quadrature(RefShapeTriangle, (npts = n,))
        @test npts(qr) == n
    end
    true
end
end
using .mrs5
mrs5.test()

module mrs6
using Elfel.RefShapes: RefShapeTetrahedron, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeTetrahedron) == 3
    qr = quadrature(RefShapeTetrahedron)
    @test npts(qr) == 1
    @test isapprox(param_coords(qr), [0.25 0.25 0.25])
    @test isapprox(weights(qr), [0.16666666666666666])
    qr = quadrature(RefShapeTetrahedron, (npts = 4,))
    @test npts(qr) == 4
    @test isapprox(param_coords(qr),  [0.1381966 0.1381966 0.1381966;
        0.5854102 0.1381966 0.1381966;
        0.1381966  0.5854102 0.1381966;
        0.1381966 0.1381966 0.5854102] )
    @test isapprox(weights(qr), [0.041666666666666664, 0.041666666666666664, 0.041666666666666664, 0.041666666666666664])
    qr = quadrature(RefShapeTetrahedron, (npts = 5,))
    @test npts(qr) == 5
    @test isapprox(param_coords(qr), [0.25 0.25 0.25;
        0.5 0.16666666666666666 0.16666666666666666;
        0.16666666666666666 0.5 0.16666666666666666;
        0.16666666666666666 0.16666666666666666 0.5;
        0.16666666666666666 0.16666666666666666 0.16666666666666666])
    @test isapprox(weights(qr), [-0.13333333333333333, 0.075, 0.075, 0.075, 0.075] )
end
end
using .mrs6
mrs6.test()

module mrs7
using Elfel.RefShapes: RefShapeInterval, manifdim
using Elfel.RefShapes: quadrature, npts, param_coords, weights
using Test
function test()
    @test manifdim(RefShapeInterval) == 1
    result_true = exp(1.0) - exp(-1.0)
    reference_results = [2.0
    , 2.3426960879097307
    , 2.350336928680011
    , 2.3504020921563744
    , 2.350402386462827
    , 2.350402387286035
    , 2.350402387287601
    , 2.3504023872876028
    , 2.3504023872876036
    , 2.350402387287604   ]
    results = Float64[]
    for N in 1:10
        qr = quadrature(RefShapeInterval, (order = N,))
        push!(results, sum(exp.(param_coords(qr)) .* weights(qr)))
    end
    @test isapprox(results, reference_results)
    end
end
using .mrs7
mrs7.test()


text: Base.promote_rule(::Type{<:Combiner},StorageT::Type{<:CuDense}) = StorageT


text:
function ShowPSUWindow(psu, psu_conf, rev_state_dict, refresh_cnt, base)
	CImGui.Begin("PST3201 Power supply")
	# DRAW DISPLAY INFO
	draw_psu_info("CH1", " ", psu_conf.C1.volt_meas, "V",
					psu_conf.C1.curr_meas, "A")
	draw_psu_info("CH2", " ", psu_conf.C2.volt_meas, "V",
					psu_conf.C2.curr_meas, "A")
	draw_psu_info("CH3", " ", psu_conf.C3.volt_meas, "V",
					psu_conf.C3.curr_meas, "A")

	CImGui.Text("")

	draw_psu_info("SET", "   ", "CH1", "    ",
					"CH2", "    ","CH3")
	draw_psu_info("Volt.", psu_conf.C1.volt_set, "V",
					psu_conf.C2.volt_set, "V",
					psu_conf.C3.volt_set, "V")
	draw_psu_info("Curr.", psu_conf.C1.curr_set, "A",
					psu_conf.C2.curr_set, "A",
					psu_conf.C3.curr_set, "A")
	draw_psu_info("O.V.P", psu_conf.C1.ovp_set, "V",
					psu_conf.C2.ovp_set, "V",
					psu_conf.C3.ovp_set, "V")
	draw_psu_info("O.C.P", " ", psu_conf.C1.ocp_set,
					"    ", psu_conf.C2.ocp_set,
					"    ", psu_conf.C3.ocp_set,)
	CImGui.Text("")

	# DRAW BUTTONS
	CImGui.Text("OVP RST ")
	CImGui.SameLine(), CImGui.Text("  CH1  ")
	CImGui.SameLine(), CImGui.Text(" CH2  ")
	CImGui.SameLine(), CImGui.Text("  CH3  ")

	CImGui.Button("  OCP  ") && begin
		if psu_conf.shift_btn.state == true
			reset_protections(psu);
		else # toggle OCP
			psu_conf.ocp_state = get_curr_protection(psu, psu_conf.crt_chan);
			crt_state = rev_state_dict[psu_conf.ocp_state];
			set_curr_protection(psu, psu_conf.crt_chan,crt_state);
		end
	end

	CImGui.SameLine()
	CImGui.Button(" VSET ") && begin
		psu_conf.shift_btn.state ? psu_conf.crt_chan = "C1" : psu_conf.crt_func = "voltage"
	end

	CImGui.SameLine()
	CImGui.Button("ISET ") && begin
		psu_conf.shift_btn.state ? psu_conf.crt_chan = "C2" : psu_conf.crt_func = "current"
	end

	CImGui.SameLine()
	CImGui.Button("OVPSET ") && begin
		psu_conf.shift_btn.state ? psu_conf.crt_chan = "C3" : psu_conf.crt_func = "OVP"
	end

	# Draw SHIFT
	psu_conf.shift_btn, pressed = draw_toggle_button(psu_conf.shift_btn)

	# Draw OUTPUT
	CImGui.SameLine()
	psu_conf.output_btn, pressed = draw_toggle_button(psu_conf.output_btn)
	pressed && set_outp(psu, psu_conf.crt_chan, rev_state_dict[psu_conf.output_state])

	# Draw ENTER
	CImGui.SameLine()
	CImGui.Button("ENTER") && begin
		if psu_conf.crt_func == "voltage"
			set_source_lev(psu, psu_conf.crt_chan, psu_conf.crt_value);
		elseif psu_conf.crt_func == "current"
			set_max_curr(psu, psu_conf.crt_chan, psu_conf.crt_value);
		elseif psu_conf.crt_func == "OVP"
			set_volt_protection(psu, psu_conf.crt_chan, psu_conf.crt_value);
		else
			@error "Brace yourself!"
		end
	end

	# Draw REFRESH
	CImGui.SameLine()
	psu_conf.refresh_btn, pressed = draw_toggle_button(psu_conf.refresh_btn)
	psu_conf.refresh_btn.state && @async update_psu_conf!(psu_conf, psu, refresh_cnt, base)

	# Draw INPUT BOX
	@c CImGui.InputDouble("", &psu_conf.crt_value, 0.01, 1.0, "%.4f")

	CImGui.End()
	return psu_conf
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
# Random number generators for Julia
# This is a slightly modified and reduced version of a Julia module by John D. Cook.
# http://www.johndcook.com/julia_rng.html

module Random

inverse_gamma(shape, rate) = 1/gamma(shape, rate)

function gamma(shape, rate)
    # Return a random sample from a gamma distribution
    # with pdf \propto x^{shape-1}exp(-rate*x)

    if shape <= 0.0
        error("Shape parameter must be positive, but shape = $shape")
    end
    if rate <= 0.0
        error("Rate parameter must be positive, but rate = $rate")
    end

    ## Implementation based on "A Simple Method for Generating Gamma Variables"
    ## by George Marsaglia and Wai Wan Tsang.
    ## ACM Transactions on Mathematical Software
    ## Vol 26, No 3, September 2000, pages 363-372.

    if shape >= 1.0
        d = shape - 1.0/3.0
        c = 1.0/sqrt(9.0*d)
        while true
            x = randn()
            v = 1.0 + c*x
            while v <= 0.0
                x = randn()
                v = 1.0 + c*x
            end
            v = v*v*v
            u = rand()
            xsq = x*x
            if u < 1.0 -.0331*xsq*xsq || log(u) < 0.5*xsq + d*(1.0 - v + log(v))
                return d*v/rate
            end
        end
    else
        g = gamma(shape+1.0, 1.0)
        w = rand()
        return g*(w^(1.0/shape))/rate
    end
end

## return a random sample from a chi square distribution
## with the specified degrees of freedom
function chi_square(dof)
    return gamma(0.5*dof, 0.5)
    # This was corrected from gamma(0.5, 2.0*dof) in the original code
end

## return a random sample from a beta distribution
function beta(a, b)
    if (a <= 0) || (b <= 0)
        error("Beta parameters must be positive")
    end

    ## There are more efficient methods for generating beta samples.
    ## However such methods are a little more efficient and much more complicated.
    ## For an explanation of why the following method works, see
    ## http://www.johndcook.com/distribution_chart.html#gamma_beta

    u = gamma(a, 1.0)
    v = gamma(b, 1.0)
    return u / (u + v)
end

end


text:
# y ~ @gaussian_noise_around x
macro gaussian_noise_around(x)
    :(normal($(esc(x)), 1.0))
end

# y = @traced_gaussian_noise_around(x, :y)
macro traced_gaussian_noise_around(x, y)
    :(
        Gen.@trace(
            normal( $(esc(x)) , 1.0),
            $y
        )
    )
end

# @add_gaussian_noise x y  # means y ~ normal(x, 1.)
macro add_gaussian_noise(x, y)
    Expr(:call, esc(:~), esc(y), :(normal($(esc(x)), 1.0)))
    #:($(esc(y)) $(esc(:~)) normal($(esc(x)), 1.0))
end

macro set_x_to_normal_tilde(mean)
    Expr(:call, esc(:~), esc(:x), :(normal($(esc(mean)), 1)))
    #:($(esc(:x)) $(esc(:~)) normal($(esc(mean)), 1))
end

macro set_x_to_normal_trace(mean)
    :($(esc(:x)) = Gen.@trace(normal($(esc(mean)), 1), $(QuoteNode(:x))))
end

macro set_x_to_normal_traceexpr(mean)
    :($(esc(:x)) = $(Expr(:gentrace,
        :(normal($(esc(mean)), 1.)),
        QuoteNode(:x)
    )))
end

@testset "macros in dynamic DSL" begin
    @gen function foo1()
        x ~ exponential(1)
        y ~ @gaussian_noise_around x # y ~ normal(x, 1)
    end
    trace = simulate(foo1, ())
    @test trace[:y] != trace[:x]

    @gen function foo2()
        x ~ exponential(1)
        y = Gen.@trace(@gaussian_noise_around(x), :y) # y = @trace(normal(x, 1), :y)
    end
    trace = simulate(foo2, ())
    @test trace[:y] != trace[:x]

    @gen function foo3()
        x ~ exponential(1)
        y = @traced_gaussian_noise_around x :y # y = @trace(normal(x, 1), :y)
    end
    trace = simulate(foo3, ())
    @test trace[:y] != trace[:x]

    @gen function foo4()
        x ~ exponential(1)
        @add_gaussian_noise x y # y ~ normal(x, 1)
    end
    trace = simulate(foo4, ())
    @test trace[:y] != trace[:x]

    @gen function foo5()
        y ~ exponential(1)
        @set_x_to_normal_tilde y # x ~ normal(y, 1)
    end
    trace = simulate(foo5, ())
    @test trace[:y] != trace[:x]

    @gen function foo6()
        y ~ exponential(1)
        @set_x_to_normal_trace y # x = Gen.@trace(normal(y, 1), :y)
    end
    trace = simulate(foo6, ())
    @test trace[:y] != trace[:x]

    @gen function foo7()
        y ~ exponential(1)
        @set_x_to_normal_traceexpr y # x = Gen.@trace(normal(y, 1), :y) USING EXPR(:GENTRACE)
    end
    trace = simulate(foo7, ())
    @test trace[:y] != trace[:x]
end

############################################
# Define static functions at global scope: #
############################################

@gen (static) function foo1()
    x ~ exponential(1)
    y ~ @gaussian_noise_around x # y ~ normal(x, 1)
end
@gen (static) function foo2()
    x ~ exponential(1)
    y = Gen.@trace(@gaussian_noise_around(x), :y) # y = @trace(normal(x, 1), :y)
end
@gen (static) function foo3()
    x ~ exponential(1)
    y = @traced_gaussian_noise_around x :y # y = @trace(normal(x, 1), :y)
end
@gen (static) function foo4()
    x ~ exponential(1)
    @add_gaussian_noise x y # y ~ normal(x, 1)
end
@gen (static) function foo5()
    y ~ exponential(1)
    @set_x_to_normal_tilde y # x ~ normal(y, 1)
end
@gen (static) function foo6()
    y ~ exponential(1)
    @set_x_to_normal_trace y # x = Gen.@trace(normal(y, 1), :y)
end
@gen (static) function foo7()
    y ~ exponential(1)
    @set_x_to_normal_traceexpr y # x = Gen.@trace(normal(y, 1), :y) USING EXPR(:GENTRACE)
end

@load_generated_functions()

@testset "macros in static DSL" begin
    trace = simulate(foo1, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo2, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo3, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo4, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo5, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo6, ())
    @test trace[:y] != trace[:x]

    trace = simulate(foo7, ())
    @test trace[:y] != trace[:x]
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using Flux.Optimise:ADAM, update!
using Zygote

"""
@with_kw struct pplm
    method::String="BoW"
    perturb::String="hidden"     # hidden or past -> hidden support BoW only without gradient based change
    bow_list::Vector{String}=["military"]
    discrim::String="toxicity"
    embed_size::Int=768
    target_class_id=-1
    file_name::String=""
    path::String=""
    stepsize::Float32=0.01
    max_length::Int=100
    num_iterations::Int=2        # more the number of iterations, more updates, more time to update
    top_k::Int=50
    top_p::Float32=0.8
    temperature::Float32=1.1
    fusion_gm_scale::Float32=0.9
    fusion_kl_scale::Float32=0.01
    cuda::Tuple{Bool, Int64}=(CUDA.has_cuda(), device_id)
    window_length::Int=0          # window length 0 corresponds to infinite length
    gamma::Float32=1.5
end

Struct to contain all the hyperparameters required for PPLM based generation.

Example:

```julia
# bow model
args = pplm(method="BoW", bow_list=["legal"], num_iterations=3)

# discriminator model
args = pplm(method="Discrim", perturb="past", stepsize=0.02)
```

"""
@with_kw struct pplm
    method::String="BoW"
    perturb::String="hidden"     # hidden or past -> hidden support BoW only without gradient based change
    bow_list::Vector{String}=["military"]
    discrim::String="toxicity"
    embed_size::Int=768
    target_class_id=-1
    file_name::String=""
    path::String=""
    stepsize::Float32=0.01
    max_length::Int=100
    num_iterations::Int=2        # more the number of iterations, more updates, more time to update
    top_k::Int=50
    top_p::Float32=0.8
    temperature::Float32=1.1
    fusion_gm_scale::Float32=0.9
    fusion_kl_scale::Float32=0.01
    cuda::Tuple{Bool, Int64}=(CUDA.has_cuda(), device_id)
    window_length::Int=0          # window length 0 corresponds to infinite length
    gamma::Float32=1.5
end

# this is gradient based perturbation technique, supports only BoW
"""
    perturb_probs(probs, tokenizer, args)

Perturb probabilities `probs` based on provided Bag of Words list (as given with `args`). This function is supported only for BoW model.
"""
function perturb_probs(probs, tokenizer, args)
    global device
    opt=ADAM(args.stepsize)
    bow_indices, _ = PPLM.get_bow_indices(args.bow_list, tokenizer)
    _, bow_ohe = PPLM.build_bow_ohe(bow_indices, tokenizer)
    bow_all = bow_ohe[1]
    for i in 2:length(bow_ohe)
        bow_all = vcat(bow_all, bow_ohe)
    end
    bow_all = bow_all |> device

    function bow_loss(probs)
        loss_words = bow_all*probs
        loss = -log(sum(loss_words))
        loss
    end

    probs_ = deepcopy(probs)
    ps = params(probs_)
    kl_scale = args.fusion_kl_scale
    for i in 1:args.num_iterations
        _, back = Zygote.pullback(ps) do
               loss = bow_loss(probs_)
               loss_2 = kl_scale * Flux.kldivergence(probs_, probs)
               loss + loss_2
        end
        grads = back(1f0)
        update!(opt, ps, grads)
    end
    probs_ = max.(probs_, 0)

    return probs_
end

# need some abstration here
"""
    perturb_hidden_bow(hidden, model, tokenizer, args)

Perturb hidden states `hidden` based on provided Bag of Words list (as given with `args`). The perturbation is primarily based on the gradient calculated over losses evaluated over desired Bag of Words and KL Divergence from original token.

Also checkout [`perturb_hidden_discrim`](@ref PPLM.perturb_hidden_discrim)
"""
function perturb_hidden_bow(hidden, model, tokenizer, args)
    global device
    bow_indices, _ = get_bow_indices(args.bow_list, tokenizer)
    _, bow_ohe = build_bow_ohe(bow_indices, tokenizer)

    bow_all = bow_ohe[1]
    for i in 2:length(bow_ohe)
        bow_all = vcat(bow_all, bow_ohe)
    end
    bow_all = bow_all |> device

    kl_scale=args.fusion_kl_scale
    opt = ADAM(args.stepsize)
    p = temp_softmax(model.lm_head(hidden)[:, end, :]; t=args.temperature)

    function bow_loss(probs)
        loss_words = bow_all*probs
        loss = -log(sum(loss_words))
        loss
    end

    new_hidden = deepcopy(hidden)
    ps = params(new_hidden)
    for i in 1:args.num_iterations
        _, back = Zygote.pullback(ps) do
                logits = model.lm_head(new_hidden)[:, end, :]
                probs = temp_softmax(logits; t=args.temperature)
                loss_1 = bow_loss(probs)
                loss_2 = kl_scale * Flux.kldivergence(probs, p)
                loss_1 + loss_2
            end
        grads = back(1f0)
        update!(opt, ps, grads)
    end
    return new_hidden
end

"""
    perturb_past_bow(model, prev, past, original_probs, args)

Perturb past key values `prev` based on provided Bag of Words list (as given with `args`). The perturbation is primarily based on the gradient calculated over losses evaluated over desired Bag of Words and KL Divergence from original token.

Also checkout [`perturb_past_discrim`](@ref PPLM.perturb_past_discrim)
"""
function perturb_past_bow(model, prev, past, original_probs, tokenizer, args)
    global device
    bow_indices, _ = PPLM.get_bow_indices(args.bow_list, tokenizer)
    _, bow_ohe = PPLM.build_bow_ohe(bow_indices, tokenizer)
    bow_all = bow_ohe[1]
    for i in 2:length(bow_ohe)
        bow_all = vcat(bow_all, bow_ohe)
    end
    bow_all =bow_all|> device

    kl_scale=args.fusion_kl_scale
    opt = ADAM(args.stepsize)

    function bow_loss(probs)
        loss_words = bow_all*probs
        loss = -log(sum(loss_words))
        loss
    end

    pert_past = deepcopy(past) |> device
    ps = params(pert_past)
    for i in 1:args.num_iterations
        _, back = Zygote.pullback(ps) do
                output = model(prev; past_key_values=pert_past,
                                        output_attentions=false,
                                        output_hidden_states=true,
                                        use_cache=true);
                hidden = output.hidden_states[end]
                logits = model.lm_head(hidden)[:, end, :]
                probs = PPLM.temp_softmax(logits; t=args.temperature)
                loss_1 = bow_loss(probs)
                loss_2 = kl_scale * Flux.kldivergence(probs, original_probs)        # throwing error
                loss_1 + loss_2
            end
        grads = back(1f0)
        update!(opt, ps, grads)
    end
    return pert_past
end

# need some abstration here
"""
    perturb_hidden_discrim(hidden, model, tokenizer, args)

Perturb hidden states `hidden` based on provided Discriminator (as given with `args`). The perturbation is primarily based on the gradient calculated over losses evaluated over desired Discriminator attribute and KL Divergence from original token.

Also checkout [`perturb_hidden_bow`](@ref PPLM.perturb_hidden_bow)
"""
function perturb_hidden_discrim(hidden, model, tokenizer, args)
    global device
    classifier, config_metadata = ClassifierHead(;load_from_pretrained=true, discrim=args.discrim, path=args.path)
    classifier = classifier|> device
    if args.target_class_id ==-1
        y_label = config_metadata.default_class
    else
        y_label = args.target_class_id
    end
    y_one_hot = Flux.onehotbatch([y_label], 1:config_metadata["class_size"]) |> device

    kl_scale=args.fusion_kl_scale
    opt = ADAM(args.stepsize)
    p = temp_softmax(model.lm_head(hidden)[:, end, :]; t=args.temperature)

    new_hidden = deepcopy(hidden)
    ps = params(new_hidden)
    for i in 1:args.num_iterations
        _, back = Zygote.pullback(ps) do
                hidden_sum = sum(new_hidden; dims=2)
                logits_1 = classifier(hidden_sum)
                if length(classifier.linear_layer.bias)==1
                    y_label = y_label |> device
                    loss_1 = Flux.logitbinarycrossentropy(logits_1, y_label)
                else
                    loss_1 = Flux.logitcrossentropy(logits_1, y_one_hot)
                end
                logits = model.lm_head(new_hidden)[:, end, :]
                probs = temp_softmax(logits; t=args.temperature)
                loss_2 = kl_scale * Flux.kldivergence(probs, p)

                loss_1 + loss_2
            end
        grads = back(1f0)
        update!(opt, ps, grads)
    end
    return new_hidden
end

# TODO add code to support horizontal length, futuristic perturbation
"""
    perturb_past_discrim(model, prev, past, original_probs, args)

Perturb past key values `prev` based on provided Discriminator (as given with `args`). The perturbation is primarily based on the gradient calculated over losses evaluated over desired Discriminator attribute and KL Divergence from original token.

Also checkout [`perturb_past_bow`](@ref PPLM.perturb_past_bow)
"""
function perturb_past_discrim(model, prev, past, original_probs, args)
    global device
    classifier, config_metadata = ClassifierHead(;load_from_pretrained=true, discrim=args.discrim, path=args.path)
    classifier = classifier|> device
    if args.target_class_id ==-1
        y_label = config_metadata.default_class
    else
        y_label = args.target_class_id
    end
    y_one_hot = Flux.onehotbatch([y_label], 1:config_metadata["class_size"]) |> device

    kl_scale=args.fusion_kl_scale
    opt = ADAM(args.stepsize)

    pert_past = deepcopy(past) |> device
    ps = params(pert_past)
    for i in 1:args.num_iterations
        _, back = Zygote.pullback(ps) do
                output = model(prev; past_key_values=pert_past,
                                        output_attentions=false,
                                        output_hidden_states=true,
                                        use_cache=true);
                hidden = output.hidden_states[end]
                hidden_sum = sum(hidden; dims=2)
                logits_1 = classifier(hidden_sum)
                if length(classifier.linear_layer.bias)==1
                    y_label = y_label |> device
                    loss_1 = Flux.logitbinarycrossentropy(logits_1, y_label)
                else
                    loss_1 = Flux.logitcrossentropy(logits_1, y_one_hot)
                end
                logits = model.lm_head(hidden)[:, end, :]
                probs = temp_softmax(logits; t=args.temperature)
                loss_2 = kl_scale * Flux.kldivergence(probs, original_probs)        # throwing error
                loss_1 + loss_2
            end
        grads = back(1f0)
        update!(opt, ps, grads)
    end
    return pert_past
end


text:
#=
  From https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  We now consider yet another way to model dice, using fair ones only. This representation
  allows for convenient use of the results in arithmetic expressions, e.g., to add up the
  results from several dice. We query for the probabilities of the possible sums we can get
  from two dice given that the first number is even and the second odd.
  """

  See ~/blog/rolling_dice3.blog
      ~/psi/rolling_dice3.psi
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function rolling_dice3()

    # d1 ~ DiscreteUniform(1,6)
    # d2 ~ DiscreteUniform(1,6)
    num_dice = 2
    dice = tzeros(2)
    for d in 1:num_dice
        dice[d] ~ DiscreteUniform(1,6)
    end

    s = sum(dice)
    # s ~ Dirac(sum(dice))
    # Dirac(s == 7)

    # Nice!
    function odd(d)
        return d % 2 == 1
    end

    # Nicer!
    function even(d)
        return !odd(d)
    end

    true ~ Dirac(even(dice[1]))
    true ~ Dirac(odd(dice[2]))

    return s
end

model = rolling_dice3()

num_chains = 2

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, MH(), 10_000)
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 1000, num_chains)
# chains = sample(model, PG(15), 1_000)

# chains = sample(model, IS(), MCMCThreads(), 1000, num_chains)
# chains = sample(model, IS(), 10_000)

# chains = sample(model, SMC(), 10_000)

# chains = sample(model, HMC(0.1,10), 10_000) # nope

# chains = sample(model, NUTS(1000,0.65), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, Gibbs(MH(:zlabels),NUTS(1000,0.65,:m,:b,:sigma)), MCMCThreads(), 40_000, num_chains)

display(chains)

show_var_dist_pct(chains,Symbol("dice[1]")) # Note: :"dice[1]" DON't work!
show_var_dist_pct(chains,Symbol("dice[2]"))
# show_var_dist_pct(chains,Symbol("s"))

println("\nd1 + d2 (s):")
show_var_dist_pct(generated_quantities(model,chains),20)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using LLLplus, Documenter, MacroTools

makedocs(
         sitename="LLLplus.jl",
         pages = [
                  "Home" => "index.md",
                  "Functions" => "functions.md"],
         format = Documenter.HTML(prettyurls = haskey(ENV, "CI")))


text:
module PrintFileTree

export printfiletree

"""
    printfiletree()
    printfiletree(root)

Like the unix utility `tree` (https://linux.die.net/man/1/tree).

Prints complete recursive directory structure of root and all its contents.
"""
function printfiletree(root=".")
    println(root);
    d,f = printfiletree_helper(root)
    println("\n$d directories, $f files")
end
function printfiletree_helper(root, depth=0, opendirs=[true], dirscount=fill(0), filescount=fill(0))
    files = readdir(root)
    for (i,f) in enumerate(files)
        startswith(f, ".") && continue
        lastitem = (i == length(files))
        lastitem && (opendirs[end] = false)
        for p in opendirs[1:end-1]
            print(p ? "│   " : "    ")
        end
        println("$(lastitem ? "└" : "├")── " * f) # path
        path = joinpath(root, f)

        if isdir(path)
            dirscount[] += 1
            push!(opendirs, true)
            printfiletree_helper(path, depth+1, opendirs, dirscount, filescount)
            pop!(opendirs)
        else
            filescount[] += 1
        end
    end
    dirscount[], filescount[]
end

end # module


text:
using Plots, Printf
# pyplot()
do_viz = true

@views function sia_2D()
    # physics
    lx, ly = 10.0, 10.0
    n      = 3.0
    nt     = 5000
    # numerics
    nx, ny = 128, 128
    nout   = 100
    # Derived numerics
    dx, dy = lx/nx, ly/ny
    xc     = LinRange(dx/2, lx-dx/2, nx)
    yc     = LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qx     = zeros(nx-1,ny-2)
    qy     = zeros(nx-2,ny-1)
    # Initial condition
    b      = 0.5*ones(nx  ,ny  )
    H      = exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2)
    rad    = (xc.-lx./2.0).^2 .+(yc.-ly./2.0)'.^2; b[rad.>lx/4] .= -0.5
    # Time loop
    for it = 1:nt
        if (it==11) global t0 = Base.time() end
        qx  .= .-(0.5.*(H[1:end-1,2:end-1].+H[2:end,2:end-1])).^n .*diff(H[:,2:end-1], dims=1)./dx
        qy  .= .-(0.5.*(H[2:end-1,1:end-1].+H[2:end-1,2:end])).^n .*diff(H[2:end-1,:], dims=2)./dy
        dt   = min(dx^2,dy^2)/maximum(H)^n/4.1/4.0
        H[2:end-1,2:end-1] .= max.(0.0, H[2:end-1,2:end-1] .+ dt.*(.-diff(qx, dims=1)./dx .-diff(qy, dims=2)./dy .+ b[2:end-1,2:end-1]))
        # Visualise
        if mod(it,nout)==0 && do_viz
            p1 = heatmap(xc, yc, H', xlabel="lx", ylabel="ly", title="shallow ice, it=$it", aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), c=:viridis)
            p2 = plot(xc, H[:,Int(round(ny/2))], xlabel="lx", ylabel="height", xlims=(xc[1], xc[end]), ylims=(0., 1.2), legend=false, framestyle=:box)
            display(plot(p1, p2, layout=(2, 1)))
        end
    end
    time_s = (Base.time()-t0)
    @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(nt-10)), sigdigits=2))
    return
end

@time sia_2D()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
# In v0.7+ the idea is that `convert`` is a shortcut to a "safe" subset of the constructors
# for a type, that guarantees the resulting type. The constructor is the central machinery
# to instantiate types and convert between instances. (For parametric types, unless overridden,
# you get an implicit internal constructor without parameters, so no need to define that externally)

Base.convert(::Type{T}, l::T) where T<:Sublat = l
Base.convert(::Type{T}, l::Sublat) where T<:Sublat = T(l)

Base.convert(::Type{T}, l::T) where T<:Bravais = l
Base.convert(::Type{T}, l::Bravais) where T<:Bravais = T(l)

# Constructors for conversion
Sublat{E,T,V}(s::Sublat, name = s.name) where {E,T,V<:Vector} =
    Sublat([padright(site, zero(T), Val(E)) for site in s.sites], name)

# We need this to promote different sublats into common dimensionality and type to combine
# into a lattice, while neglecting orbital dimension
Base.promote(ss::Sublat{E,T}...) where {E,T} = ss
function Base.promote_rule(::Type{Sublat{E1,T1,Vector{SVector{E1,T1}}}}, ::Type{Sublat{E2,T2,Vector{SVector{E2,T2}}}}) where {E1,E2,T1,T2}
    E´ = max(E1, E2)
    T´ = promote_type(T1, T2)
    return Sublat{E´, T´, Vector{SVector{E´,T´}}}
end

Bravais{E,L,T}(b::Bravais) where {E,L,T} = Bravais(padtotype(b.matrix, SMatrix{E,L,T}))

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using Pkg

Pkg.activate(@__DIR__)

using HTTP
import JSON3
using Dates

@info "starting listener"
HTTP.listen("0.0.0.0", 8080, reuseaddr=true) do http
   HTTP.setheader(http, "Server" => "Julia-HTTP")
   HTTP.setheader(http, "Date" => Dates.format(Dates.now(), Dates.RFC1123Format) * " GMT")
   if endswith(http.message.target, "/plaintext")
      HTTP.setheader(http, "Content-Type" => "text/plain")
      HTTP.setstatus(http, 200)
      startwrite(http)
      write(http, "Hello, World!")
   elseif endswith(http.message.target, "/json")
      HTTP.setheader(http, "Content-Type" => "application/json")
      HTTP.setstatus(http, 200)
      startwrite(http)
      JSON3.write(http, (;message = "Hello, World!"))
   else
       HTTP.setstatus(http, 404)
       startwrite(http)
       write(http, "Not Found")
   end
end

max_stars_repo_name: tommilligan/FrameworkBenchmarks

max_issues_repo_name: tommilligan/FrameworkBenchmarks

max_forks_repo_name: tommilligan/FrameworkBenchmarks


text:
################################################################################
#
# RBC_VFI.jl solves the RBC model using value function iteration
#
# Shane McMiken
# June 2021
# Boston College
#
################################################################################

#--------------------------------#
#         PREAMBLE
#--------------------------------#

using Plots
using Distributions

#--------------------------------#
#     RBC MODEL
#--------------------------------#
#
# MODEL EQUATIONS
# 1. 1 - β [ (c^η/cp^η) (1 + rp - δ) ]
# 2. w - h^ψ c^η
# 3. (1-δ) k + i - kp
# 4. c + i - y
# 5. y - a * k^α * h^(1-α)
# 6. w - (1-α)*(y/h)
# 7. r - α * (y/k)
# 8. log(ap) - ρ*log(a)

# model parameters
β   = 0.95  # discount rate
η   = 0.8   # relative risk aversion
ψ   = 0.2   # inverse Frisch elasticity
α   = 0.3   # capital share of income
δ   = 0.05  # depreciation rate
ρ   = 0.9   # persistence parameter
σ   = 0.02  # Standard deviation of TFP shocks

# steady state
abar  = 1.0
rbar  = β^(-1)-(1-δ)
khbar = (rbar/α)^(1/(α-1))
chbar = abar*khbar^α-δ*khbar
hbar  = ((1-α)*abar*khbar^α/(abar*khbar^α-δ*khbar)^η)^(1/(ψ+η))
kbar  = khbar*hbar
cbar  = chbar*hbar
wbar  = (1-α)*abar*khbar^α
ybar  = kbar^α*hbar^(1-α)
ibar  = δ*kbar

#--------------------------------#
#  GRIDS and TRANSITION MATRIX
#--------------------------------#
# Steady state value function
Vbar = 1.0/(1.0-β)*(cbar^(1.0-η)/(1-η) - hbar^(1.0+ψ)/(1.0+ψ))

# Hours grid - in levels
nh = 150
hmin = 0.7*hbar
hmax = 1.3*hbar
vecsize = nh
step = (hmax - hmin)/(vecsize-1)
hgrid = [hmin + (i-1)*step for i∈1:nh]

# Capital grid - in logs
nk = 150
kmin = 0.7*kbar
kmax = 1.3*kbar
vecsize = nk
step = (kmax - kmin)/(vecsize-1)
kgrid  = [kmin + (i-1)*step for i∈1:nk]

# Technology grid - in logs - Tauchen (1986)
cover = 3
na = 2*cover + 1
m  = 1.5

function MyTauchen(σ_a::Float64,ρ_a::Float64,na::Int64,m::Float64)

  vecsize = na
  σ_y = sqrt(σ_a^2 / (1-(ρ_a^2))) # Stationary distribution variance
  step = 2 * m * σ_y / (vecsize-1)
  agrid  = [-m*σ_y + (i-1)*step for i∈1:na]

  # Calculate transition probabilities, AR(1) - Tauchen
  P = zeros(na, na)

  mm = agrid[2] - agrid[1]; # re-compute step
  for j = 1:na
    for k = 1:na
      if(k == 1) # first column
        P[j, k] = cdf(Normal(), (agrid[k] - ρ_a*agrid[j] + (mm/2))/σ_a)
      elseif(k == na) # second-last column
        P[j, k] = 1 - cdf(Normal(), (agrid[k] - ρ_a*agrid[j] - (mm/2))/σ_a)
      else
        P[j, k] = cdf(Normal(), (agrid[k] - ρ_a*agrid[j] + (mm/2))/σ_a) - cdf(Normal(), (agrid[k] - ρ_a*agrid[j] - (mm/2))/σ_a)
      end
    end #k
  end #j

  agrid = exp.(agrid) # Exponentiate grid

  return agrid, P #return productivity grid and Markov transition matrix
end

agrid, P = MyTauchen(σ,ρ,na,m)

#--------------------------------#
#    VALUE FUNCTION ITERATION
#--------------------------------#
# RECURSIVE PROBLEM:
#
# V(k,a) = max_{l,k_(+1)} log(c) - χ(l) + β*E[V(k_(+1),a_(t+1)]}
#
# subject to
#
# c + k - (1-δ)*k = a*k^α*l^(1-α

struct modelstate
  ind::Int64
  na::Int64
  nk::Int64
  nh::Int64
  P::Array
  agrid::Vector{Float64}
  kgrid::Vector{Float64}
  hgrid::Vector{Float64}
  β::Float64 # Discount rate
  η::Float64 # Relative risk aversion
  ψ::Float64 # inverse Frisch elasticity
  δ::Float64 # Depreication rate
  α::Float64 # Capital share of income
  V::Array
end

# computes value function given model state
function value(currentState::modelstate)
  ind     = currentState.ind
  na      = currentState.na
  nk      = currentState.nk
  nh      = currentState.nh
  P       = currentState.P
  agrid   = currentState.agrid
  kgrid   = currentState.kgrid
  hgrid   = currentState.hgrid
  β       = currentState.β
  η       = currentState.η
  ψ       = currentState.ψ
  δ       = currentState.δ
  α       = currentState.α
  VV      = currentState.V

  # Get states
  kk      = convert(Int, floor((ind-0.05)/na))+1
  aa      = convert(Int, floor(mod(ind-0.05,na))+1)

  # states
  kc      = kgrid[kk] # current capital
  ac      = agrid[aa] # current technology

  # hold capital kc constant, compute continuation value conditional on "ac" technology
  cont    = P[aa,:]' * VV
  cont    = vec(repeat(cont, nk)) # value is constant in the dimension of h-choice

  # compute model objects
  kkgr    = vec(repeat(kgrid',nh))
  hhgr    = vec(repeat(hgrid,nk,1))

  output      = ac.*kc.^α.*hhgr.^(1.0-α)
  investment  = kkgr .- (1.0-δ)*kc
  consumption = output - investment

  # compute new value funciton
  VVp = zeros(nh*nk)
  for (i,cons) in enumerate(consumption)
    if cons <= 0
      VVp[i] = -Inf
    else
      VVp[i] = cons^(1.0-η)/(1.0-η) - (hhgr[i])^(1.0+ψ)/(1.0+ψ) + β*cont[i]
    end
  end

  # find maximum
  Vnew, tempidx = findmax(vec(VVp))

  return Vnew, tempidx
end

#--------------------------------#
#           MAIN LOOP
#--------------------------------#

function VFI(Vbar::Float64, agrid::Array, kgrid::Array, hgrid::Array, P::Array; α = α, β = β, η = η, ψ = ψ, δ = δ, crit = Inf, maxiter = 1000, iter = 0, tol = 1.0e-8, nfix = 10)

  na = size(agrid,1)
  nh = size(hgrid,1)
  nk = size(kgrid,1)
  V_init = Vbar*ones(na,nk)
  V_next = Array{Float64,1}
  index = ones(na*nk,)

  # useful grids
  kagr = vec(repeat(kgrid',na))
  aagr = repeat(agrid,nk,1)
  hhgr = repeat(hgrid,nk,1)
  kkgr = vec(repeat(kgrid',nh))

  kidx    = collect(1:1:nh)
  kidx    = repeat(kidx,1,nh)'

  while (crit > tol && iter < maxiter)

    V_next = copy(V_init)

    if mod(iter,nfix) == 0 # speed-up: optimization every 10 iterations

      for ind = 1:(na*nk)
        State   = modelstate(ind,na,nk,nh,P,agrid,kgrid,hgrid,β,η,ψ,δ,α,V_next)
        tempVal, tempIndex  = value(State)

        V_next[ind] = tempVal
        index[ind] = tempIndex
      end

      else # on mod(iter,nfix) != 0, update the value function only.

        # Re-use conjectured index
        idx = convert(Array{Int64,1}, index)

        # E(V) condition on A(t) and each possible K(t+1)
        EVp = P*V_init

        # E(V) conditioned on A(t) and optimal K(t+1)
        ktmp = kidx[idx]
        evp_k = zeros(na*nk)
        for ind = 1:(na*nk)
          kk = convert(Int, floor((ind-0.05)/na))+1
          aa = convert(Int, floor(mod(ind-0.05,na))+1)
          evp_k[ind] = EVp[aa,ktmp[ind]]
        end

        #Recompute value using conjectured value function
        output = aagr.*kagr.^α.*hhgr[idx].^(1.0-α)
        investment = kkgr[idx] .- (1.0-δ)*kagr
        consumption = output - investment

        # Recompute value function under conjectured policy
        V_next = consumption.^(1.0-η)./(1.0-η) - hhgr[idx].^(1.0+ψ)./(1.0+ψ) + β*evp_k
        V_next = reshape(V_next,na,nk)
    end

    crit = maximum(abs.(V_next .- V_init))
    V_init = copy(V_next)
    iter += 1

    print("iteration: ",iter, ", critical val: ", crit*1.0e8," times 10^8", "\n")
  end #while

  return V_next, index, crit
end

@time Valuefunction, index, crit = VFI(Vbar,agrid,kgrid,hgrid,P)

#--------------------------------#
#    PLOT POLICY FUNCTIONS
#--------------------------------#

index   = convert(Array{Int64,1}, index)
pols    = reshape(index,(na,nk)) # policies

alow    = pols[cover-1,:] # fix productivity at low
amid    = pols[cover+1,:] # fix productivity at steady state
ahigh   = pols[cover+2,:] # fix productivity high

# Policies for low,mid,high realisations of "a"
hpol_low  = floor.(mod.(alow.-0.05,nh)) .+ 1.0
hpol_mid  = floor.(mod.(amid.-0.05,nh)) .+ 1.0
hpol_high = floor.(mod.(ahigh.-0.05,nh)) .+ 1.0
kpol_low  = floor.((alow.-0.05)./nk) .+ 1.0
kpol_mid  = floor.((amid.-0.05)./nk) .+ 1.0
kpol_high = floor.((ahigh.-0.05)./nk) .+ 1.0

# convert to integer arrays
hpol_low = convert(Array{Int64},hpol_low)
hpol_mid = convert(Array{Int64},hpol_mid)
hpol_high = convert(Array{Int64},hpol_high)
kpol_low = convert(Array{Int64},kpol_low)
kpol_mid = convert(Array{Int64},kpol_mid)
kpol_high = convert(Array{Int64},kpol_high)

# make plots
p1 = plot(kgrid,hgrid[hpol_low], ylims = (hmin,hmax), ylabel="low a", title="h policy")
p3 = plot(kgrid,hgrid[hpol_mid], ylims = (hmin,hmax), ylabel="steady state a")
p5 = plot(kgrid,hgrid[hpol_high], ylims = (hmin,hmax), ylabel="high a")

p2 = plot(kgrid,kgrid[kpol_low], ylims = (kmin,kmax), title="k policy")
p4 = plot(kgrid,kgrid[kpol_mid], ylims = (kmin,kmax))
p6 = plot(kgrid,kgrid[kpol_high], ylims = (kmin,kmax))

VFI_policyfunctions = plot(p1,p2,p3,p4,p5,p6, layout = (3,2), legend = false)
display(VFI_policyfunctions)
savefig(VFI_policyfunctions,"./VFI_policyfunc.png")


text:
# Julia wrapper for header: /usr/include/scip/pub_pricer.h
# Automatically generated using Clang.jl wrap_c

function SCIPpricerComp(elem1, elem2)
    ccall((:SCIPpricerComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpricerCompName(elem1, elem2)
    ccall((:SCIPpricerCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpricerGetData(pricer)
    ccall((:SCIPpricerGetData, libscip), Ptr{SCIP_PRICERDATA}, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerSetData(pricer, pricerdata)
    ccall((:SCIPpricerSetData, libscip), Cvoid, (Ptr{SCIP_PRICER}, Ptr{SCIP_PRICERDATA}), pricer, pricerdata)
end

function SCIPpricerGetName(pricer)
    ccall((:SCIPpricerGetName, libscip), Cstring, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetDesc(pricer)
    ccall((:SCIPpricerGetDesc, libscip), Cstring, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetPriority(pricer)
    ccall((:SCIPpricerGetPriority, libscip), Cint, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetNCalls(pricer)
    ccall((:SCIPpricerGetNCalls, libscip), Cint, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetNVarsFound(pricer)
    ccall((:SCIPpricerGetNVarsFound, libscip), Cint, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetSetupTime(pricer)
    ccall((:SCIPpricerGetSetupTime, libscip), Cdouble, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerGetTime(pricer)
    ccall((:SCIPpricerGetTime, libscip), Cdouble, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerIsActive(pricer)
    ccall((:SCIPpricerIsActive, libscip), UInt32, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerIsDelayed(pricer)
    ccall((:SCIPpricerIsDelayed, libscip), UInt32, (Ptr{SCIP_PRICER},), pricer)
end

function SCIPpricerIsInitialized(pricer)
    ccall((:SCIPpricerIsInitialized, libscip), UInt32, (Ptr{SCIP_PRICER},), pricer)
end