text:
module ParametricCurves

using ForwardDiff
using StaticArrays
using LinearAlgebra

abstract type AbstractCurve{N} end

abstract type AbstractParametricCurve{N} <: AbstractCurve{N} end

function derivatives(curve::AbstractParametricCurve{N}, n::Int) where {N}
    # default to AD
    return derivatives(curve, Val(n))
end

derivatives(curve::AbstractParametricCurve{N}) where {N} = derivatives(curve, 1)

function derivatives(curve::AbstractParametricCurve{N}, ::Val{n}) where {N,n}
    # concrete implementation of AbstractParametricCurve can define this function
    # so that AD wont be used.
    # this can be done by defining your own ParametricCurve type.
    return _derivatives(curve.coordinates, Val(n))
end
Base.getproperty(curve::AbstractParametricCurve, s::Symbol) = _getproperty(curve, Val(s))
_getproperty(curve::AbstractParametricCurve, ::Val{S}) where {S} = getfield(curve, S)
_getproperty(curve::AbstractParametricCurve, ::Val{:x}) = curve.coordinates[1]
_getproperty(curve::AbstractParametricCurve, ::Val{:y}) = curve.coordinates[2]
_getproperty(curve::AbstractParametricCurve, ::Val{:z}) = curve.coordinates[3]
_getproperty(curve::AbstractParametricCurve{2}, ::Val{:z}) = error("2D curve has no z coordinate.")

function _derivatives(fs, ::Val{1})
    return map(f -> (t -> ForwardDiff.derivative(f, t)), fs)
end
function _derivatives(fs, ::Val{n}) where {n}
    return map(f -> (t -> ForwardDiff.derivative(f, t)), _derivatives(fs, Val(n - 1)))
end
# function _derivatives(fs, n::Int)
#     return _derivatives(fs, Val(n))
# end

@inline function tangent(curve::AbstractParametricCurve{N}, t::T) where {N,T}
    # compute the tangent vector at the parameter location t
    dval = map((f) -> f(t), derivatives(curve))
    #convert tuple to SVector
    dval = SVector{N,T}(dval)
    return dval / norm(dval)
end

@inline function normal(curve::AbstractParametricCurve{2}, t::T) where {T}
    # in2D, can rotate tangent to get normal
    tvec = tangent(curve, t)
    return SVector{2,T}(-tvec[2], tvec[1])
end

@inline function normal(curve::AbstractParametricCurve{N}, t::T) where {N,T}
    # compute the tangent vector at the parameter location t
    dval = map((f) -> f(t), derivatives(curve, Val(2)))
    #convert tuple to SVector
    dval = SVector{N,T}(dval)
    return dval / norm(dval)
end

@generated function radius(curve::AbstractParametricCurve{N}, t::Real) where {N}
    ex = :(curve.coordinates[1](t)^2)
    for i = 2:N
        ex = :(curve.coordinates[$i](t)^2 + $ex)
    end
    return :(sqrt($ex))
end

@inline function cartesian(curve::AbstractParametricCurve{N}, t::T) where {N,T}
    x = map(f -> f(t), curve.coordinates)
    return SVector{N,T}(x)
end

"""
    explicit2parametric(f)
Convert an explicit 2D curve of the form `y=f(x)` to its parametric form.
"""
@inline function explicit2parametric(f)
    return (t -> t, t -> f(t))
end

struct ParametricCurve{N,C} <: AbstractParametricCurve{N}
    coordinates::C
    function ParametricCurve(coordinates)
        N = length(coordinates)
        return new{N,typeof(coordinates)}(coordinates)
    end
end

export ParametricCurve, tangent, normal, radius, cartesian, explicit2parametric

end # end module


text:
using GainPatterns
using Base.Test

#################################################################
# testing angular_error function
#################################################################
@test angular_error(5, 355) == 10
@test angular_error(355, 5) == 10

@test angular_error(5, -5) == 10
@test angular_error(-5, 5) == 10

@test angular_error(340, 17) == 37
@test angular_error(17, 340) == 37

@test angular_error(-17, 340) == 3
@test angular_error(340, -17) == 3

#################################################################
# testing angular_error_rel function
#################################################################
@test angular_error_rel(5, 355) == -10
@test angular_error_rel(355, 5) == 10

@test angular_error_rel(5, -5) == -10
@test angular_error_rel(-5, 5) == 10

@test angular_error_rel(340, 17) == 37
@test angular_error_rel(17, 340) == -37

@test angular_error_rel(-17, 340) == -3
@test angular_error_rel(340, -17) == 3

v1 = [5, 355, 5, -5, 340, 17, -17, 340]
v2 = [355, 5, -5, 5, 17, 340, 340, -17]
vec_ans = [-10, 10, -10, 10, 37, -37, -3, 3]
@test angular_error_rel(v1,v2) == vec_ans

ext: jl


text:
int while(){ return 1; }

int main() { return 1; }

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using Gadfly, DataFrames

#demonstration of BLE channel hopping noise avoidance
function gaussian(x, mu, sigmaSquared)
    sigma = sqrt(sigmaSquared)
    1 / (sigma * sqrt(2 * pi)) * exp(-1.0 * ((x - mu) ^ 2) / (2.0 * (sigmaSquared) ))
end

WiFiFreqs = [2412, 2437, 2462]
BLEFreqs = vcat([2404 : 2 : 2424], [2428 : 2 : 2478])
x = [2400 : 1 : 2500]

WiFiPlots = [2.0 * gaussian(ix, iy, 40.0) for ix=x, iy=WiFiFreqs]

BLEValues = [ ((findfirst(BLEFreqs, ix) > 0) ? 0.127 : 0/0) for ix=x]
BLEChannels = [findfirst(BLEFreqs, ix) for ix=x]
BLEChannelStrs = [(ix < 10)?"0$ix":"$ix" for ix=BLEChannels]
df = DataFrame(frequency=x, value=BLEValues, channel=BLEChannels)

plot(
    layer(x=x, y=WiFiPlots[:, 1], Geom.line, Theme(default_color=color("#95a5b5"), line_width=7px)),
    layer(x=x, y=WiFiPlots[:, 2], Geom.line, Theme(default_color=color("#95a5b5"), line_width=7px)),
    layer(x=x, y=WiFiPlots[:, 3], Geom.line, Theme(default_color=color("#95a5b5"), line_width=7px)),
    layer(df, x="frequency", y="value", color="channel", Geom.bar,
          Theme(default_color=color("#ff8585"), default_point_size=5px)),
     Scale.y_sqrt
)


text:
###############################################################################
#
# Initialize Canopy4RT
#
###############################################################################
"""
    create_canopy_rt(FT; nLayer::Int = 20, LAI::Number = FT(3))

Create [`Canopy4RT`](@ref), given
- `FT` Floating number type
- `nLayer` Number of canopy layers
- `LAI` Leaf area index
"""
function create_canopy_rt(FT; nLayer::Int = 20, LAI::Number = FT(3))
    return Canopy4RT{FT}(LAI    = LAI   ,
                         nLayer = nLayer)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# abstract type AbstractRule{T} end

# abstract type ConditionalRule{T} <: AbstractRule{T} end
# struct SimpleRule{T} <: AbstractReplacementRule{T}
#     r::T
#     priority::Int
# end

# struct UpwindRule{T, T2}
#     r::T
#     priority::Int
#     condition::T2
# end
# struct RuleSet{T, T2}
#     derivrules::T
#     valrules::T2
# end

# function RuleSet(rules::Vector{T}, conditional_rules::Vector{C}) where {T<:SimpleRule, C<:ConditionalRule}
#     priorities = vcat(map(r -> r.priority, rules), map(r -> r.priority, conditional_rules))

#     for (i,r) in enumerate(vcat(rules, conditional_rules))

#     end

# ModelingToolkit.substitute(expr, rule::AbstractRule{T}) where T = substitute(expr, rule.r)

"""
`interpolate_discrete_param`

Interpolate gridpoints by taking the average of the values of the discrete points, or if the offset is outside the grid, extrapolate the value with dx.
"""
@inline function interpolate_discrete_param(i, s, itap, x, bpc)

    return s.grid[x][i+itap]+s.dxs[x]*.5

end

"""
`cartesian_nonlinear_laplacian`

Differential(x)(expr(x)*Differential(x)(u(x)))

Given an internal multiplying expression `expr`, return the correct finite difference equation for the nonlinear laplacian at the location in the grid given by `II`.

The inner derivative is discretized with the half offset centered scheme, giving the derivative at interpolated grid points offset by dx/2 from the regular grid.

The outer derivative is discretized with the centered scheme, giving the nonlinear laplacian at the grid point `II`.
For first order returns something like this:
`d/dx( a du/dx ) ~ (a(x+1/2) * (u[i+1] - u[i]) - a(x-1/2) * (u[i] - u[i-1]) / dx^2`

For 4th order, returns something like this:
```
first_finite_diffs = [a(x-3/2)*finitediff(u, i-3/2),
                      a(x-1/2)*finitediff(u, i-1/2),
                      a(x+1/2)*finitediff(u, i+1/2),
                      a(x+3/2)*finitediff(u, i+3/2)]

dot(central_finite_diff_weights, first_finite_diffs)
```

where `finitediff(u, i)` is the finite difference at the interpolated point `i` in the grid.

And so on.
"""
function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace{N}, b, depvars, x, u) where N
    # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
    # See scheme 1, namely the term without the 1/r dependence. See also #354 and #371 in DiffEqOperators, the previous home of this package.
    ndims(u,s) == 0 && return Num(0)
    jx = j, x = (x2i(s, u, x), x)
    @assert II[j] != 1 "The nonlinear laplacian is only defined on the interior of the grid, it is unsupported in boundary conditions."
    @assert II[j] != length(s, x) "The nonlinear laplacian is only defined on the interior of the grid, it is unsupported in boundary conditions."

    D_inner = derivweights.halfoffsetmap[Differential(x)]
    inner_interpolater = derivweights.interpmap[x]

    # Get the outer weights and stencil. clip() essentially removes a point from either end of the grid, for this reason this function is only defined on the interior, not in bcs#
    cliplen = length(s, x) - 1

    outerweights, outerstencil = get_half_offset_weights_and_stencil(D_inner, II-unitindex(N,j), s, b, u, jx, cliplen)

    # Get the correct weights and stencils for this II
    inner_deriv_weights_and_stencil = [get_half_offset_weights_and_stencil(D_inner, I, s, b, u, jx) for I in outerstencil]
    interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, I, s, b, u, jx) for I in outerstencil]

    # map variables to symbolically inerpolated/extrapolated expressions
    map_vars_to_interpolated(stencil, weights) = [v => dot(weights, s.discvars[v][stencil]) for v in depvars]

    # Map parameters to interpolated values. Using simplistic extrapolation/interpolation for now as grids are uniform
    #TODO: make this more efficient
    map_params_to_interpolated(stencil, weights) = vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(stencil, (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])])

    # Take the inner finite difference
    inner_difference = [dot(inner_weights, s.discvars[u][inner_stencil]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil]

    # Symbolically interpolate the multiplying expression

    interpolated_expr = map(interp_weights_and_stencil) do (weights, stencil)
        Num(substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_params_to_interpolated(stencil, weights)))
    end

    # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference
    return dot(outerweights, inner_difference .* interpolated_expr)
end

"""
`spherical_diffusion`

Based on https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf

See scheme 1 in appendix A. The r = 0 case is treated in a later appendix
"""
function spherical_diffusion(innerexpr, II, derivweights, s, b, depvars, r, u)
    # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
    D_1 = derivweights.map[Differential(r)]
    D_2 = derivweights.map[Differential(r)^2]

    #TODO!: Update this to use indvars of the pde
    # What to replace parameter x with given I
    _rsubs(x, I) = x => s.grid[x][I[s.x2i[x]]]
    # Full rules for substituting parameters in the inner expression
    rsubs(I) = vcat([v => s.discvars[v][I] for v in depvars], [_rsubs(x, I) for x in s.x̄])
    # Discretization func for u
    ufunc_u(v, I, x) = s.discvars[v][I]

    # 2nd order finite difference in u
    exprhere = Num(substitute(innerexpr, rsubs(II)))
    # Catch the r ≈ 0 case
    if Symbolics.unwrap(substitute(r, _rsubs(r, II))) ≈ 0
        D_2_u = central_difference(D_2, II, s, b, (s.x2i[r], r), u, ufunc_u)
        return 3exprhere*D_2_u # See appendix B of the paper
    end
    D_1_u = central_difference(D_1, II, s, b, (s.x2i[r], r), u, ufunc_u)
    # See scheme 1 in appendix A of the paper

    return exprhere*(D_1_u/Num(substitute(r, _rsubs(r, II))) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, b, depvars, r, u))
end

@inline function generate_cartesian_rules(II, s, depvars, derivweights, pmap, indexmap, terms)
    central_ufunc(u, I, x) = s.discvars[u][I]
    return reduce(vcat, [reduce(vcat, [[(Differential(x)^d)(u) => central_difference(derivweights.map[Differential(x)^d], Idx(II, s, u, indexmap), s, pmap.map[operation(u)][x], (x2i(s,u,x),x), u, central_ufunc) for d in (let orders = derivweights.orders[x]; orders[iseven.(orders)] end)] for x in params(u, s)]) for u in depvars])
end

@inline function upwind_difference(expr, d::Int, II::CartesianIndex{N}, s::DiscreteSpace{N}, b, depvars, derivweights, (j,x), u, central_ufunc, indexmap) where N
    # TODO: Allow derivatives in expr
    expr = substitute(expr, valmaps(s, u, depvars, Idx(II, s, depvar(u, s), indexmap), indexmap))
    IfElse.ifelse(expr > 0,
                  expr*upwind_difference(d, II, s, b, derivweights, (j,x), u, central_ufunc, true),
                  expr*upwind_difference(d, II, s, b, derivweights, (j,x), u, central_ufunc, false))
end

@inline function generate_winding_rules(II, s, depvars, derivweights, pmap, indexmap, terms)
    wind_ufunc(v, I, x) = s.discvars[v][I]
    # for all independent variables and dependant variables
    rules = vcat(#Catch multiplication
                reduce(vcat, [reduce(vcat, [[@rule *(~~a, $(Differential(x)^d)(u), ~~b) => upwind_difference(*(~a..., ~b...), d, Idx(II, s, u, indexmap), s, pmap.map[operation(u)][x], depvars, derivweights, (x2i(s,u,x),x), u, wind_ufunc, indexmap) for d in (let orders = derivweights.orders[x]; orders[isodd.(orders)] end)] for x in params(u, s)]) for u in depvars]),

                #Catch division and multiplication, see issue #1
                reduce(vcat, [reduce(vcat, [[@rule /(*(~~a, $(Differential(x)^d)(u), ~~b), ~c) => upwind_difference(*(~a..., ~b...)/~c, d, Idx(II, s, u, indexmap), s, pmap.map[operation(u)][x], depvars, derivweights, (x2i(s,u,x), x), u, wind_ufunc, indexmap) for d in (let orders = derivweights.orders[x]; orders[isodd.(orders)] end)] for x in params(u, s)]) for u in depvars])
                )

    wind_rules = []

    # wind_exprs = []
    for t in terms
        for r in rules
            if r(t) !== nothing
                push!(wind_rules, t => r(t))
            end
        end
    end

    return vcat(wind_rules, vec(mapreduce(vcat, depvars) do u
        mapreduce(vcat, params(u, s)) do x
            j = x2i(s,u,x)
            let orders = derivweights.orders[x]
                oddorders = orders[isodd.(orders)]
                # for all odd orders
                if length(oddorders) > 0
                    map(oddorders) do d
                        (Differential(x)^d)(u) => upwind_difference(d, Idx(II, s, u, indexmap), s, pmap.map[operation(u)][x], derivweights, (j,x), u, wind_ufunc, true)
                    end
                else
                    []
                end
            end
        end
    end))

end

@inline function generate_nonlinlap_rules(II, s, depvars, derivweights, pmap, indexmap, terms)
    rules = reduce(vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c...,cartesian_nonlinear_laplacian(*(a..., b...), Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][x], depvars, x, u), ~d...) for x in params(u, s)]) for u in depvars])

    rules = vcat(rules, reduce(vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(a..., b...), Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][x], depvars, x, u) for x in params(u, s)]) for u in depvars]))

    rules = vcat(rules, reduce(vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u)/~a)) => cartesian_nonlinear_laplacian(1/~a, Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][x], depvars, x, u) for x in params(u, s)]) for u in depvars]))

    nonlinlap_rules = []
    for t in terms
        for r in rules
            if r(t) !== nothing
                push!(nonlinlap_rules, t => r(t))
            end
        end
    end
    return nonlinlap_rules
end

@inline function generate_spherical_diffusion_rules(II, s, depvars, derivweights, pmap, indexmap, terms)
    rules = reduce(vcat, [vec([@rule *(~~a, 1/(r^2), ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) => *(~a..., spherical_diffusion(*(~c..., ~d..., ~e...), Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][r], depvars, r, u), ~b...)
            for r in params(u,s)]) for u in depvars])

    rules = vcat(rules, reduce(vcat, [vec([@rule /(*(~~a, $(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e)), ~~b), (r^2)) => *(~a..., ~b..., spherical_diffusion(*(~c..., ~d..., ~e...), Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][r], depvars, r, u))
    for r in params(u,s)]) for u in depvars]))

    rules = vcat(rules, reduce(vcat, [vec([@rule /(($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), (r^2)) => spherical_diffusion(*(~c..., ~d..., ~e...), Idx(II, s, u, indexmap), derivweights, s, pmap.map[operation(u)][r], depvars, r, u)
    for r in params(u, s)]) for u in depvars]))

    spherical_diffusion_rules = []
    for t in terms
        for r in rules
            if r(t) !== nothing
                push!(spherical_diffusion_rules, t => r(t))
            end
        end
    end
    return spherical_diffusion_rules
end

"""
`generate_finite_difference_rules`

Generate a vector of finite difference rules to dictate what to replace variables in the `pde` with at the gridpoint `II`.

Care is taken to make sure that the rules only use points that are actually in the discretized grid by progressively up/downwinding the stencils when the gridpoint `II` is close to the boundary.

There is a genral catch all ruleset that uses the cartesian centered difference scheme for derivatives, and simply the discretized variable at the given gridpoint for particular variables.

There are of course more specific schemes that are used to improve stability/speed/accuracy when particular forms are encountered in the PDE. These rules are applied first to override the general ruleset.

##Currently implemented special cases are as follows:
    - Spherical derivatives
    - Nonlinear laplacian uses a half offset centered scheme for the inner derivative to improve stability
    - Spherical nonlinear laplacian.
    - Up/Downwind schemes to be used for odd ordered derivatives multiplied by a coefficient, downwinding when the coefficient is positive, and upwinding when the coefficient is negative.

Please submit an issue if you know of any special cases which impact stability or accuracy that are not implemented, with links to papers and/or code that demonstrates the special case.
"""
function generate_finite_difference_rules(II, s, depvars, pde, derivweights, pmap, indexmap)

    terms = split_terms(pde, s.x̄)

    # Standard cartesian centered difference scheme
    central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, pmap, indexmap, terms)

    # Nonlinear laplacian scheme
    nonlinlap_rules = generate_nonlinlap_rules(II, s, depvars, derivweights, pmap, indexmap, terms)

    # Because winding needs to know about multiplying terms, we can't split the terms into additive and multiplicative terms.
    winding_rules = generate_winding_rules(II, s, depvars, derivweights, pmap, indexmap, terms)

    # Spherical diffusion scheme
    spherical_diffusion_rules = generate_spherical_diffusion_rules(II, s, depvars, derivweights, pmap, indexmap, split_additive_terms(pde))

    return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(winding_rules), vec(central_deriv_rules_cartesian))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV, Statistics, Dates, Plots, Measures
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.19"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("temperatures.csv"))
	brisbane = data.Brisbane
	goldcoast = data.GoldCoast

	diff = brisbane - goldcoast
	dates = [Date(
				Year(data.Year[i]),
				Month(data.Month[i]),
				Day(data.Day[i])
			) for i in 1:nrow(data)]

	fortnightRange = 250:263
	brisFortnight = brisbane[fortnightRange]
	goldFortnight = goldcoast[fortnightRange]
end;

# ╔═╡ 47673b4e-250b-11eb-16f2-2f9e1f103b04
describe(data)

# ╔═╡ 22fd2570-250b-11eb-02af-73597d381ffb
begin
	default(xlabel="Time", ylabel="Temperature")
	default(label=["Brisbane" "Gold Coast"])

	p1 = plot(dates, [brisbane goldcoast],
			c=[:blue :red])
	p2 = plot(dates[fortnightRange], [brisFortnight goldFortnight],
			c=[:blue :red], m=(:dot, 5, Plots.stroke(1)))
	p3 = plot(dates, diff,
			c=:black, ylabel="Temperature Difference",legend=false)
	p4 = histogram(diff, bins=-4:0.5:6,
			ylims=(0,140), legend = false,
			xlabel="Temperature Difference", ylabel="Frequency")
	plot(p1,p2,p3,p4, size = (800,500), margin = 5mm)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.19"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═47673b4e-250b-11eb-16f2-2f9e1f103b04
# ╠═22fd2570-250b-11eb-02af-73597d381ffb
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
# Board struct
struct Board
  n :: Int
  lines :: Vector{Line}
  cells :: Matrix{Int}
end

# Creates board that can hold the given lines.
function Board(lines::Vector{Line})
  n = maximum(map(l -> max(l.startPt[1],l.startPt[2],l.endPt[1],l.endPt[2]),lines))
  return Board(n,lines,zeros(Int,(n,n)))
end

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

# Line struct
struct Line
  startPt :: CartesianIndex{2}
  endPt   :: CartesianIndex{2}
  points  :: Vector{CartesianIndex{2}}
end

# Creates line from the two endpoints.
function Line(p1::CartesianIndex{2}, p2::CartesianIndex{2})
  slope = (p2-p1)
  npts = max(abs(slope[1]), abs(slope[2]))
  slope = CartesianIndex(slope[1] ÷ npts,slope[2] ÷ npts)
  pts = [p1 + i*slope for i = 0:npts]
  return Line(p1, p2, pts)
end

# Check line properties
isHorizontal(l::Line) :: Bool = (l.startPt[1] == l.endPt[1])
isVertical(l::Line)   :: Bool = (l.startPt[2] == l.endPt[2])
isFlat(l::Line)       :: Bool = isHorizontal(l) || isVertical(l)

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

function parsePoint(s::String) :: CartesianIndex{2}
  coords = map(x -> parse(Int,x),split(s, ',', keepempty=false))
  return CartesianIndex(coords[1], coords[2])
end

function parseLine(s::String) :: Line
  pts = split(s, "->", keepempty=false)
  p1 = parsePoint(String(pts[1]))
  p2 = parsePoint(String(pts[2]))
  return Line(p1,p2)
end

function readInput() :: Vector{Line}
  f = open("./inputs/day5.txt")

  lines :: Vector{Line} = []
  for line in eachline(f)
    push!(lines, parseLine(line))
  end
  return lines
end

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

# Draw a single line into a matrix
function add!(line::Line, M::Matrix{Int})
  for p in line.points
    M[p] += 1
  end
end

# Draw all lines in a board
function drawLines!(b::Board)
  map(l -> add!(l,b.cells), b.lines)
end

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

lines = readInput()

# Problem 1
lines1 = deepcopy(filter(isFlat, lines)) # Only consider flat lines
board1 = Board(lines1)
drawLines!(board1)
res = count( map( x -> (x > 1) ? true : false, board1.cells))

# Problem 2
lines2 = deepcopy(lines)
board2 = Board(lines2)
drawLines!(board2)
res = count( map( x -> (x > 1) ? true : false, board2.cells))


text:
module Model

using Gen
using ..MangroveBase: ReviewInfo, Rating, RATINGS

@gen function subset(items::Vector{String})::Vector
  filter(item -> @trace(bernoulli(0.5), item), items)
end

@dist quality() = uniform_discrete(1, RATINGS)

@gen function triangle_categorical(peak::Real)::Float64
  raw = [(5 - abs(x - peak)) for x in 1:5]
  @trace(categorical(raw ./ sum(raw)), :rating)
end

@gen function mangrove_model(subjects::Set{String}, reviewers::Set{String})
  for sub in subjects
    true_quality = @trace(quality(), :quality => sub)
    for kid in @trace(subset(collect(reviewers)), :selected => sub)
      @trace(
        triangle_categorical(true_quality),
        ReviewInfo(sub, kid)
      )
    end
  end
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
"""
    shortest_path(g::OSMGraph,
                  origin,
                  destination;
                  weights=g.weights,
                  restrictions=g.indexed_restrictions,
                  algorithm::Symbol=:dijkstra,
                  save_dijkstra_state::Bool=false,
                  heuristic::Union{Function,Nothing}=nothing
                  )::Vector{U} where {U <: Integer}

Calculates the shortest path between two OpenStreetMap node ids.

# Arguments
- `g::OSMGraph`: Graph container.
- `origin`: Origin OpenStreetMap node or node id.
- `destination`: Destination OpenStreetMap node or node id.
- `weights`: Matrix of node to node edge weights, defaults to `g.weights`. If a custom weights matrix
    is being used with algorithm set to `:astar`, make sure that a correct heuristic is being used.
- `restrictions `: Turn restictions dictionary. Pass in `nothing` to not use
    restrictions, defaults to `g.indexed_restrictions`.
- `algorithm::Symbol=:dijkstra`: Shortest path algorithm, either. `:dijkstra` or `:astar`.
- `save_dijkstra_state::Bool=false`: Option to cache dijkstra parent states from a single source.
- `heuristic::Union{Function,Nothing}=nothing`: Option to use custom astar heuristic. If `nothing`,
    default haversine distance will be used if `g.weight_type` is `:distance`, otherwise an estimated
    time will be used by dividing the harversine distance by the speed at this node (specified in the
    `maxspeed` tag).

# Return
- `Vector{U}`: Array of OpenStreetMap node ids making up the shortest path.
"""
function shortest_path(g::OSMGraph,
                       origin::U,
                       destination::U;
                       weights=g.weights,
                       restrictions=g.indexed_restrictions,
                       algorithm::Symbol=:dijkstra,
                       save_dijkstra_state::Bool=false,
                       heuristic::Union{Function,Nothing}=nothing
                       )::Vector{U} where {U <: Integer}

    o_index = g.node_to_index[origin]
    d_index = g.node_to_index[destination]

    if isassigned(g.dijkstra_states, o_index)
        # State already exists
        parents = g.dijkstra_states[o_index]

    elseif algorithm == :dijkstra
        if save_dijkstra_state
            parents = dijkstra(g.graph, o_index, distmx=weights, restrictions=restrictions)
            g.dijkstra_states[o_index] = parents
        else
            parents = dijkstra(g.graph, o_index; goal=d_index, distmx=weights, restrictions=restrictions)
        end

    elseif algorithm == :astar
        h = heuristic === nothing ? default_heuristic(g) : heuristic
        parents = astar(g.graph, o_index, goal=d_index, distmx=weights, heuristic=h, restrictions=restrictions)

    else
        throw(ErrorException("No such algorthm $algorithm, pick from `:dijkstra` or `:astar`"))
    end

    path = path_from_parents(parents, o_index, d_index)

    return [g.index_to_node[i] for i in path]
end

function shortest_path(g::OSMGraph,
                       origin::Node{U},
                       destination::Node{U};
                       kwargs...)::Vector{U} where {U <: Integer}
    return shortest_path(g, origin.id, destination.id; kwargs...)
end

"""
Returns the default heuristic function used in astar shorest path calculation.
"""
function default_heuristic(g::OSMGraph, distance_function::Symbol=:haversine)::Function
    if g.weight_type == :distance
        return (u, v) -> distance(g.node_coordinates[u], g.node_coordinates[v], distance_function)
    elseif g.weight_type == :time || g.weight_type == :lane_efficiency
        return (u, v) -> distance(g.node_coordinates[u], g.node_coordinates[v], distance_function) / g.nodes[g.index_to_node[u]].tags["maxspeed"]
    end
end

"""
    path_from_parents(parents::Vector{U}, origin::V, destination::V)::Vector where {U <: Integer, V <: Integer}

Extracts array of shortest path node ids from a single source dijkstra parent state.

# Arguments
- `parents::Vector{U}`: Array of dijkstra parent states.
- `origin::U`: Origin OpenStreetMap node id.
- `destination::U`: Destination OpenStreetMap node id.

# Return
- `Vector{U}`: Array of OpenStreetMap node ids making up the shortest path.
"""
function path_from_parents(parents::Vector{U}, origin::V, destination::V)::Vector where {U <: Integer,V <: Integer}
    if parents[destination] == 0
        throw(ErrorException("Path does not exist between origin node index $origin and destination node index $destination"))
    end

    src = destination
    path = []

    while src != 0 # parent of origin is always 0
        push!(path, src)
        src = parents[src]
    end

    return reverse(path)
end

"""
    weights_from_path(g::OSMGraph{U,T,W}, path::Vector{T}; weights=g.weights)::Vector{W} where {U <: Integer,T <: Integer,W <: Real}

Extracts edge weights from a path using the weight matrix stored in `g.weights` unless
a different matrix is passed to the `weights` kwarg.

# Arguments
- `g::OSMGraph`: Graph container.
- `path::Vector{T}`: Array of OpenStreetMap node ids.
- `weights=g.weights`: the matrix that the edge weights are extracted from. Defaults to `g.weights`.

# Return
- `Vector{W}`: Array of edge weights, distances are in km, time is in hours.
"""
function weights_from_path(g::OSMGraph{U,T,W}, path::Vector{T}; weights=g.weights)::Vector{W} where {U <: Integer,T <: Integer,W <: Real}
    return [weights[g.node_to_index[path[i]], g.node_to_index[path[i + 1]]] for i in 1:length(path) - 1]
end

"""
    total_path_weight(g::OSMGraph{U,T,W}, path::Vector{T}; weights=g.weights)::W where {U <: Integer,T <: Integer,W <: Real}

Extract total edge weight along a path.

# Arguments
- `g::OSMGraph`: Graph container.
- `path::Vector{T}`: Array of OpenStreetMap node ids.
- `weights=g.weights`: the matrix that the edge weights are extracted from. Defaults to `g.weights`.

# Return
- `sum::W`: Total path edge weight, distances are in km, time is in hours.
"""
function total_path_weight(g::OSMGraph{U,T,W}, path::Vector{T}; weights=g.weights)::W where {U <: Integer,T <: Integer,W <: Real}
    sum::W = zero(W)
    for i in 1:length(path) - 1
        sum += weights[g.node_to_index[path[i]], g.node_to_index[path[i + 1]]]
    end
    return sum
end


text:
module LlvmMca

export code_mca, @code_mca

using InteractiveUtils, Base.Meta
import InteractiveUtils: gen_call_with_extracted_types_and_kwargs

macro code_mca(expr...)
    gen_call_with_extracted_types_and_kwargs(__module__, :code_mca, expr)
end

#= TODO

It seems to inline or include some base functions, like log2. I would
like to strip them out before passing the assembly to mca.

See also: https://github.com/JuliaLang/julia/issues/28046

=#

function code_mca(func, argtypes; save::Bool=false)
    io = IOBuffer()
    code_native(io, func, argtypes)

    s = String(take!(io))
    s = replace(s, r"^;.*\n"m => "")
    tmp, tmpio = mktemp()
    if save
        println(tmpio, s)
        println(stderr, "Saved temporary output to $tmp")
    end

    open(`llvm-mca`, "w", stdout) do cin
        println(cin, s)
    end
end

"""    @region_mca <code-block>

tries to insert llvm-mca code region markers.
"""
macro region_mca(expr)
    quote
        ccall(:fesetround, Int32, (Int32,), 0)
        t1 = Base.llvmcall("""
call void asm sideeffect "pause", "~{memory}"()
ret void
    """, Cvoid, Tuple{})
        $(esc(expr))
        ccall(:fesetround, Int32, (Int32,), 0)
        t2 = Base.llvmcall("""
call void asm sideeffect "pause", "~{memory}"()
ret void
    """, Cvoid, Tuple{})
    end
end

end


text:
include("interact.jl")

##

run_main     = true

display_loop = true

##

hidden_size   = 32

batch_size    = 50
learning_rate = 1
hm_epochs     = 1

generate_converge = true

##

rbm = nothing

main() = begin

    global rbm

    rbm = RBM(in_size, hidden_size)

    train(rbm=rbm,batch_size=batch_size,learning_rate=learning_rate,hm_epochs=hm_epochs)
    #train2(hidden_size=32,batch_size=10,learning_rate=1,hm_epochs=1)

    # generate(model)

    while display_loop

        display(generate(rbm,converge=generate_converge))

        sleep(1)

    end

    results = []

    # for _ in 1:10
    #
    #     model = RBM(in_size, 1)
    #
    #     update_weights!(model, batch_grads(model, choices(data_train,1)), 1)
    #
    #     push!(results, generate(model)) # ; display(results[end]) ; display(generate(model))
    #
    # end

results
end

##

test_basic() = begin

    rbm = RBM(10, 12)

    datapoint = binarize_data.(rand(1,10))

    inp = datapoint

    prev_visibles = nothing
    prev_hiddens = nothing

    for i in 1:5

        rbm(inp)
        rbm()

        inp = rbm.visibles

        println("iteration $(i)")
        @show rbm.visibles
        @show rbm.hiddens
        @show energy(rbm)

        prev_visibles == rbm.visibles ? println("visibles stable.") : prev_visibles = rbm.visibles
        prev_hiddens == rbm.hiddens ? println("hiddens stable.") : prev_hiddens = rbm.hiddens

    end

end ; #test_basic()

##

test_bestsize() = begin

    hm_avg = 100

    hidden_sizes = 1:12

    results = [0.0 for _ in hidden_sizes]

    for ii in 1:hm_avg

        for (i,hs) in enum(hidden_sizes)

            results[i] = (results[i] * (ii-1) + train(hidden_size=hs,hm_epochs=1,do_print=false)[end][end][end]) / ii

        end

        println("iteration: $(ii), results:")
        for (hs, res) in zip(hidden_sizes, results)
            println("\thidden_size: $(hs), loss: $(res)")
        end
        println("minimum loss: $(hidden_sizes[argmin(results)]) - $(argmin(results))")

    end

    println("Final Results:")
    for (hs,res) in zip(hidden_sizes, results)
        println("\t$(hs): $(res)")
    end

end ; #test_bestsize()

##

; run_main ? results = main() : ()


text:
using Documenter
using DocumenterTools
using MINLPLib

makedocs(
    format = Documenter.HTML(analytics = "UA-367975-10", mathengine = Documenter.MathJax()),
    sitename = "MINLPLib",
    pages = [
        "Libraries Summary" => "index.md",
        "Libraries" => [
            "bcp" => "bcp.md",
            "global" => "global.md",
            "ibm" => "ibm.md",
            "inf" => "inf.md",
            "minlp" => "minlp.md",
            "minlp2" => "minlp2.md",
            "morg" => "morg.md",
            "mpec" => "mpec.md",
            "mult3" => "mult3.md",
            "mult4" => "mult4.md",
            "poly" => "poly.md",
            "prince" => "prince.md",
            "qcqp" => "qcqp.md",
            "qcqp2" => "qcqp2.md",
            "qcqp3" => "qcqp3.md",
            "PODLib" => "PODLib.md"
            ],
        "Issues" => "issues.md"
    ]
)

deploydocs(
    repo = "github.com/lanl-ansi/MINLPLib.jl.git",
)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# smoke test
let
  @fixture function f() :f_result end
  @fixture function g(f) f end

  @pytest function(f)
    @test f == :f_result
  end
  @pytest function(g)
    @test g == :f_result
  end
end

# correct syntax, chaining of fixtures, basic usage flow
let
  @fixture function f() return [:f_result] end
  @fixture function g(f)
    [[:g_result]; f]
  end
  @fixture function h(g, f) return [[:h_result]; f; g] end
  @fixture function k(h, f) return [[:k_result]; h; f] end

  @pytest function(f, g, h, k)
    @test f == [:f_result]
    @test g == [:g_result; :f_result]
    @test h == [:h_result; :f_result; :g_result; :f_result]
    @test k == [:k_result; :h_result; :f_result; :g_result; :f_result; :f_result]
  end
end

# test-level (default) fixtures executed once per test
let
  @fixture function f() rand() end
  @fixture function g(f) f end
  @fixture function h(f, g) return [f; g] end

  the_fs = []

  @pytest function(f)
    push!(the_fs, f)
  end

  @pytest function(f, g, h)
    push!(the_fs, f)
    @test f == g
    @test h == [f, f]
  end

  @pytest function(h)
    push!(the_fs, h[1])
    @test h[1] == h[2]
  end

  @test allunique(the_fs) || "$the_fs should be different"
end

# no-fixture tests, execution of body anyway
let
  called = false
  @pytest function() called = true end
  @test called
end

# correct scoping within and outside fixture/test bodies
let
  @fixture function f() :f_result end
  @fixture function g(f)
    @test f == :f_result
    f
  end
  @fixture function h(g)
    @test typeof(f) == PyTest.Fixture
    @test g == :f_result
    g
  end

  @test typeof(f) == typeof(g) == typeof(h) == PyTest.Fixture

  @pytest function(h)
    @test typeof(f) == PyTest.Fixture
    @test h == :f_result
  end
end

# fixture defined in different module
using different_module
let
  @pytest function(different_module_f, different_module_g)
    @test different_module_f == :dmf_result
    @test different_module_g == [:dmf_result, :dmg_result]
  end
end

# redefinition of fixture
let
  @fixture function f() :one end
  @fixture function f() :two end
  @pytest function(f) @test f == :two end
end

# use of :fixture_name in fixture body?
# not sure if relevant?
let
  @fixture function f() :f_result end
  @fixture function g(f) :f end
  @pytest function(f, g)
    @test f == :f_result
    @test g == :f
  end
end

# use of imported symbols in fixture/test body
using different_module
let
  @fixture function f() @test some_function() == :some_function_result end
  @pytest function(f) @test some_function() == :some_function_result end
end

# teardown gets called in right order
let
  active_objects = []
  @fixture function f()
    @yield push!(active_objects, :f_result)  # next lines are teardown
    pop!(active_objects)
  end
  @fixture function g(f) push!(active_objects, :g_result) end
  @fixture function h(f, g)
    @yield push!(active_objects, :h_result)  # ditto
    push!(active_objects, :h_down)
  end
  @pytest function(h)
    @test active_objects == [:f_result, :g_result, :h_result]
  end
  @test active_objects == [:f_result, :g_result, :h_down]
end

# tears down even if intermediate fixture is noop
let
  deleted = []
  @fixture function f()
    @yield nothing
    push!(deleted, :f_result)
  end
  @fixture function g(f) end
  @pytest function(g) error end
  @test deleted == [:f_result]
end

# tears down also on exception in test
let
  deleted = []
  @fixture function f()
    @yield nothing
    push!(deleted, :f_result)
  end
  @fixture function g(f)
    @yield nothing
    push!(deleted, :g_result)
  end
  @pytest function(g) error end
  @test deleted == [:f_result, :g_result]
end

# tears down also on exception in fixture
let
  deleted = []
  @fixture function f()
    @yield nothing
    push!(deleted, :f_result)
  end
  @fixture function g(f)
    @yield nothing
    push!(deleted, :g_result)
  end
  @fixture function h(g) end
  @pytest function(h) error end
  @test deleted == [:f_result, :g_result]
end


text:
# Compute fast Walsh-Hadamard transforms (in natural, dyadic, or sequency order)
# using FFTW, by interpreting them as 2x2x2...x2x2 DFTs.  We follow Matlab's
# convention in which ifwht is the unnormalized transform and fwht has a 1/N
# normalization (as opposed to using a unitary normalization).

"""
The `Hadamard` module provides functions to compute fast Walsh-Hadamard transforms in Julia,
for arbitrary dimensions and arbitrary power-of-two transform sizes, with the three standard
orderings: natural (Hadamard), dyadic (Paley), and sequency (Walsh) ordering.
See the `fwht`, `fwht_natural`, and `fwht_dyadic` functions, along with their inverses
`ifwht` etc.

There is also a function `hadamard(n)` that returns Hadamard matrices for known sizes
(including non powers of two).
"""
module Hadamard
export fwht, ifwht, fwht_natural, ifwht_natural, fwht_natural!, ifwht_natural!, fwht_dyadic, ifwht_dyadic, hadamard

using FFTW, LinearAlgebra
import FFTW: set_timelimit, dims_howmany, unsafe_execute!, cFFTWPlan, r2rFFTWPlan, PlanPtr, fftwNumber, ESTIMATE, NO_TIMELIMIT, R2HC
import AbstractFFTs: normalization, complexfloat

# A power-of-two dimension to be transformed is interpreted as a
# 2x2x2x....x2x2  multidimensional DFT.  This function transforms
# each individual power-of-two dimension n into the corresponding log2(n)
# DFT dimensions.  If bitreverse is true, then the output is in
# bit-reversed (transposed) order (which may not work in-place).
function hadamardize(dims::Array{Int,2}, bitreverse::Bool)
    ntot = 0
    for i = 1:size(dims,2)
        n = dims[1,i]
        if !ispow2(n)
            throw(ArgumentError("non power-of-two Hadamard-transform length"))
        end
        ntot += trailing_zeros(n)
    end
    hdims = Array{Int}(undef, 3,ntot)
    j = 0
    for i = 1:size(dims,2)
        n = dims[1,i]
        is = dims[2,i]
        os = dims[3,i]
        log2n = trailing_zeros(n)
        krange = j+1:j+log2n
        for k in krange
            hdims[1,k] = 2
            hdims[2,k] = is
            hdims[3,k] = os
            is *= 2
            os *= 2
        end
        if bitreverse
            hdims[3,krange] = reverse(hdims[3,krange])
        end
        j += log2n
    end
    return hdims
end

const libfftw = isdefined(FFTW, :libfftw) ? FFTW.libfftw : FFTW.libfftw3
const libfftwf = isdefined(FFTW, :libfftwf) ? FFTW.libfftwf : FFTW.libfftw3f

for (Tr,Tc,fftw,lib) in ((:Float64,:ComplexF64,"fftw",libfftw),
                         (:Float32,:ComplexF32,"fftwf",libfftwf))
    @eval function Plan_Hadamard(X::StridedArray{$Tc,N}, Y::StridedArray{$Tc,N},
                                 region, flags::Unsigned, timelimit::Real,
                                 bitreverse::Bool) where {N}
        set_timelimit($Tr, timelimit)
        dims, howmany = dims_howmany(X, Y, [size(X)...], region)
        dims = hadamardize(dims, bitreverse)
        plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib[]),
                     PlanPtr,
                     (Int32, Ptr{Int}, Int32, Ptr{Int},
                      Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32),
                     size(dims,2), dims, size(howmany,2), howmany,
                     X, Y, FFTW.FORWARD, flags)
        set_timelimit($Tr, NO_TIMELIMIT)
        if plan == C_NULL
            if $(FFTW.fftw_provider == "mkl")
                error("MKL is not supported — reconfigure FFTW.jl to use FFTW")
            else
                error("FFTW could not create plan") # shouldn't normally happen
            end
        end
        return cFFTWPlan{$Tc,FFTW.FORWARD,X===Y,N}(plan, flags, region, X, Y)
    end

    @eval function Plan_Hadamard(X::StridedArray{$Tr,N}, Y::StridedArray{$Tr,N},
                                 region, flags::Unsigned, timelimit::Real,
                                 bitreverse::Bool) where {N}
        set_timelimit($Tr, timelimit)
        dims, howmany = dims_howmany(X, Y, [size(X)...], region)
        dims = hadamardize(dims, bitreverse)
        kind = Array{Int32}(undef, size(dims,2))
        kind .= R2HC
        plan = ccall(($(string(fftw,"_plan_guru64_r2r")),$lib[]),
                     PlanPtr,
                     (Int32, Ptr{Int}, Int32, Ptr{Int},
                      Ptr{$Tr}, Ptr{$Tr}, Ptr{Int32}, UInt32),
                     size(dims,2), dims, size(howmany,2), howmany,
                     X, Y, kind, flags)
        set_timelimit($Tr, NO_TIMELIMIT)
        if plan == C_NULL
            if $(FFTW.fftw_provider == "mkl")
                error("MKL is not supported — reconfigure FFTW.jl to use FFTW")
            else
                error("FFTW could not create plan") # shouldn't normally happen
            end
        end
        return r2rFFTWPlan{$Tr,(map(Int,kind)...,),X===Y,N}(plan, flags, region, X, Y)
    end
end

############################################################################
# Define ifwht (inverse/unnormalized) transforms for various orderings

# Natural (Hadamard) ordering:

function ifwht_natural(X::StridedArray{<:fftwNumber}, region)
    Y = similar(X)
    p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, false)
    unsafe_execute!(p)
    return Y
end

function ifwht_natural(X::StridedArray{<:Number}, region)
    Y = float(X)
    p = Plan_Hadamard(Y, Y, region, ESTIMATE, NO_TIMELIMIT, false)
    unsafe_execute!(p)
    return Y
end

# Dyadic (Paley, bit-reversed) ordering:

function ifwht_dyadic(X::StridedArray{<:fftwNumber}, region)
    Y = similar(X)
    p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, true)
    unsafe_execute!(p)
    return Y
end

function ifwht_dyadic(X::StridedArray{<:Number}, region)
    return ifwht_dyadic(float(X), region)
end

############################################################################
# Sequency (Walsh) ordering:

# ifwht along a single dimension d of X
function ifwht(X::Array{T}, region) where {T<:fftwNumber}
    Y = ifwht_dyadic(X, region)

    # Perform Gray-code permutation of Y (TODO: in-place?)
    if isempty(region)
        return Y
    elseif ndims(Y) == 1
        return [ Y[1 + ((i >> 1) ⊻ i)] for i = 0:length(Y)-1 ]
    else
        sz = [size(Y)...]
        tmp = Array{T}(undef, maximum(sz[region])) # storage for out-of-place perm.
        for d in region
            # somewhat ugly loops to do 1d permutations along dimension d
            na = prod(sz[d+1:end])
            n = sz[d]
            nb = prod(sz[1:d-1])
            sa = nb * n
            for ia = 0:na-1
                for ib = 1:nb
                    i0 = ib + sa * ia
                    for i = 0:n-1
                        tmp[i+1] = Y[i0 + nb * ((i >> 1) ⊻ i)]
                    end
                    for i = 0:n-1
                        Y[i0 + nb * i] = tmp[i+1]
                    end
                end
            end
        end
        return Y
    end
end

# handle 1d case of strided arrays (loops in multidim case are too annoying)
function ifwht(X::StridedVector{<:fftwNumber}, region)
    Y = ifwht_dyadic(X, region)

    # Perform Gray-code permutation of Y (TODO: in-place?)
    if isempty(region)
        return Y
    else
        return [ Y[1 + ((i >> 1) ⊻ i)] for i = 0:length(Y)-1 ]
    end
end

# fallback for subarrays
function ifwht(X::StridedArray{<:fftwNumber}, region)
    return ifwht(copy(X), region)
end

# fallback for other types
function ifwht(X::StridedArray{<:Number}, region)
    return ifwht(float(X), region)
end

############################################################################
# in-place transforms, currently for natural ordering only.

"""
    ifwht_natural!(X, dims=1:ndims(X))

Similar to `ifwht_natural`, but works in-place on the input array `X`.
"""
function ifwht_natural!(X::StridedArray{<:fftwNumber}, region=1:ndims(X))
    p = Plan_Hadamard(X, X, region, ESTIMATE, NO_TIMELIMIT, false)
    unsafe_execute!(p)
    return X
end

"""
    fwht_natural!(X, dims=1:ndims(X))

Similar to `fwht_natural`, but works in-place on the input array `X`.
"""
fwht_natural!(X::StridedArray{<:fftwNumber}, region=1:ndims(X)) =
    rmul!(ifwht_natural!(X, region), normalization(X,region))

############################################################################
# Forward transforms (normalized by 1/N as in Matlab) and transforms
# without the region argument:

for f in (:ifwht_natural, :ifwht_dyadic, :ifwht)
    g = Symbol(string(f)[2:end])
    @eval begin
        $f(X) = $f(X, 1:ndims(X))
        $g(X) = rmul!($f(X), normalization(X,1:ndims(X)))
        $g(X,r) = rmul!($f(X,r), normalization(X,r))
    end
end

############################################################################
# Docstrings — since they are all very similar, we
# assign them in a loop.

for inverse in (true, false), ord in (:natural, :dyadic, :sequency)
    name = (inverse ? "i" : "") * "fwht"
    if ord != :sequency
        name = string(name, '_', ord)
    end
    ordstr = ord == :natural ? "natural (Hadamard)" :
             ord == :dyadic ? "dyadic (Paley, bit-reversed)" :
             "sequency"
    docstring = """    $name(X, dims=1:ndims(X))

Return the$(inverse ? " inverse" : "") fast Walsh-Hadamard transform (WHT) of
the array `X` along the dimensions `dims` (an integer,
tuple, or array of dimensions, defaulting to all dimensions).
Only power-of-two sizes (along the transformed dimensions)
are supported.  The result is returned in the $ordstr ordering.

Our WHT is normalized so that the forward transform has a `1/n`
coefficient (where `n` is the product of the transformed dimensions)
and the inverse WHT has no scaling factor.
"""
    @eval @doc $docstring $(Symbol(name))
end

############################################################################
# Utilities to work with a precomputed cache of known Hadamard matrices
# of various sizes, produced by util/fetchhadamard.jl from Sloane's web page
# and stored as BitMatrices

function readcache(cachefile::AbstractString)
    B = BitMatrix[]
    open(cachefile, "r") do io
        while !eof(io)
            k = convert(Int, ntoh(read(io, Int64)))
            b = BitArray(undef, k, k)
            b.chunks .= ntoh.(read!(io, b.chunks))
            push!(B, b)
        end
    end
    sort!(B, by = h -> (isodd(size(h,1)>>2), size(h,1)))
    return B
end
readcache() = readcache(joinpath(dirname(@__FILE__), "..", "src", "cache.dat"))

"""
    printsigns([io, ] A)

Print a table of `+` and `-` characters indicating the
signs of the entries of `A` (with `0` for zero entries).
The output stream `io` defaults to `stdout`.
"""
function printsigns(io::IO, A::AbstractMatrix{<:Real})
    m, n = size(A)
    println(io, m, "×", n, " sign matrix from ", typeof(A))
    for i = 1:m
        for j = 1:n
            print(io, A[i,j] > 0 ? '+' : A[i,j] < 0 ? '-' : '0')
        end
        println(io)
    end
end
printsigns(A) = printsigns(stdout, A)

function frombits(B::BitMatrix)
    A = Matrix{Int8}(B)
    for i = 1:length(A)
        A[i] = A[i] == 0 ? -1 : 1
    end
    return A
end

const hadamards = BitMatrix[] # read lazily below

############################################################################
# Create Int8 order-n Hadamard matrix (in natural order), by factorizing
# n into a product of known Hadamard sizes (if possible)

const H2 = Int8[1 1; 1 -1]

"""
    hadamard(n)

Return a Hadamard matrix of order `n`.  The known Hadamard
matrices up to size 256 are currently supported (via a lookup table),
along with any size that factorizes into products of these known
sizes and/or powers of two.

The return value is a matrix of `Int8` values. If you want to do
further matrix computations with this matrix, you may want to
convert to `Float64` first via `float(hadamard(n))`.

You can pretty-print a Hadamard matrix as a table of `+` and `-`
(indicating the signs of the entries) via `Hadamard.printsigns`,
e.g. `Hadamard.printsigns(hadamard(28))` for the 28×28 Hadamard matrix.
"""
function hadamard(n::Integer)
    n < 1 && throw(ArgumentError("size n=$n should be positive"))
    n0 = n
    H = reshape(Int8[1], 1,1)
    if !ispow2(n)
        isempty(hadamards) && append!(hadamards, readcache())
        for i = length(hadamards):-1:1
            k = size(hadamards[i], 1)
            if rem(n, k) == 0
                Hk = frombits(hadamards[i])
                while true
                    H = kron(H, Hk)
                    n = div(n, k)
                    rem(n, k) != 0 && break
                end
            end
        end
    end
    if !ispow2(n)
        n >>= trailing_zeros(n)
        throw(ArgumentError("unknown Hadamard factor $n of $n0"))
    end
    # Note: it would be faster to do a "power-by-squaring" like algorithm
    # here, where we repeatedly double the order via kron(H, H), but
    # it's not clear to me that we care much about performance here.
    while n > 1
        H = kron(H, H2)
        n >>= 1
    end
    return H
end

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

end # module


text:
# ------------ Choice sites ------------ #

@inline function (ctx::ScoreContext)(call::typeof(rand),
                                     addr::T,
                                     d::Distribution{K}) where {T <: Address, K}
    visit!(ctx, addr)
    has_top(ctx.select, addr) || error("ScoreError: constrained selection must provide constraints for all possible addresses in trace. Missing at address $addr.")
    val = get_top(ctx.select, addr)
    increment!(ctx, logpdf(d, val))
    return val
end

# ------------ Learnable ------------ #

@inline function (ctx::ScoreContext)(fn::typeof(learnable), addr::Address)
    visit!(ctx, addr)
    has_top(ctx.params, addr) && return get_top(ctx.params, addr)
    error("Parameter not provided at address $addr.")
end

# ------------ Call sites ------------ #

@inline function (ctx::ScoreContext)(c::typeof(rand),
                                     addr::T,
                                     call::Function,
                                     args...) where T <: Address
    visit!(ctx, addr)
    ps = get_subparameters(ctx, addr)
    ss = get_subselection(ctx, addr)
    ret, w = score(ss, ps, call, args...)
    increment!(ctx, w)
    return ret
end


text:
using NLsolve: OnceDifferentiable, NewtonCache, newton_
using LineSearches

struct NLsolveNewton{T, FT, DT, CT, ST, LT} <: AbstractNewtonSolver{T}
    x::Vector{T}
    f::Vector{T}
    J::Matrix{T}

    cache::CT

    F!::FT
    DF::DT

    line_search::ST
    linear_solver::LT

    params::NonlinearSolverParameters{T}
    status::NonlinearSolverStatus{T}

    function NLsolveNewton(x::AbstractVector{T}, f::AbstractVector{T}, J::AbstractMatrix{T},
                    F!::FT, DF::DT, cache::CT, line_search::ST, linear_solver::LT) where {T,FT,DT,CT,ST,LT}

        nls_params = NonlinearSolverParameters(T)
        nls_status = NonlinearSolverStatus{T}(length(x))

        new{T,FT,DT,CT,ST,LT}(x, f, J, cache, F!, DF, line_search, linear_solver, nls_params, nls_status)
    end
end

function NLsolveNewton(x::AbstractVector{T}, F!::Function; J!::Union{Function,Nothing}=nothing) where {T}
    linear_solver = getLinearSolver(x)

    f  = linear_solver.b
    df = linear_solver.A
    f! = (f, x) -> F!(x, f)

    if J! == nothing
        if get_config(:jacobian_autodiff)
            DF = OnceDifferentiable(f!, x, f, df; autodiff=:forward, inplace=true)
        else
            DF = OnceDifferentiable(f!, x, f, df; autodiff=:finite, inplace=true)
        end
    else
        j! = (df, x) -> J!(x, df)
        DF = OnceDifferentiable(f!, j!, x, f, df; inplace=true)
    end

    NLsolveNewton(x, f, df, F!, DF, NewtonCache(DF), LineSearches.Static(), linear_solver)
end

function linsolve!(s::NLsolveNewton, x, A, b)
    copyto!(s.linear_solver.A, A)
    copyto!(s.linear_solver.b, b)
    factorize!(s.linear_solver)
    solve!(s.linear_solver)
    copyto!(x, s.linear_solver.b)
end

function solve!(s::NLsolveNewton; n::Int=0)
    local nmax::Int = n > 0 ? nmax = n : s.params.nmax

    res=newton_(s.DF, s.x, s.params.stol, s.params.atol, nmax, false, false, false,
                s.line_search, (x, A, b) -> linsolve!(s, x, A, b), s.cache)

    s.x .= res.zero
    s.status.i  = res.iterations
    s.status.rₐ = res.residual_norm

    nothing
end


text:
module Trees
export PatternTree, PatternStep, PatternCheck,
       PatternLeaf, PatternNode, PatternGate,
       PatternRoot, PatternHead, Binding,

       ExprHead, SlurpHead,

       nodehead, bindings, insert!, newnode!,
       newleaf!, depth

#-----------------------------------------------------------------------------------
# Type definitions
#-----------------------------------------------------------------------------------

abstract type PatternTree  end
abstract type PatternStep  end
abstract type PatternCheck end
abstract type PatternHead  end

mutable struct PatternRoot <: PatternTree
  child::PatternTree
  PatternRoot() = new()
end

struct PatternNode{PH <: PatternHead} <: PatternTree
  head     :: PH
  step     :: PatternStep
  children :: Vector{PatternTree}
  bindings :: Set{Symbol}
  depth    :: Int
end

struct PatternLeaf <: PatternTree
end

struct Binding <: PatternCheck
  name::Symbol
end

mutable struct PatternGate <: PatternTree
  check    :: PatternCheck
  bindings :: Set{Symbol}
  depth    :: Int
  child    :: PatternTree

  PatternGate(check::Any,     depth) = new(check, Set{Symbol}(), depth)
  PatternGate(check::Binding, depth) = new(check, Set{Symbol}([check.name]), depth)
end

abstract type SlurpHead <: PatternHead end

struct ExprHead <: PatternHead
  sym::Symbol
end

const SingleChildNode = Union{PatternGate, PatternRoot}

#-----------------------------------------------------------------------------------
# Functions
#-----------------------------------------------------------------------------------

function nodehead(node::PatternNode)
  isa(node.head, ExprHead) ? node.head.sym : (:slurp)
end

bindings(leaf::PatternLeaf) = Set{Symbol}()
bindings(gate::PatternGate) = gate.bindings
bindings(node::PatternNode) = node.bindings

depth(node::PatternRoot) = 0
depth(gate::PatternGate) = gate.depth
depth(node::PatternNode) = node.depth

function makenode(head, step, depth)
  children   = PatternTree[]
  bindings   = Set{Symbol}()
  slurpdepth = isa(head, SlurpHead) ? depth+1 : depth

  PatternNode(head, step, children, bindings, slurpdepth)
end

import Base: insert!

function insert!(parent::PatternNode, child)
  push!(parent.children, child)
  union!(parent.bindings, bindings(child))
end

function insert!(parent::PatternGate, child)
  parent.child = child
  union!(parent.bindings, bindings(child))
end

function insert!(parent::PatternRoot, child)
  parent.child = child
end

function newnode!(head, step, parent::PatternTree)
  node = makenode(head, step, depth(parent))
  insert!(parent, node)
  return node
end

function newleaf!(check, parent::PatternTree)
  leaf = PatternLeaf()
  gate = PatternGate(check, depth(parent))
  insert!(gate, leaf)
  insert!(parent, gate)
  return leaf
end

#-----------------------------------------------------------------------------------
end


text:
function Integrate(fun::Function, xl::Point{Float64}, xu::Point{Float64}, xl_p::Point{Float64}, xu_p::Point{Float64}, cal::Int64, iterations::Int64, dim::Float64, time_stamp::Float64)

	# First Step: cal evaluations of the integrand using uniform distributed random number
	rng = MersenneTwister(0)
	#la dimensione non va bene
	pt = [Point{Float64}(rand(rng), rand(rng), rand(rng), undef, undef, undef)() for n = 1:Int64(floor(3. * cal))]
	pt_p = [Point{Float64}(rand(rng), rand(rng), rand(rng), undef, undef, undef)() for n = 1:Int64(floor(3. * cal))]

	#Deve farmi la differenza degli r, theta e phi e calcolarmi le x,y,z di conseguenza
	# Coordinate transformation. ∫ₐᵇ f(x) dx  = ∫₀¹ f(y *(b-a) +a )*(b-a) dy
	p   = ( pt .* Ref(xu-xl) ) .+ Ref(xl)
	p_p = ( pt_p .* Ref(xu_p -xl_p) ) .+ Ref(xl)
	# Jacobian of the translation
	J_t = jac_trasl(xl, xu, xl_p, xu_p)

	density_unif = 1. / (Int64(floor(3. * calls)))

	# Integrand evaluations
	f = [fun(p[i], p_p[i], time_stamp) * J_t * density_unif  for i = 1:Int64(floor(3. * cal)) ]
	f_abs = abs.(f)

	# Perform the summation of intergrand and its absolute value
	Integral      = sum(f)
	Integral_abs  = sum(f_abs)

	# Check value - \Int_{-1}^{1} gauss_3D(x) dx
	#Int_True = 0.00378269
	#println("0. Integral = ", Integral)

	# Estimation of the new weight function
	g = (f_abs ./ Integral_abs)

	#println("Check PDF Integral = ",sum(g))
	# I have to study hard to find a more suitable method to estimate the binning
	NBIN = 3

	# PDF approximation with stepper function
	w_m	 = zeros(Int64(NBIN^dim))
	Int_box  = zeros(Int64(NBIN^dim))
	Int_abs_box = zeros(Int64(NBIN^dim))
	dix = diy = diz = div = diw = dik = 1.0/NBIN
	Integral_avg = 0

	# Box_R = Array{Float64}(undef, 0)
	# Box_T = Array{Float64}(undef, 0)
	# Box_P = Array{Float64}(undef, 0)
	#
	# Box_R_P = Array{Float64}(undef, 0)
	# Box_T_P = Array{Float64}(undef, 0)
	# Box_P_P = Array{Float64}(undef, 0)

	# HIT = Array{Float64}(undef, 0)

	chi_square = 0

	for iter = 1:iterations
		# Evaluate the weigth on each BOX
		for idx_g = 1:length(g)
			box_id_x = Int64(ceil(pt[idx_g].r / dix)) -1
			box_id_y = Int64(ceil(pt[idx_g].θ / diy)) -1
			box_id_z = Int64(ceil(pt[idx_g].ϕ / diz)) -1

			box_id_v = Int64(ceil(pt_p[idx_g].r / dix)) -1
			box_id_w = Int64(ceil(pt_p[idx_g].θ / diy)) -1
			box_id_k = Int64(ceil(pt_p[idx_g].ϕ / diz)) -1
			# Qui si puo` forse scrivere meglio
			ID = (box_id_x) + (box_id_y)*NBIN + (box_id_z)*NBIN*NBIN + (box_id_v)*NBIN*NBIN*NBIN + (box_id_w)*NBIN*NBIN*NBIN*NBIN + (box_id_k)*NBIN*NBIN*NBIN*NBIN*NBIN

			w_m[ID+1] += g[idx_g]
		end

		# Fino a qui mi sembra tutto in ordine. Gli indici mi tornano tutti.

		# Check the weights distribution on each box and check if the
		# PDF requirements is satisfied \Int_{\Omega} pdf(x) dx = 1

		#println("Densita` in BOX")
		#for i in 1:length(w_m)
			#println(w_m[i])
		#end
		#println("Check PDF Integral = ", sum(w_m))

		# Evaluate the Integral using a non uniform random number distribution
		total_hits = 0
		pt_2_all   = Array{Point{Float64}}(undef, 0)
		pt_2_p_all = Array{Point{Float64}}(undef, 0)
		f_2_all    = Array{Float64,1}(undef, 0)

		for ix = (1:NBIN)
			for iy = (1:NBIN)
				for iz = (1:NBIN)
					for iv = (1:NBIN)
						for iw = (1:NBIN)
							for ik = (1:NBIN)
								p_start = Point{Float64}((ix-1) * dix , (iy-1) * diy, (iz-1) * diz, undef, undef, undef)
								p_final = Point{Float64}(ix * dix, iy * diy, iz * diz, undef, undef, undef)
								p_p_start = Point{Float64}((iv-1) * div , (iw-1) * diw, (ik-1) * dik, undef, undef, undef)
								p_p_final = Point{Float64}(iv * div, iw * diw, ik * dik, undef, undef, undef)

								# Si puo` semplificare la scrittura?
								ID_BOX = ( (ix-1) + ( (iy-1) * NBIN ) + ( (iz-1) * NBIN * NBIN ) +  ( (iv-1) * NBIN * NBIN * NBIN ) + ( (iw-1) * NBIN * NBIN * NBIN * NBIN ) + ( (ik-1) * NBIN * NBIN * NBIN * NBIN * NBIN)    ) + 1

								# Uniform density
								#density    = 1/(NBIN*NBIN*NBIN)

								# Non uniform weight PDF
								density    = w_m[ID_BOX]

								hit_in_box = Int64(floor(density * (3. * cal)))

								total_hits += hit_in_box
								#println("Hit = ", hit_in_box)
								if hit_in_box != 0

									x_ran = (p_final.r-p_start.r)
									y_ran = (p_final.θ-p_start.θ)
									z_ran = (p_final.ϕ-p_start.ϕ)

									v_ran = (p_p_final.r-p_p_start.r)
									w_ran = (p_p_final.θ-p_p_start.θ)
									k_ran = (p_p_final.ϕ-p_p_start.ϕ)

									pt_2 = [Point{Float64}(x_ran*rand(rng)+p_start.r, y_ran*rand(rng)+p_start.θ, z_ran*rand(rng)+p_start.ϕ, undef,undef,undef)() for n = 1:hit_in_box]
									pt_p_2 = [Point{Float64}(v_ran*rand(rng)+p_p_start.r, w_ran*rand(rng)+p_p_start.θ, k_ran*rand(rng)+p_p_start.ϕ,undef,undef,undef)() for n = 1:hit_in_box]

									p_2 = ( pt_2 .* Ref(xu-xl) ) .+ Ref(xl)
									p_p_2 = ( pt_p_2 .* Ref(xu_p-xl_p) ) .+ Ref(xl_p)
									# Aggiustare quel /density!!!!!! Non e` conforme con quello sopra!!!!
									f_2 = [fun(p_2[i], p_p_2[i], time_stamp) * J_t / density  for i = 1:hit_in_box ]

									append!(pt_2_all, pt_2)
									append!(pt_2_p_all, pt_p_2)

									append!(f_2_all, f_2)

									# append!(Box_R, ix)
									# append!(Box_T, iy)
									# append!(Box_P, iz)
									#
									# append!(Box_R_P, iv)
									# append!(Box_T_P, iw)
									# append!(Box_P_P, ik)
									# append!(HIT,hit_in_box)

								else
									continue
								end
							end
						end
					end
				end
			end
		end
		Integral += sum(f_2_all) / (total_hits * NBIN^6 )
		Integral_abs = sum(abs.(f_2_all))
		g = (abs.(f_2_all) ./ Integral_abs)
		w_m = zeros(Int64(NBIN^6))

		pt = pt_2_all
		pt_p = pt_2_p_all

		#println(iter,". Integral = ", Integral / (iter + 1)  )
		Integral_avg = Integral / (iter+1)
		chi_square += (sum(f_2_all) / (total_hits * NBIN^6 ) - Integral_avg)^2

		# println("#######################################################################")
		# println( "r_min = ",     ((findmin(Box_R)[1]-1)/NBIN * (xu.r-xl.r)) + xl.r   )
		# println( "r_max = ",    ((findmax(Box_R)[1])/NBIN * (xu.r-xl.r)) + xl.r  )
		#
		# println( "θ_min = ",   ((findmin(Box_T)[1]-1)/NBIN * (xu.θ-xl.θ)) + xl.θ   )
		# println( "θ_max = ",    ((findmax(Box_T)[1])/NBIN * (xu.θ-xl.θ)) + xl.θ  )
		#
		# println( "ϕ_min = ",    ((findmin(Box_P)[1]-1)/NBIN * (xu.ϕ-xl.ϕ)) + xl.ϕ   )
		# println( "ϕ_min = ",   ((findmax(Box_P)[1])/NBIN * (xu.ϕ-xl.ϕ)) + xl.ϕ  )
		#
		# println("################################## PRIMED #############################")
		#
		# println( "r_min = ",   ((findmin(Box_R_P)[1]-1)/NBIN * (xu_p.r-xl_p.r)) + xl_p.r   )
		# println( "r_max = ",    ((findmax(Box_R_P)[1])/NBIN * (xu_p.r-xl_p.r)) + xl_p.r  )
		#
		# println( "θ_min = ",   ((findmin(Box_T_P)[1]-1)/NBIN * (xu_p.θ-xl_p.θ)) + xl_p.θ   )
		# println( "θ_max = ",   ((findmax(Box_T_P)[1])/NBIN * (xu_p.θ-xl_p.θ)) + xl_p.θ  )
		#
		# println( "ϕ_min = ",   ((findmin(Box_P_P)[1]-1)/NBIN * (xu_p.ϕ-xl_p.ϕ)) + xl_p.ϕ   )
		# println( "ϕ_max = ",   ((findmax(Box_P_P)[1])/NBIN * (xu_p.ϕ-xl_p.ϕ)) + xl_p.ϕ  )
		# println("#######################################################################")
		#
		#
		# Box_R = Array{Float64}(undef, 0)
		# Box_T = Array{Float64}(undef, 0)
		# Box_P = Array{Float64}(undef, 0)
		#
		# Box_R_P = Array{Float64}(undef, 0)
		# Box_T_P = Array{Float64}(undef, 0)
		# Box_P_P = Array{Float64}(undef, 0)

	end

	return Integral_avg, chi_square #, HIT
end

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
module CrystalStructures

using   Printf
using   StaticArrays
using   LinearAlgebra
import  Crystalline

# Number of space groups
# OEIS: A006227
const NO_SGS = (2, 17, 230, 4985, 222097)
# Number of space groups, with chiral copies not distinct
# OEIS: A004029
const NO_SGS_PAIRED = (2, 17, 219, 4783, 222018, 28927915)

# List of elements
# You can convert a vector of element numbers to a vector of element names with ELEMENTS[v]
const ELEMENTS =
[
    "H",  "He", "Li", "Be", "B",  "C",  "N",  "O",  "F",  "Ne",
    "Na", "Mg", "Al", "Si", "P",  "S",  "Cl", "Ar", "K",  "Ca",
    "Sc", "Ti", "V",  "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
    "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y",  "Zr",
    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
    "Sb", "Te", "I",  "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
    "Lu", "Hf", "Ta", "W",  "Re", "Os", "Ir", "Pt", "Au", "Hg",
    "Tl", "Pb", "Bi", "Po", "As", "Rn", "Fr", "Ra", "Ac", "Th",
    "Pa", "U",  "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
    "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
    "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"
]

# TODO: can we modify this to make it unnecessary?
const ELEMENT_LOOKUP = Dict{String, Int}([(ELEMENTS[n] => n) for n in 1:length(ELEMENTS)])

# Metallic radii of the elements in angstroms
# Probably will be useful for sphere packing related things
const RADII_METALLIC_ANGSTROMS = Dict{String, Float64}(
    "Li" => 1.52,
    "Be" => 1.12,
    "Na" => 1.86,
    "Mg" => 1.60,
    "Al" => 1.43,
    "K"  => 2.27,
    "Ca" => 1.97,
    "Sc" => 1.62,
    "Ti" => 1.47,
    "V"  => 1.34,
    "Cr" => 1.28,
    "Mn" => 1.27,
    "Fe" => 1.26,
    "Co" => 1.25,
    "Ni" => 1.24,
    "Cu" => 1.28,
    "Zn" => 1.34,
    "Ga" => 1.35,
    "Rb" => 2.48,
    "Sr" => 2.15,
    "Y"  => 1.80,
    "Zr" => 1.60,
    "Nb" => 1.46,
    "Mo" => 1.39,
    "Tc" => 1.36,
)

include("vectors.jl")
include("datagrids.jl")
include("crystals.jl")
include("reciprocal.jl")
include("filetypes.jl")

# Types
export  AbstractCrystal, CrystalStructure, AbstractCrystalData, CrystalStructureWithData,
        AbstractKPoints, KPointGrid, KPointList
# Functions
export  formula, readXYZ, writeXYZ, readXSF3D, writeXSF3D

end # module


text:
@testset "Randgraphs" begin
    r1 = Graph(10,20)
    r2 = DiGraph(5,10)

    @test nv(r1) == 10
    @test ne(r1) == 20
    @test nv(r2) == 5
    @test ne(r2) == 10
    @test eltype(r1) == Int
    @test eltype(r2) == Int

    @test eltype(Graph(0x5, 0x2)) == eltype(Graph(0x5, 2)) == UInt8
    for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt, Int]
        @test eltype(Graph{T}(5,2)) == T
        @test eltype(DiGraph{T}(5,2)) == T
    end

    @test Graph(10,20,seed=3) == Graph(10,20,seed=3)
    @test DiGraph(10,20,seed=3) == DiGraph(10,20,seed=3)
    @test Graph(10,20,seed=3) == erdos_renyi(10,20,seed=3)
    @test ne(Graph(10,40,seed=3)) == 40
    @test ne(DiGraph(10,80,seed=3)) == 80

    er = erdos_renyi(10, 0.5)
    @test nv(er) == 10
    @test is_directed(er) == false
    er = erdos_renyi(10, 0.5, is_directed=true)
    @test nv(er) == 10
    @test is_directed(er) == true

    er = erdos_renyi(10, 0.5, seed=17)
    @test nv(er) == 10
    @test is_directed(er) == false

    ws = watts_strogatz(10,4,0.2)
    @test nv(ws) == 10
    @test ne(ws) == 20
    @test is_directed(ws) == false

    ws = watts_strogatz(10, 4, 0.2, is_directed=true)
    @test nv(ws) == 10
    @test ne(ws) == 20
    @test is_directed(ws) == true

    ba = barabasi_albert(10, 2)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, 2)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 4, 2)
    @test nv(ba) == 10
    @test ne(ba) == 12
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 17
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 17
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 4, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 4, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 12
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 4, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 24
    @test is_directed(ba) == true

    fm = static_fitness_model(20, rand(10))
    @test nv(fm) == 10
    @test ne(fm) == 20
    @test is_directed(fm) == false

    fm = static_fitness_model(20, rand(10), rand(10))
    @test nv(fm) == 10
    @test ne(fm) == 20
    @test is_directed(fm) == true

    sf = static_scale_free(10, 20, 2.0)
    @test nv(sf) == 10
    @test ne(sf) == 20
    @test is_directed(sf) == false

    sf = static_scale_free(10, 20, 2.0, 2.0)
    @test nv(sf) == 10
    @test ne(sf) == 20
    @test is_directed(sf) == true

    rr = random_regular_graph(5, 0)
    @test nv(rr) == 5
    @test ne(rr) == 0
    @test is_directed(rr) == false

    rd = random_regular_digraph(10,0)
    @test nv(rd) == 10
    @test ne(rd) == 0
    @test is_directed(rd)

    rr = random_regular_graph(6, 3, seed=1)
    @test nv(rr) == 6
    @test ne(rr) == 9
    @test is_directed(rr) == false

    rr = random_regular_graph(1000, 50)
    @test nv(rr) == 1000
    @test ne(rr) == 25000
    @test is_directed(rr) == false
    for v in vertices(rr)
        @test degree(rr, v) == 50
    end

    rr = random_configuration_model(10, repmat([2,4] ,5), seed=3)
    @test nv(rr) == 10
    @test ne(rr) == 15
    @test is_directed(rr) == false
    num2 = 0; num4 = 0
    for v in vertices(rr)
        d = degree(rr, v)
        @test  d == 2 || d == 4
        d == 2 ? num2 += 1 : num4 += 1
    end
    @test num4 == 5
    @test num2 == 5

    rr = random_configuration_model(1000, zeros(Int,1000))
    @test nv(rr) == 1000
    @test ne(rr) == 0
    @test is_directed(rr) == false

    rr = random_configuration_model(3, [2,2,2], check_graphical=true)
    @test nv(rr) == 3
    @test ne(rr) == 3
    @test is_directed(rr) == false

    rd = random_regular_digraph(1000, 4)
    @test nv(rd) == 1000
    @test ne(rd) == 4000
    @test is_directed(rd)
    @test std(outdegree(rd)) == 0

    rd = random_regular_digraph(1000, 4, dir=:in)
    @test nv(rd) == 1000
    @test ne(rd) == 4000
    @test is_directed(rd)
    @test std(indegree(rd)) == 0

    rr = random_regular_graph(10, 8, seed=4)
    @test nv(rr) == 10
    @test ne(rr) == 40
    @test is_directed(rr) == false
    for v in vertices(rr)
        @test degree(rr, v) == 8
    end

    rd = random_regular_digraph(10, 8, dir=:out, seed=4)
    @test nv(rd) == 10
    @test ne(rd) == 80
    @test is_directed(rd)

    g = stochastic_block_model(2., 3., [100,100])
    @test  4.3 < mean(degree(g)) < 5.8
    g = stochastic_block_model(3., 4., [100,100,100])
    @test  10.3 < mean(degree(g)) < 11.8

    function generate_nbp_sbm(numedges, sizes)
        density = 1
        # print(STDERR, "Generating communites with sizes: $sizes\n")
        between = density * 0.90
        intra = density * -0.005
        noise = density * 0.00501
        sbm = nearbipartiteSBM(sizes, between, intra, noise)
        edgestream = make_edgestream(sbm)
        g = Graph(sum(sizes), numedges, edgestream)
        return sbm, g
    end

    function test_sbm(sbm, bp)
        @test sum(sbm.affinities) != NaN
        @test all(sbm.affinities .> 0)
        @test sum(sbm.affinities) != 0
        @test all(bp .>= 0)
        @test all(bp .!= NaN)
    end

    numedges = 100
    sizes = [10, 10, 10, 10]

    n = sum(sizes)
    sbm, g = generate_nbp_sbm(numedges, sizes)
    @test ne(g) >= 0.9numedges
    bc = blockcounts(sbm, g)
    bp = blockfractions(sbm, g) ./ (sizes * sizes')
    ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
    test_sbm(sbm, bp)
    @test norm(collect(ratios)) < 0.25

    sizes = [200, 200, 100]
    internaldeg = 15
    externaldeg = 6
    internalp = Float64[internaldeg/i for i in sizes]
    externalp = externaldeg/sum(sizes)
    numedges = internaldeg + externaldeg #+ sum(externaldeg.*sizes[2:end])
    numedges *= div(sum(sizes), 2)
    sbm = StochasticBlockModel(internalp, externalp, sizes)
    g = Graph(sum(sizes), numedges, sbm)
    @test ne(g) >= 0.9numedges
    @test ne(g) <= numedges
    @test nv(g) == sum(sizes)
    bc = blockcounts(sbm, g)
    bp = blockfractions(sbm, g) ./ (sizes * sizes')
    test_sbm(sbm, bp)
    ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
    @test norm(collect(ratios)) < 0.25

    # check that average degree is not too high
    # factor of two is cushion for random process
    @test mean(degree(g)) <= 4//2*numedges/sum(sizes)
    # check that the internal degrees are higher than the external degrees
    # 5//4 is cushion for random process.
    @test all(sum(bc-diagm(diag(bc)), 1) .<= 5//4 .* diag(bc))

    sbm2 = StochasticBlockModel(0.5*ones(4), 0.3, 10*ones(Int,4))
    sbm  = StochasticBlockModel(0.5, 0.3, 10, 4)
    @test sbm == sbm2
    sbm.affinities[1,1] = 0
    @test sbm != sbm2

    kg = @inferred kronecker(5,5)
    @test nv(kg) == 32
    @test is_directed(kg)
end


text:
# Ported from Radford Neal's R code, with a few thinggs missing

# Arguments:
#
#   x0    Initial point
#   g     Function returning the log of the probability density (plus constant)
#   w     Size of the steps for creating interval (default 1)
#   m     Limit on steps
#   lower Lower bound on support of the distribution (default -Inf)
#   upper Upper bound on support of the distribution (default +Inf)
#   gx0   g(x0)
#
function slice_sampler(x0::Float64, g::Function, w::Float64, m::Int64, lower::Float64, upper::Float64, gx0::Float64)

  if w <= 0
    error("Negative w not allowed")
  end
  if m <= 0
    error("Limit on steps must be positive")
  end
  if upper < lower
    error("Upper limit must be above lower limit")
  end

  # Determine the slice level, in log terms.

  logy::Float64 = gx0 - rand(Exponential(1.0))

  # Find the initial interval to sample from.

  u::Float64 = rand() * w
  L::Float64 = x0 - u
  R::Float64 = x0 + (w-u)

  # Expand the interval until its ends are outside the slice, or until
  # the limit on steps is reached.

  J::Int64 = floor(rand() * m)
  K::Int64 = (m-1) - J

  while J > 0
    if L <= lower || g(L)::Float64 <= logy
      break
    end
    L -= w
    J -= 1
  end

  while K > 0
    if R >= upper || g(R)::Float64 <= logy
      break
    end
    R += w
    K -= 1
  end

  # Shrink interval to lower and upper bounds.

  L = L < lower ? lower : L
  R = R > upper ? upper : R

  # Sample from the interval, shrinking it on each rejection.

  x1::Float64 = 0.0  # need to initialize it in this scope first
  gx1::Float64 = 0.0
  while true
    x1 = rand() * (R-L) + L
    gx1 = g(x1)::Float64
    if gx1 >= logy
      break
    end
    if x1 > x0
      R = x1
    else
      L = x1
    end
  end

  return x1,gx1
end

function slice_sampler(x0::Float64, g::Function, w::Float64, m::Int64, lower::Float64, upper::Float64)
  gx0 = g(x0)
  slice_sampler(x0,g,w,m,lower,upper,gx0)
end

function slice_sampler(x0::Float64, g::Function, gx0::Float64)
  slice_sampler(x0,g,.5,10000,-Inf,Inf,gx0)
end

function slice_sampler(x0::Float64, g::Function)
  slice_sampler(x0,g,.5,10000,-Inf,Inf)
end

# Interface using Density type

slice_sampler(x0::Float64, d::Density, w::Float64, m::Int64, lower::Float64, upper::Float64, gx0::Float64) = slice_sampler(x0::Float64, d.f, w::Float64, m::Int64, lower::Float64, upper::Float64, gx0::Float64)

slice_sampler(x0::Float64, d::Density, gx0::Float64) = slice_sampler(x0::Float64, d.f, gx0::Float64)

function slice_sampler(x0::Float64, d::Density, opts::Options)
  gx0 = g(x0)
  slice_sampler(x0::Float64, d.f, opts[:w],opts[:m],opts[:lower],opts[:upper],gx0::Float64)
end


text:
abstract type AbstractAltPolicy end

Base.string(a::AbstractAltPolicy) = string(a.key)

"""
```
mutable struct AltPolicy
```

Types defining an alternative policy rule.

### Fields

- `key::Symbol`: alternative policy identifier

- `eqcond::Function`: a version of `DSGE.eqcond` which computes the equilibrium
  condition matrices under the alternative policy. Like `DSGE.eqcond`, it should
  take in one argument of mutable struct `AbstractDSGEModel` and return the `Γ0`,
  `Γ1`, `C`, `Ψ`, and `Π` matrices.

- `solve::Function`: a version of `DSGE.solve` which solves the model under the
  alternative policy. Like `DSGE.solve`, it should take in one argument of mutable
  struct `AbstractDSGEModel` and return the `TTT`, `RRR`, and `CCC` matrices.

- `forecast_init::Function`: a function that initializes forecasts under the
  alternative policy rule. Specifically, it accepts a model, an `nshocks` x
  `n_forecast_periods` matrix of shocks to be applied in the forecast, and a
  vector of initial states for the forecast. It must return a new matrix of
  shocks and a new initial state vector. If no adjustments to shocks or initial
  state vectors are necessary under the policy rule, this field may be omitted.

- `color::Colorant`: color to plot this alternative policy in. Defaults to blue.

- `linestyle::Symbol`: line style for forecast plots under this alternative
  policy. See options from `Plots.jl`. Defaults to `:solid`.
"""
mutable struct AltPolicy <: AbstractAltPolicy
    key::Symbol
    eqcond::Function
    solve::Function
    forecast_init::Function
    color::Colorant
    linestyle::Symbol
end

function AltPolicy(key::Symbol, eqcond_fcn::Function, solve_fcn::Function;
                   forecast_init::Function = identity,
                   color::Colorant = RGB(0., 0., 1.),
                   linestyle::Symbol = :solid)

    AltPolicy(key, eqcond_fcn, solve_fcn, forecast_init, color, linestyle)
end

"""
```
get_altpolicy(key::Symbol)
```

returns the `AltPolicy` for a alternative policy identified by `key`.
"""
function get_altpolicy(key::Symbol)
    altpol = if key == :flexible_ait
        flexible_ait()
    elseif key == :zero_rate
        zero_rate()
    elseif key == :zlb_rule
        zlb_rule()
    elseif key == :default_policy
        default_policy()
    elseif key == :taylor_rule
        taylor_rule()
    elseif key == :smooth_ait_gdp_alt
        smooth_ait_gdp_alt()
    elseif key == :smooth_ait_gdp
        smooth_ait_gdp()
    elseif key == :ait
        ait()
    elseif key == :ngdp
        ngdp()
    elseif key == :rw_zero_rate
        rw_zero_rate()
    elseif key == :rw
        rw()
    else
        error("AltPolicy $(key) has not been added to `get_altpolicy` as having a method for replacing equilibrium conditions")
    end
end

"""
```
mutable struct EqcondEntry
```

Type to hold the entries in the regime_eqcond_info dictionary for
alternative policies, regime switching, and imperfect awareness.
"""
mutable struct EqcondEntry
    alternative_policy::Union{AltPolicy, Missing}
    weights::Union{Array{Float64, 1}, Missing}
end

function EqcondEntry()
    return EqcondEntry(missing, missing)
end

function EqcondEntry(altpol::AltPolicy)
    return EqcondEntry(altpol, missing)
end

"""
```
mutable struct MultiPeriodAltPolicy
```

Types defining an alternative policy rule.

### Fields

- `key::Symbol`: alternative policy identifier

- `regime_eqcond_info::AbstractDict{Int64, EqcondEntry}`: a dictionary mapping
    regimes to equilibrium conditions which replace the default ones in a
    given regime.

- `gensys2::Bool`: if true, the multi-period alternative policy needs
    to call `gensys2` instead of `gensys` to work.

- `temporary_altpolicy_names::Union{Vector{Symbol}, Nothing}`: specifies
    the names of temporary policies which may occur, e.g. `[:zero_rate]`
    if a temporary ZLB is implemented using the `zero_rate` `AltPolicy`.

- `temporary_altpolicy_length::Int64`: the temporary alternative policy's length
    which is used to determine which regimes need to be solved using `gensys2`
    instead of `gensys`.

- `infoset::Union{Vector{UnitRange{Int64}}, Nothing}`: either a vector specifying
    the information set used for expectations in the measurement equation or
    `nothing` to indicate myopia in expectations across regimes.

- `perfect_credibility_identical_transitions::Union{Dict{Int64, Int64}, Nothing}`: if not a `Nothing`,
    then this field specifies regimes which have identical transition equations
    in the case of perfect credibility.

- `identical_eqcond_regimes::Union{Dict{Int64, Int64}, Nothing}`: if not a `Nothing`,
    then this field specifies regimes which have identical equilibrium conditions.

- `forecast_init::Function`: a function that initializes forecasts under the
  alternative policy rule. Specifically, it accepts a model, an `nshocks` x
  `n_forecast_periods` matrix of shocks to be applied in the forecast, and a
  vector of initial states for the forecast. It must return a new matrix of
  shocks and a new initial state vector. If no adjustments to shocks or initial
  state vectors are necessary under the policy rule, this field may be omitted.

- `color::Colorant`: color to plot this alternative policy in. Defaults to blue.

- `linestyle::Symbol`: line style for forecast plots under this alternative
  policy. See options from `Plots.jl`. Defaults to `:solid`.

"""
mutable struct MultiPeriodAltPolicy{S <: AbstractDict{Int64, EqcondEntry}} <: AbstractAltPolicy
    key::Symbol
    n_regimes::Int64
    regime_eqcond_info::S
    gensys2::Bool
    temporary_altpolicy_names::Union{Vector{Symbol}, Nothing}
    temporary_altpolicy_length::Int64
    infoset::Union{Vector{UnitRange{Int64}}, Nothing}
    identical_eqcond_regimes::Union{Dict{Int64, Int64}, Nothing}
    perfect_credibility_identical_transitions::Union{Dict{Int64, Int64}, Nothing}
    forecast_init::Function
    color::Colorant
    linestyle::Symbol
end

function MultiPeriodAltPolicy(key::Symbol, n_regimes::Int64, regime_eqcond_info::AbstractDict{Int64, EqcondEntry};
                              gensys2::Bool = false,
                              temporary_altpolicy_names::Union{Vector{Symbol}, Nothing} = nothing,
                              temporary_altpolicy_length::Int =  # default assumes the last regime is a
                              length(regime_eqcond_info) - 1, # lift-off regime, and all others are temporary ones.
                              infoset::Union{Vector{UnitRange{Int64}}, Nothing} = nothing,
                              identical_eqcond_regimes::Union{Dict{Int64, Int64}, Nothing} = nothing,
                              perfect_credibility_identical_transitions::Union{Dict{Int64, Int64}, Nothing} = nothing,
                              forecast_init::Function = identity,
                              color::Colorant = RGB(0., 0., 1.),
                              linestyle::Symbol = :solid)

    MultiPeriodAltPolicy(key, n_regimes, regime_eqcond_info, gensys2, temporary_altpolicy_names,
                         temporary_altpolicy_length, infoset, identical_eqcond_regimes,
                         perfect_credibility_identical_transitions, forecast_init, color, linestyle)
end

"""
```
setup_permanent_altpol!(m::AbstractDSGEModel, altpol::AltPolicy;
    cond_type::Symbol = :none, remove_extra_eqcond_entry::Bool = true)
```

sets up `m` to implement a permanent alternative policy via exogenous regime-switching.
This function assumes that the settings `:date_forecast_start` and/or
`:date_conditional_end` are properly set. The latter setting only matters
if the keyword `cond_type` is set to `:semi` or `:full`.

The keyword `remove_extra_eqcond_entry` removes any `EqcondEntry`
entries in the setting `:regime_eqcond_info` after the
regime implementing the alternative policy, i.e. regimes in the forecast horizon.
If `false`, we will perform this step.
"""
function setup_permanent_altpol!(m::AbstractDSGEModel, altpol::AltPolicy;
                                 cond_type::Symbol = :none, remove_extra_eqcond_entry::Bool = true)

    # Figure out the date of the regime switch to the alt policy
    altpol_date = cond_type == :none ? date_forecast_start(m) : iterate_quarters(date_conditional_end(m), 1)

    # Add regime-switching setting
    if !haskey(get_settings(m), :regime_switching) ||
        (haskey(get_settings(m), :regime_switching) && !get_setting(m, :regime_switching))
        @warn "The setting :regime_switching is either missing or false. Changing it to true"
        m <= Setting(:regime_switching, true)
    end

    # Update regime dates
    if haskey(get_settings(m), :regime_dates)
        n_regs      = length(get_setting(m, :regime_dates)) # calculate length here rather than n_regimes in case
                                                            # setup_regime_switching_inds! has not been called already
        counter     = 0
        reg_missing = true                                  # if false, then altpol date is already in regime_dates
        reg_dates   = get_setting(m, :regime_dates)
        for reg in 1:n_regs
            if reg_dates[reg] > altpol_date
                break
            elseif reg_dates[reg] == altpol_date
                reg_missing = false
                break
            end
            counter += 1
        end
        if counter == n_regs # Case where all regimes occur in history
            altpol_reg = n_regs + 1
            reg_dates[altpol_reg] = altpol_date
        elseif reg_missing
            altpol_reg = counter + 1 # add 1 to get the regime for the alternative policy
            for reg in altpol_reg:(n_regs - 1)
                reg_dates[reg + 1] = reg_dates[reg] # shift all dates up by one
            end
            reg_dates[n_regs + 1] = reg_dates[n_regs] # add extra regime for new date
            reg_dates[altpol_reg] = altpol_date
        else
            altpol_reg = counter + 1 # just need to set the regime number
        end
    else
        altpol_reg = 2
        m <= Setting(:regime_dates, Dict{Int64, Date}(1 => date_presample_start(m), altpol_reg => altpol_date))
    end

    # Call setup_regime_switching_inds! to update indices
    setup_regime_switching_inds!(m; cond_type = cond_type)

    # Update regime_eqcond_info
    m <= Setting(:replace_eqcond, true)
    if haskey(get_settings(m), :regime_eqcond_info)
        reg_eqcond_info = get_setting(m, :regime_eqcond_info)
        reg_eqcond_info[altpol_reg] = EqcondEntry(altpol)
        if remove_extra_eqcond_entry # Remove any extra EqcondEntry values
            for reg in (altpol_reg + 1):get_setting(m, :n_regimes) # go to n_regimes b/c called setup_regime_switching_inds! already
                delete!(reg_eqcond_info, reg)
            end
        end
    else
        m <= Setting(:regime_eqcond_info, Dict{Int, EqcondEntry}(altpol_reg => EqcondEntry(altpol)))
    end

    return m
end