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

using Markdown
using InteractiveUtils

# ╔═╡ e9fdca3a-2d74-4a42-acda-2b289172a1d4
begin
	import Pkg
	Pkg.activate(tempname())
	Pkg.add("DynamicPolynomials")
	Pkg.add("NLopt")
	Pkg.add("StaticArrays")
	Pkg.add("GLMakie")
	Pkg.add("Makie")
	Pkg.add("Combinatorics")
end

# ╔═╡ 9956f874-e481-4e3e-97fc-150aef8b69fa
begin
	using DynamicPolynomials
	using BenchmarkTools
	import NLopt
	using StaticArrays
	using Makie, GLMakie
	using Combinatorics
end

# ╔═╡ 0bf0c2a6-7487-49a9-8adf-73381ea2e9fc
md"""
We denote by ``Π_n^d`` the space of ``n``-variate Polyoniams of degree at most ``d``.
To construct polynomials we use `DynamicPolynomials.jl`.
"""

# ╔═╡ af1b5146-5b3a-425d-95a7-6d0532bcca2f
md"""
The canonical basis is obtained by calculating the non-negative integral solutions to the euqation
```math
x_1 + … + x_n \le d.
```
"""

# ╔═╡ 6c51dda0-5c69-4f58-8348-d62245b1eaeb
md"""
These solutions can be found using the `Combinatorics` package via `multiexponents(n,d)` (`d` must be successively increased).
"""

# ╔═╡ 802e1cd2-fceb-4ee7-b178-62b05ced5409
function non_negative_ineq_solutions(deg, n_vars)
	collect(Iterators.flatten( ( collect(multiexponents( n_vars, d )) for d = 0 : deg ) ))
end

# ╔═╡ c94bdc4e-2c95-42ae-a7dd-4172da278585
function get_poly_basis( deg, n_vars)
	exponents = non_negative_ineq_solutions(deg, n_vars )
	polys = let
		@polyvar x[1:n_vars]
		[ prod(x.^e) for e in exponents ]
	end
	return polys
end

# ╔═╡ b70e7d30-6856-4011-8a59-4727c690634f
md"""
We are going to use the canonical basis to determine a **poised set** of points.
This does in fact work with any polynomial basis for ``Π_n^d``. \
In the process of doing so, we also modify (a copy of?) the basis so that it becomes the Lagrange basis for the returned point set.
"""

# ╔═╡ 42ac42a6-d11e-4c03-bb18-0276d6a319fa
md"The Larange basis is formed by normalizing and orthogonalizing with respect to the point set:"

# ╔═╡ 688d4cf3-4e4d-45de-a43c-141646f01b00
function orthogonalize_polys( poly_arr, x, i )
	# normalize i-th polynomial with respect to `x`
	p_i = poly_arr[i] / poly_arr[i](x)

	# orthogoalize
	return [ j != i ? poly_arr[j] - poly_arr[j](x) * p_i : p_i for j = eachindex( poly_arr ) ]
end

# ╔═╡ e8136cef-9235-4ad8-bfac-a9e176923de7
md"""
We use Algorithm 6.2 and Algorithm 6.3 from the book
"Introduction to Derivative-Free Optimization"
by Conn et. al. \
Algorithm 6.2 makes the set poised (suited for interpolation) and returns the corresponding Lagrange basis.
Algorithm 6.3 takes the poised set and the Lagrange basis and tries to make it ``Λ``-poised. ``Λ`` must be greater 1 and a smaller value makes the set more suited for good models.
"""

# ╔═╡ f48b8d93-e98c-4f30-960c-203a8077b6b6
"""
    get_poised_set( basis, points; solver = :LN_BOBYQA, max_solver_evals = -1 )

Compute a point set suited for polynomial interpolation.

Input:
* `basis`: A vector of polynomials constituting a basis for the polynomial space.
* `points`: (optional) A set of candidate points to be tried for inclusion into the poised set.
* `solver`: NLopt solver to use. Should be derivative-free.
* `max_solver_evals`: Maximum number of evaluations in each optimization run.

Return:
* `poised_points :: Vector{T}` where `T` is either a `Vector{F}` or an `SVector{n_vars, F}` and `F` is the precision of the points in `points`, but at least `Float32`.
* `lagrange_basis :: Vector{<:AbstractPolynomialLike}`: The Lagrange basis corresponding to `poised_points`.
* `point_indices`: An array indicating which points from `points` are also in `poised_points`. A positive entry corresponds to the index of a poised point in `points`. If a poised point is new, then the entry is `-1`.
"""
function get_poised_set( basis, points :: AbstractArray{T} = Vector{Float32}[];
		solver = :LN_BOBYQA, max_solver_evals = -1 ) where {
		T <: AbstractArray{<:Real}
	}

	p = length(basis)
	@assert p > 0 "`basis` must not be an empty array."

	vars = variables( basis[end] )
	n_vars = length(vars)
	@assert n_vars > 0 "The number of variables must be positive."

	if max_solver_evals < 0
		max_solver_evals = 2000 * n_vars
	end

	F = promote_type( eltype( T ), Float32 )
	P_type = n_vars > 100 ? Vector{Vector{F}} : Vector{SVector{n_vars, F}}
	ZERO_TOL = min(eps(F) * 100, eps(Float16) * 10)

	# indicates which points from points have been accepted
	point_indices = fill(-1, p)
	not_accepted_indices = collect( eachindex( points ) )
	# return array of points that form a poised set
	poised_points = P_type(undef, p)

	new_basis = basis
	for i = 1 : p
		_points = points[not_accepted_indices]

		# find the point that maximizes the i-th polynomial
		# if the polynomial is constant, then the first remaining point is used (j = 1)
		l_max, j = if isempty(_points)
			0.0, 0
		else
			findmax( abs.( [ new_basis[i]( x ) for x in _points ] ) )
		end

		if l_max > ZERO_TOL
			# accept the `j`-th point from `_points`
			poised_points[i] = _points[j]
			## indicate what the actual point index was
			point_indices[i] = not_accepted_indices[j]
			## delete from further consideration
			deleteat!(not_accepted_indices, j)
		else
			# no point was suitable to add to the set
			# trying to find the maximizer for a | l_i(x) |
			opt = NLopt.Opt( solver, n_vars )
			opt.lower_bounds = zeros(F, n_vars )
            opt.upper_bounds = ones(F, n_vars )
            opt.maxeval = max_solver_evals
            opt.xtol_rel = 1e-3
            opt.max_objective = (x,g) -> abs( new_basis[i](x) )

			x₀_tmp = [ rand(F, n_vars) for i = 1 : 50 * n_vars ]
            x₀ = x₀_tmp[argmax( abs.(new_basis[i].(x₀_tmp)) ) ]

			_, ξ, ret = NLopt.optimize(opt, x₀)

			poised_points[i] = ξ
		end

		new_basis = orthogonalize_polys( new_basis, poised_points[i], i )
	end

	return poised_points, new_basis, point_indices
end

# ╔═╡ 85d345a0-2d3f-46ed-8b39-f88f700d545b
"""
    make_set_lambda_poised( basis, points;
        LAMBDA = 1.5, solver = :LN_BOBYQA, max_solver_evals = -1, max_loops = -1, skip_indices = [1,] )

Make the output of `get_poised_set` even better suited for interpolation.

Input:
* `basis`: A vector of polynomials constituting a Lagrange basis for the polynomial space.
* `points`: The vector of points belonging to the Lagrange basis.
* `LAMBDA :: Real > 1`: Determines the quality of the interpolation.
* `solver`: NLopt solver to use. Should be derivative-free.
* `max_solver_evals`: Maximum number of evaluations in each optimization run.
* `max_loops`: Maximum number of loops that try to make the set Λ-poised.
* `skip_indices`: Inidices of points to discard last.

Return:
* `poised_points :: Vector{T}` where `T` is either a `Vector{F}` or an `SVector{n_vars, F}` and `F` is the precision of the points in `points`, but at least `Float32`.
* `lagrange_basis :: Vector{<:AbstractPolynomialLike}`: The Lagrange basis corresponding to `poised_points`.
* `point_indices`: An array indicating which points from `points` are also in `poised_points`. A positive entry corresponds to the index of a poised point in `points`. If a poised point is new, then the entry is `-1`.
"""
function make_set_lambda_poised( basis, points :: AbstractArray{T};
		LAMBDA :: Real = 1.5, solver = :LN_BOBYQA, max_solver_evals = -1,
		max_loops = -1, skip_indices = [1,] ) where {
		T <: AbstractArray{<:Real}
	}

	@assert length(basis) == length(points) "Polynomial array `basis` and point array `points` must have the same length."

	if length(points) > 0
		n_vars = length(points[1])
		@assert n_vars > 0 "The number of variables must be positive."

		F = promote_type( eltype( T ), Float32 )
		P_type = n_vars > 100 ? Vector{Vector{F}} : Vector{SVector{n_vars, F}}

		if max_loops < 0
			max_loops = length(basis) * 100
		end

		if max_solver_evals < 0
			max_solver_evals = 2000 * n_vars
		end

		new_basis = basis
		new_points = P_type(points)
		point_indices = collect(eachindex(new_points))

		for k = 1 : max_loops
			iₖ = -1
			xₖ = points[1]
			for (i, polyᵢ) in enumerate(new_basis)
				opt = NLopt.Opt( solver, n_vars )
				opt.lower_bounds = zeros(F, n_vars)
				opt.upper_bounds = ones(F, n_vars)
				opt.maxeval = max_solver_evals
				opt.xtol_rel = 1e-3
				opt.max_objective = (x,g) -> abs( polyᵢ( x ) )

				x₀_tmp = [ rand(F, n_vars) for i = 1 : 50 * n_vars ]
				x₀ = x₀_tmp[argmax( abs.(new_basis[i].(x₀_tmp)) ) ]

				abs_lᵢ, xᵢ, _ = NLopt.optimize(opt, x₀)

				if abs_lᵢ > LAMBDA
					iₖ = i
					xₖ = xᵢ
					if i ∉ skip_indices
						# i is not prioritized we can brake here
						break
					end#if
				end#if
			end#for

			if iₖ > 0
				# perform a point swap
				new_points[iₖ] = xₖ
				point_indices[iₖ] = -1
				# adapt coefficients of lagrange basis
				new_basis = orthogonalize_polys( new_basis, xₖ, iₖ )
			else
				# we are done, the set is lambda poised
				break
			end#if
		end#for

		return new_points, new_basis, point_indices
	else
		return points, basis, collect(eachindex(points))
	end

end

# ╔═╡ c417a8e4-2564-451a-8105-58ae5f713720
md"For the sake of convenience we combine both functions:"

# ╔═╡ a57e64c4-cd79-485a-b21d-97591784cd1c
function get_lambda_poised_set( basis, points; solver1 = :LN_BOBYQA, solver2 = :LN_BOBYQA, max_solver_evals1 = -1, max_solver_evals2 = -1, LAMBDA = 1.5, max_lambda_loops = -1 )
	lagrange_points, lagrange_basis, lagrange_indices = get_poised_set(
		basis, points; solver = solver1, max_solver_evals = max_solver_evals1 )
	lambda_points, lambda_basis, lambda_indices = make_set_lambda_poised(
		lagrange_basis, lagrange_points; LAMBDA, max_loops = max_lambda_loops,
		solver = solver2, max_solver_evals = max_solver_evals2 )
	combined_indices = [ i < 0 ? i : lagrange_indices[j] for (j,i) in enumerate( lambda_indices ) ]
	return lambda_points, lambda_basis, combined_indices
end

# ╔═╡ ffc8cd79-9d85-4f04-81e1-3ae0193395b9
md"Let's have a look at what the points look like:"

# ╔═╡ f72d72de-119b-4c4c-a0d7-38b8ac76c2d6
begin
	basis = get_poly_basis(2,2)
	custom_points =  [ ones(Float32, 2), ones(Float32,2)]

	#lambda_points, lambda_basis, c_indices = get_lambda_poised_set( basis,custom_points)
	lambda_points, lambda_basis, c_indices = get_poised_set( basis,custom_points)
	c_indices
end

# ╔═╡ 9659c50c-c30f-4e78-a80a-eaee23209a4d
scatter(Tuple.(lambda_points))

# ╔═╡ Cell order:
# ╠═e9fdca3a-2d74-4a42-acda-2b289172a1d4
# ╠═9956f874-e481-4e3e-97fc-150aef8b69fa
# ╟─0bf0c2a6-7487-49a9-8adf-73381ea2e9fc
# ╠═af1b5146-5b3a-425d-95a7-6d0532bcca2f
# ╠═6c51dda0-5c69-4f58-8348-d62245b1eaeb
# ╠═802e1cd2-fceb-4ee7-b178-62b05ced5409
# ╠═c94bdc4e-2c95-42ae-a7dd-4172da278585
# ╠═b70e7d30-6856-4011-8a59-4727c690634f
# ╟─42ac42a6-d11e-4c03-bb18-0276d6a319fa
# ╠═688d4cf3-4e4d-45de-a43c-141646f01b00
# ╠═e8136cef-9235-4ad8-bfac-a9e176923de7
# ╠═f48b8d93-e98c-4f30-960c-203a8077b6b6
# ╠═85d345a0-2d3f-46ed-8b39-f88f700d545b
# ╟─c417a8e4-2564-451a-8105-58ae5f713720
# ╠═a57e64c4-cd79-485a-b21d-97591784cd1c
# ╟─ffc8cd79-9d85-4f04-81e1-3ae0193395b9
# ╠═f72d72de-119b-4c4c-a0d7-38b8ac76c2d6
# ╠═9659c50c-c30f-4e78-a80a-eaee23209a4d


text:
using CorpusLoaders
using CorpusLoaders: @NestedVector
using Test
using Base.Iterators
using MultiResolutionIterators
using DataDeps

@testset "basic use $category" for category in ["train_pos", "train_neg", "test_pos", "test_neg"]
    imdb_gen = load(IMDB(category))
    docs = collect(take(imdb_gen, 5))

    @test typeof(docs) == @NestedVector(String, 3)

    words = flatten_levels(docs, (!lvls)(IMDB, :words))|>full_consolidate
    @test typeof(words) == Vector{String}

    sentences = flatten_levels(docs, lvls(IMDB, :documents))|>full_consolidate
    @test typeof(sentences) == @NestedVector(String, 2)
    @test sum(length.(sentences)) == length(words)

    @test typeof(flatten_levels(docs, (!lvls)(IMDB, :documents))|>full_consolidate) == Vector{String}
    @test typeof(flatten_levels(docs, (!lvls)(IMDB, :sentences))|>full_consolidate) == Vector{String}

end

ext: jl


text:
# File to be loaded into eFEMpart

"""
FEMForwardEuler()

Executes Forward Euler time stepping to solve the linear equation Mx=b
and also computes the mass matrix M and b results from a matrix-vector product
using the linear operator. Uses LU factorization on the mass matrix
for the first step, then each step after that costs only O(n^2)

  Right now written to accept :Poisson2D OperatorType
"""
function FEMForwardEuler(prob,mesh,param,U0,Tfinal,dt;
                         printVTK=false,printTotal=100)
  t            = 0.0
  Ndiv = ceil(Tfinal/dt/printTotal)

  # create mass matrix
  M = MassMatrix(mesh)

  # perform LU decomp on mass matrix
  Mlu = lufact(M)

  # obtain weak form of initial condition
  U = forceassemble(mesh,U0)
  sol = ScalarSolution(U)

  # generate derivative matrices
  LinOp = MatrixGen(mesh,prob,param)

  # Apply boundary conditions
  BCApply!(LinOp,mesh,prob,param,prob.OperatorType)

  A = LinOp.Op
  F = LinOp.rhs

  i = 1; j=1
  # step in time
  while t<Tfinal
    # compute RHS
    RHS = (M-dt*A)*U + dt*F
    # use LU factorization to compute solution
    U = Mlu\RHS

    # print problem
    if printVTK
      if mod(j,Ndiv)==0
        s = string(prob.name)
        fn = Path(string(s,"_",i))
        sol.u = U
        vtksave(mesh,prob,param,sol,fn)
        i += 1
      end
    end
    # step time
    t += dt
    j += 1
  end

  return U
end

"""
  progressBar(t,timer1,timer2,timeAvg,counter,
              printTotal,printCounter,printSkip)

Calculates and displays progress bar information
"""
function progressBar(initial,printTotal,printCounter,printSkip)
  now = time()
  timer = now - initial
  percentage = Int(floor(printCounter/printTotal*100))

  progressBarDisplay(percentage,printCounter,printTotal,timer)

  return timer
end

"""
  progressBarDisplay()

Displays progress bar
"""
function progressBarDisplay(percentage,printCounter,printTotal,timer)
  TERMINALLENGTH = 80

  percent = string("\r  [",percentage,"%]")
  skip1  = 14-length(percent)
  exportFigure = string("exported figure ",printCounter,"/",printTotal)
  skip2  = 42-(length(percent)+skip1+length(exportFigure))
  timeTaken = string("stopwatch: ",sToHMS(timer),"\r")

  iterChar(TERMINALLENGTH," ")

  print(percent)

  space()
  iterChar(skip1,"-")
  space()
  print(exportFigure)

  space()
  iterChar(skip2,"-")
  space()
  print(timeTaken)

  if printCounter==printTotal
    println()
  end
end

"""
  iterChar(iterations::Int,character::String)

Prints string character iterations number of times. Used in progress bar
"""
function iterChar(iterations::Int,character::String)
  for i=1:iterations
    print(character)
  end
end

"""
  space()

Prints single space
"""
function space()
  print(" ")
end

"""
  scalarIC(xy,dNodes,dBCarr=0)

creates initial condition array for scalar equation using dirichlet Boundary nodes. If no dirichlet Boundary nodes are given, just returns zero vector
"""
function scalarIC(xy,dNodes,dBCarr)
  N = length(xy)
  ICarr = zeros(N)
  if length(dNodes) != 0
    ICarr[dNodes] = dBCarr
  end

  return ICarr
end

scalarIC(xy) = zeros(Float64,length(xy))

"""
  sToHMS(s::Float64)

Converts Float64 in seconds to HH:MM:SS string format
"""
function sToHMS(s)
  minutes = 0.0
  hours   = 0.0
  timestr = ""

  # compute times
  hours   = Int(floor(s/3600))
  minutes = Int(floor(s/60 - 60*hours))
  seconds = Int(round(s - 3600*hours - 60*minutes,digits=0))

  hours < 10   ? hstr = "0" : hstr = ""
  minutes < 10 ? mstr = "0" : mstr = ""
  seconds < 10 ? sstr = "0" : sstr = ""

  if hours > 0
    timestr=string(timestr,hstr,hours,"h:")
  end
  if (minutes>0) || (hours > 0)
    timestr=string(timestr,mstr,minutes,"m:")
  end

  timestr = string(timestr,sstr,seconds,"s")

  return timestr
end


text:
"""
State machine for these tests:

@startuml
state TransitionTestsMachine {
    state Start
    state Dest

    Start : entry / set entered
    Start : exit / set exited
    Dest : entry / set entered
    Dest : exit / set exited
}
state Unreachable

TransitionTestsMachine : entry / set entered
TransitionTestsMachine : exit / set exited
Unreachable : entry / set entered
Unreachable : exit / set exited
@enduml
"""

using Logging, Test, HSM
include("state_transition_test_utils.jl")

@transition_tests_state TransitionTestsMachine
@transition_tests_state TransitionTestsStateStart
@transition_tests_state TransitionTestsStateDest
@transition_tests_state TransitionTestsStateUnreachable

function HSM.on_entry!(state::TransitionTestsMachine)
    state.entered = true
end

function HSM.on_exit!(state::TransitionTestsMachine)
    state.exited = true
end

function HSM.on_entry!(state::TransitionTestsStateStart)
    state.entered = true
end

function HSM.on_exit!(state::TransitionTestsStateStart)
    state.exited = true
end

function HSM.on_entry!(state::TransitionTestsStateDest)
    state.entered = true
end

function HSM.on_exit!(state::TransitionTestsStateDest)
    state.exited = true
end

function HSM.on_entry!(state::TransitionTestsStateUnreachable)
    state.entered = true
end

function HSM.on_exit!(state::TransitionTestsStateUnreachable)
    state.exited = true
end

function build_state_machine()
    machine = TransitionTestsMachine(nothing)
    start_state = TransitionTestsStateStart(machine)
    dest_state = TransitionTestsStateDest(machine)
    unreachable_state = TransitionTestsStateUnreachable(nothing)

    HSM.transition_to_state!(machine, start_state)

    reset_transition_test_state!(machine)
    reset_transition_test_state!(start_state)
    reset_transition_test_state!(dest_state)
    reset_transition_test_state!(unreachable_state)

    return machine, start_state, dest_state, unreachable_state
end

#
# Transition from Start to Dest should succeed.
#
@testset "transition_to_state!() -- Target state in state machine" begin
    machine, start_state, dest_state, unreachable_state = build_state_machine()
    HSM.transition_to_state!(machine, dest_state)
    @test HSM.active_state(machine) == dest_state
    @test !machine.entered
    @test !machine.exited
    @test !start_state.entered
    @test start_state.exited
    @test dest_state.entered
    @test !dest_state.exited
    @test !unreachable_state.entered
    @test !unreachable_state.exited
end

#
# Transition from Start to Unreachable should fail.
#
@testset "transition_to_state!() -- Target state not in state machine" begin
    machine, start_state, dest_state, unreachable_state = build_state_machine()
    @test_throws HSM.HsmStateTransitionError HSM.transition_to_state!(machine, unreachable_state)
    @test HSM.active_state(machine) == start_state
    @test !machine.entered
    @test !machine.exited
    @test !start_state.entered
    @test !start_state.exited
    @test !dest_state.entered
    @test !dest_state.exited
    @test !unreachable_state.entered
    @test !unreachable_state.exited
end


text:
# load a PBF file into memory
# often scan is a better tool, because it can process larger than memory datasets,
# up to and including Planet.osm

struct OsmData
    ways::Dict{Int64, Way}
    nodes::Dict{Int64, Node}
    relations::Dict{Int64, Relation}
end

OsmData() = OsmData(Dict{Int64, Way}(), Dict{Int64, Node}(), Dict{Int64, Relation}())

function load_pbf(pbffile::AbstractString)
    osm = OsmData()

    scan_pbf(pbffile, nodes=n -> osm.nodes[n.id] = n, ways=w -> osm.ways[w.id] = w, relations=r -> osm.relations[r.id] = r)

    return osm
end


text:
using Clang, ECOS_jll

_GEN_DIR = joinpath(dirname(@__DIR__), "src", "gen")
_COMMON = joinpath(_GEN_DIR, "libecos_common.jl")

function _get_headers(include_path)
    return map(
        f -> joinpath(include_path, f),
        filter(f -> endswith(f, ".h") && f != "timer.h", readdir(include_path)),
    )
end

wc = Clang.init(
    headers = _get_headers(joinpath(ECOS_jll.artifact_dir, "include")),
    output_file = joinpath(_GEN_DIR, "libecos_api.jl"),
    common_file = _COMMON,
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "ecos",
    clang_args = String["-I" * header for header in Clang.find_std_headers()],
    clang_diagnostics = true,
)

run(wc)

function manual_replacements()
    s = read(_COMMON, String)
    for (old, new) in [
        "DBL_MAX + DBL_MAX" => "Inf",
        "ECOS_INFINITY - ECOS_INFINITY" => "NaN",
        "const PRINTTEXT = printf" => "# const PRINTTEXT = printf",
        "const MALLOC = malloc" => "# const MALLOC = malloc",
        "const FREE = free" => "# const FREE = free",
        "const idxint = Cint" => "const idxint = Int64",
    ]
        s = replace(s, old => new)
    end
    write(_COMMON, s)
    return
end

manual_replacements()

rm(joinpath(_GEN_DIR, "LibTemplate.jl"))
rm(joinpath(_GEN_DIR, "ctypes.jl"))


text:
module MT1D
using PyPlot
const μ = 4*pi*1e-7

function Z_f(freqs, ρ, h, irxlayer=1)
    n = length(ρ)
    @assert length(h) == n - 1
    [Z_f(f, ρ, h, n, irxlayer) for f in freqs]
end

function Z_f(f, ρ, h, n, irxlayer)
    # from Kaufman & Keller 1981
    ω = 2pi*f
    k = sqrt(im*ω*μ/ρ[end])
    Z = ω*μ/k
    for i = n-1:-1:irxlayer
        k = sqrt(im*ω*μ/ρ[i])
        ωμ_over_k = ω*μ/k
        Z = ωμ_over_k*coth(-im*k*h[i] + acoth(Z/ωμ_over_k))
    end
    Z
end

ρapp(freqs, Z) = abs2.(Z)./(2pi*freqs*μ)

phase(Z) = rad2deg.(angle.(Z))

function twolayer_ex(h1, ρ1, ρ2::AbstractArray; ntimesperdecade=10)
    X = 10 .^(-1:1/ntimesperdecade:5) # λ₁/h₁
    T = ((h1*X/(2pi*503)).^2)/ρ1
    Z = [twolayer_ex(h1, ρ1, rho2, T) for rho2 in ρ2]
    fig, ax = plt.subplots(2,1, sharex=true, figsize=(6,10))
    for ZZ in Z
        ax[1].loglog(X, ρapp(1 ./T, ZZ)/ρ1)
        ax[2].semilogx(X, phase(ZZ))
    end
    ax[1].set_ylim(1e-3, 1e4)
    ax[2].set_ylim(-95, 10)
    ax[1].grid()
    ax[2].grid()
    ax[2].set_xlabel(L"\lambda_1/h_1")
    ax[1].set_ylabel(L"\rho_a/\rho_1")
    ax[2].set_ylabel("Phase "*L"^\circ")
    fig.suptitle("Kaufman & Keller 1981 Fig 3.1")
    fig.tight_layout()
    nothing
end

function twolayer_ex(h1, ρ1, ρ2::S, T) where S<:Number
    ρ = [ρ1, ρ2]
    Z_f(1 ./T, ρ, [h1])
end

function plotcurve(T, Z; showfreq=false, gridalpha=0.5)
    fig, ax = plt.subplots(1, 2, sharex=true)
    plotcurve(T, Z, fig, showfreq=showfreq, gridalpha=gridalpha)
    fig
end

function plotcurve(T, Z, fig; showfreq=false, iaxis=1, gridalpha=0.5, modelalpha=0.5, lcolor="nocolor")
    f = 1 ./T
    ρₐ = ρapp(f, Z)
    ϕ  = phase(Z)
    ax = fig.axes
    xlabel, abcissa = f_or_T(f, showfreq=showfreq)
    if lcolor == "nocolor"
        ax[iaxis].semilogx(abcissa, log10.(ρₐ))
        ax[iaxis+1].semilogx(abcissa, ϕ)
    else
        ax[iaxis].semilogx(abcissa, log10.(ρₐ), alpha=modelalpha, color=lcolor)
        ax[iaxis+1].semilogx(abcissa, ϕ, alpha=modelalpha, color=lcolor)
    end
    labelaxis(xlabel, ax, iaxis, gridalpha=gridalpha)
    fig.tight_layout()
end

function f_or_T(f; showfreq=false)
    if showfreq
        abcissa = f
        xlabel = "Frequency Hz"
    else
        abcissa = 1 ./f
        xlabel = "Time s"
    end
    xlabel, abcissa
end

function labelaxis(xlabel, ax, iaxis; gridalpha=0.5)
    ax[iaxis].set_xlabel(xlabel)
    ax[iaxis].set_ylabel(L"\log_{10}\rho_{app}"*" (ohm-m)")
    ax[iaxis].grid(b=true, which="both", alpha=gridalpha)
    ax[iaxis+1].set_xlabel(xlabel)
    ax[iaxis+1].set_ylabel("Phase "*L"^\circ")
    ax[iaxis+1].grid(b=true, which="both", alpha=gridalpha)
end

function plotmodelcurve(T, ρ, z; showfreq=false, figsize=(10,4), gridalpha=0.5, logscaledepth=false, irxlayer=1)
    fig = figure(figsize=(figsize))
    s1 = subplot(131)
    s2 = subplot(132)
    s3 = subplot(133, sharex=s2)
    plotmodelcurve(T, ρ, z, fig, showfreq=showfreq, gridalpha=gridalpha, logscaledepth=logscaledepth, irxlayer=irxlayer)
    fig
end

function plotmodelcurve(T, ρ, z, fig; showfreq=false, gridalpha=0.5, logscaledepth=false, lcolor="nocolor", modelalpha=0.5, irxlayer=1)
    f = 1 ./T
    h = diff(z)
    Z = Z_f(f, ρ, h, irxlayer)
    ax = fig.axes
    zlast = diff(z)[end] + z[end]
    if lcolor == "nocolor"
        ax[1].step([ρ; ρ[end]], [z; zlast])
    else
        ax[1].step([ρ; ρ[end]], [z; zlast], color=lcolor, alpha=modelalpha)
    end
    ax[1].set_xlabel("ρ ohm-m")
    ax[1].set_ylabel("Depth m")
    y1, y2 = ax[1].get_ylim()
    y1 < y2 && ax[1].invert_yaxis()
    if logscaledepth
        ymax = max(y1, y2)
        ax[1].set_ylim(ymax, 1)
        ax[1].set_yscale("log")
    end
    ax[1].set_xscale("log")
    ax[1].grid(b=true, which="both", alpha=gridalpha)
    plotcurve(T, Z, fig, showfreq=showfreq, iaxis=2, gridalpha=gridalpha, modelalpha=modelalpha, lcolor=lcolor)
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/misc_utils/period-radius_plots.jl
## (c) 2017 Keir A. Ashby
# Takes any number of SimParams and, calculates observed and physical catalogs with summary statistics and returns:
#   1.A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
#   2.A period-radius distribution of the observed and true catalogs for each sim_param individually
#   3.An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
#   4.An overlaid comparison period-radius distribution of the observed and true catalogs for every input sim_param
#   5.The percentage of observed planets to true planets in each period-radius bin
#   6.CDF plots for radii and periods

using PyPlot
using ExoplanetsSysSim
using Combinatorics

include(joinpath(Pkg.dir(),"ExoplanetsSysSim","src","constants.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2018","abc_setup.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2018","christiansen_func.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2018","misc_utils","calc_sum_stats_pr.jl"))

#see general file definition at the top for a description of this function's purpose.
function sim_param_plots(have_sim_params::Bool; sim_arr::Array{Any,1} = [], num_param::Int64 = 0, print_rates=false) 	#if the user already has SimParams to put in, the first input will be true and the second will be an array (even if there's only one) of your SimParams. If you don't already have one, put false as your first input, and how many SimParams you'd like to generate. If you want to see the ratios of observed to kepler planets in all the bins, set print_rates=true)
  if have_sim_params == false 						#we will make SimParams for the user if she doesn't already have one
    sim_arr = setup_sim_arr(num_param) 					#returns an array of SimParams
  end
  sim_true = setup_sim_arr(1)[1] 					#generates a SimParam that will be used to get the true Kepler catalog.
  ss_arr, ss_true = setup_sum_stats(sim_arr, sim_true) 			#generates physical and observed catalogs and returns summary statistics for each SimParam and for the true Kepler catalog
  num_cats = length(ss_arr)
  n=1									#iterative integer to keep track of figures
  pr_lists = make_pr_lists(num_cats, ss_arr, ss_true)			#Makes lists of all the periods and radii for each catalog
  n = all_cats_pr(n, num_cats, pr_lists)				#A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
  n = all_cats_pr(n, num_cats, pr_lists, false)				#A period-radius distribution of the observed and true catalogs for each sim_param individually
  if num_cats > 1
    n = compare_2_obs_pr(n, num_cats, pr_lists)				#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
  end
  if num_cats > 2
    n = compare_all_obs_pr(n, num_cats, pr_lists)			#An overlaid comparison period-radius distribution of the observed and true catalogs for every input sim_param
  end
  ratios = pr_rates(pr_lists["periods_obs_1"], pr_lists["radii_obs_1"], pr_lists["periods_true"], pr_lists["radii_true"], pr_lists["p_lim_arr"], pr_lists["r_lim_arr"]) #The percentage of observed planets to true planets in each period-radius bin
  n = cdf_plot(n, pr_lists["periods_obs_1"], pr_lists["periods_true"], pr_lists["p_lim_arr"], "Periods") 	#CDF plots for periods
  n = cdf_plot(n, pr_lists["radii_obs_1"], pr_lists["radii_true"], pr_lists["r_lim_arr"], "Radii")		#CDF plots for radii
  return(ratios)							#Normally, pr_rates would just be called last so it printed ratios at the end automatically. However, cdf_plot for radii alters the radii_obs_1 list so that it doesn't work in pr_rates. I'm working on how to alter the list in cdf_plot without altering it globally.
end

function setup_sim_arr(num_param::Int64 = 1) 				#takes an integer and returns an array of that many Sim Params
  sim_arr = Any[]
  for i in 1:num_param
    sim_param = setup_sim_param_christiansen() 				#sim-param used in christiansen_func.jl
    hsu_final_rates = [0.073  0.62   4.9   9.2   11.9 10.  20.  75. 74.
		       0.118  0.376  1.42  5.35  6.3  16.  4.4  8.5 27.
		       0.089  0.236  0.76  2.41  5.04 5.3  7.5  9.4 25.
		       0.052  0.123  0.772 1.84  3.63 3.9  5.3  8.  18.
		       0.053  0.129  0.61  1.24  2.16 3.0  3.8  6.  6.7
		       0.039  0.053  0.402 0.87  1.93 1.6  3.1  3.6 6.
		       0.0105 0.047  0.336 1.16  2.59 4.3  4.8  6.2 13.
   		       0.0058 0.048  0.21  0.96  2.31 2.9  4.1  4.2 5.1
		       0.0036 0.0183 0.133 0.584 1.42 2.1  2.6  4.3 7.4
		       0.0021 0.0167 0.136 0.267 0.48 0.91 1.02 1.5 1.9
		       0.0022 0.0056 0.056 0.101 0.19 0.22 0.74 0.9 2.8
       		       0.0039 0.012  0.086 0.152 0.23 0.36 0.42 1.5 1.7
		       0.0041 0.037  0.119 0.074 0.16 0.19 0.20 0.2 0.69]*0.01	#occurence rates over all bins drawn from results in Hsu et al 2017
    hsu_radii = [0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 6, 8, 12, 16]	#limits for the radius bins
    hsu_periods = [0.5, 1.25, 2.5, 5, 10, 20, 40, 80, 160, 320]			#limits for the period bins
    update_param(sim_param, "obs_par", hsu_final_rates)
    update_param(sim_param, "r_lim_arr", hsu_radii*ExoplanetsSysSim.earth_radius)
    update_param(sim_param, "p_lim_arr", hsu_periods)
    update_param(sim_param, "num_targets_sim_pass_one", 150518)
    update_param(sim_param, "num_kepler_targets", 150518)
    push!(sim_arr,sim_param)
  end
  return(sim_arr)
end

#takes an array of Sim Params and returns summary statistics for each. It will also take one Sim Param and make summary statistics for the true Kepler Catalog
function setup_sum_stats(sim_arr::Array{Any,1}, sim_true::ExoplanetsSysSim.SimParam)
  setup_MES()
  df_star = setup_star_table_christiansen(sim_true)
  df_koi,usable_koi = read_koi_catalog(sim_true)
  cat_true = setup_actual_planet_candidate_catalog(df_star, df_koi, usable_koi, sim_true)
  ss_true = calc_sum_stats_pr(cat_true, sim_true, true)
  ss_arr = Array{Any}(length(sim_arr))
  for i in 1:length(sim_arr)
    cat_phys = generate_kepler_physical_catalog(sim_arr[i])
    cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_arr[i])
    cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys_cut, sim_arr[i])
    ss_arr[i] = calc_sum_stats_pr(cat_obs, sim_arr[i], false, cat_phys)
  end
  return(ss_arr, ss_true)
end

function make_pr_lists(num_cats, ss_arr::Array{Any,1}, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)	#Makes lists of all the periods and radii for each catalog
  p_lim_arr = [0.5, 1.25, 2.5, 5, 10, 20, 40, 80, 160, 320]				#bin boundaries for periods
  r_lim_arr = [0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 6, 8, 12, 16]		#bin boundaries for radii
  pr_lists = Dict("p_lim_arr"=> p_lim_arr, "r_lim_arr"=> r_lim_arr)
  for i in 1:num_cats
    pr_lists["periods_obs_$i"] = ss_arr[i].stat["period_obs_list"] 			#all the periods in this observable catalog
    pr_lists["radii_obs_$i"] = ss_arr[i].stat["radius_obs_list"]./earth_radius 		#all the radii in this observable catalog
    pr_lists["periods_phys_$i"] = ss_arr[i].stat["period_phys_list"] 			#all the periods in this physical catalog
    pr_lists["radii_phys_$i"] = ss_arr[i].stat["radius_phys_list"]./earth_radius 	#all the radii in this physical catalog
    if length(pr_lists["periods_phys_$i"])==0 || length(pr_lists["periods_obs_$i"])==0 	#test to make sure we have something to plot
      println("One of your catalogs is empty")
      return()
    end
  end
  pr_lists["periods_true"] = ss_true.stat["period_list"] 				#all the periods in this true catalog
  pr_lists["radii_true"] = ss_true.stat["radius_list"]./earth_radius 			#all the radii in this true catalog
  return(pr_lists)
end

#A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
function all_cats_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}}, plot_phys = true)
  for i in 1:num_cats
    figure(n)
    if plot_phys == true
      loglog(pr_lists["periods_phys_$i"], pr_lists["radii_phys_$i"], ".", ms = 0.5, color = "orange", label = "Physical")
    end
    loglog(pr_lists["periods_obs_$i"], pr_lists["radii_obs_$i"], ".", ms = 3, label = "Simulated Observed")
    loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
    title("Periods vs. Radii $i", fontsize = 18)
    axis([0.0, 350.0, 0.4, 22.0])
    xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
    yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
    xlabel("Periods (in days)", fontsize=14)
    ylabel("Radius (in Earth radii)", fontsize = 14)
    legend(loc = "lower right", fontsize = 10)
    if plot_phys == true
     savefig("all_cats_pr_$i.pdf")
    else
     savefig("obs_cats_pr_$i.pdf")
    end
    n+=1
  end
  return(n)
end

#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
function compare_2_obs_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}})
  for (a,b) in combinations(collect(1:num_cats),2)
    figure(n)
    loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
    loglog(pr_lists["periods_obs_$a"], pr_lists["radii_obs_$a"], ".", ms = 3, label = "Simulated Observed: $a")
    loglog(pr_lists["periods_obs_$b"], pr_lists["radii_obs_$b"], ".", ms = 3, label = "Simulated Observed: $b")
    title("Periods vs. Radii Comparison: $a & $b", fontsize = 18)
    axis([0.0, 350.0, 0.4, 22.0])
    xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
    yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
    xlabel("Periods (in days)", fontsize = 14)
    ylabel("Radius (in Earth radii)", fontsize = 14)
    legend(loc = "lower right", fontsize = 10)
    n+=1
#    savefig("compare_2_obs_pr_$a_&_$b.pdf")
  end
  return(n)
end

#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
function compare_all_obs_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}})
  figure(n)
  loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
  for i in 1:num_cats
    loglog(pr_lists["periods_obs_$i"], pr_lists["radii_obs_$i"], ".", ms = 3, label = "Simulated Observed: $i")
  end
  title("Periods vs. Radii Comparison: All", fontsize = 18)
  axis([0.0, 350.0, 0.4, 22.0])
  xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
  yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
  xlabel("Periods (in days)", fontsize=14)
  ylabel("Radius (in Earth radii)", fontsize = 14)
  legend(loc = "lower right", fontsize = 10)
#  savefig("compare_all_obs_pr.pdf")
  n+=1
  return(n)
end

#CDF plots for radii and periods
function cdf_plot(n::Int64, obs_list::Array{Float64,1}, true_list::Array{Float64,1}, lim_arr::Array{Float64,1}, cdf_type::String)
  figure(n)
  obs_list = trim_catalog(obs_list, lim_arr[1], lim_arr[length(lim_arr)])	#we want to restrict our catalog to stay within our bin ranges
  true_list = trim_catalog(true_list, lim_arr[1], lim_arr[length(lim_arr)])	#we want to restrict our catalog to stay within our bin ranges
  semilogx(obs_list, collect(linspace(0.0,1.0,length(obs_list))), label = "Simulated Observed")
  semilogx(true_list, collect(linspace(0.0,1.0,length(true_list))), label = "Kepler Observed")
  title("Cumulative Distribution of $cdf_type")
  xlabel("Periods (in days)")
  ylabel("Probability")
  legend(loc = "lower right", fontsize = 10)
#  savefig("cdf_$cdf_type.pdf")
  n+=1
  return(n)
end

#This is a function for testing the accuracy of our model for how many planets Kepler would see. Returns a 2X2 array  of observed/true rates in fraction or decimal form
function pr_rates(periods_obs::Array{Float64,1}, radii_obs::Array{Float64,1}, periods_true::Array{Float64,1}, radii_true::Array{Float64,1}, p_lim_arr::Array{Float64,1}, r_lim_arr::Array{Float64,1}, percent::Bool = true) #if percent is false, the rates will be printed as fractions
  rates_obs = zeros(length(p_lim_arr)-1,length(r_lim_arr)-1)			#table to hold the number of observed planets in each bin
  rates_true = zeros(length(p_lim_arr)-1,length(r_lim_arr)-1)			#table to hold the number of true planets in each bin
  ratios = Array{String}(length(p_lim_arr)-1,length(r_lim_arr)-1)		#table to hold the ratios of observed/true planets in each bin
  for i in 1:length(p_lim_arr)-1						#for each column of periods
    f(p) = p>p_lim_arr[i] && p<p_lim_arr[i+1]
    p_bins_obs = find(f,periods_obs)						#this will hold all the indices of planets that fall into the ith period bin
    r_bins_obs = zeros(length(p_bins_obs))					#this will hold the radius corresponding to each exoplanet whose index is held in p_bins_obs
    for j in 1:length(p_bins_obs)
      r_bins_obs[j]=radii_obs[p_bins_obs][j]
    end
    for k in 1:length(r_lim_arr)-1						#for each row of radii in the ith column of periods
      g(r) = r>r_lim_arr[k] && r<r_lim_arr[k+1]
      rates_obs[i,length(r_lim_arr)-k] = convert(Int64,length(find(g,r_bins_obs)))
    end
  end
#  for h in 1:length(r_lim_arr)-1 						#shows rates_obs
#    println(rates[:,h])
#  end

  for i in 1:length(p_lim_arr)-1						#for each column of periods
    f(p) = p>p_lim_arr[i] && p<p_lim_arr[i+1]
    p_bins_true = find(f,periods_true)						#this will hold all the indices of planets that fall into the ith period bin
    r_bins_true = zeros(length(p_bins_true))					#this will hold the radius corresponding to each exoplanet whose index is held in p_bins_true
    for j in 1:length(p_bins_true)
      r_bins_true[j]=radii_true[p_bins_true][j]
    end
    for k in 1:length(r_lim_arr)-1						#for each row of radii in the ith column of periods
      g(r) = r>r_lim_arr[k] && r<r_lim_arr[k+1]
      rates_true[i,length(r_lim_arr)-k] = convert(Int64,length(find(g,r_bins_true)))
    end
  end
#  for h in 1:length(r_lim_arr)-1						#shows rates_true
#    println(rates_true[:,h])
#  end

  for i in 1:length(p_lim_arr)-1
    for j in 1:length(r_lim_arr)-1
      if percent == true
	if rates_obs[i,j] == 0.0 && rates_true[i,j] == 0.0
	  ratios[i,j] = "NaN"
        elseif rates_true[i,j] == 0.0
	  ratios[i,j] = "Inf"
	else
	  ratios[i,j] = string(floor(Int64,rates_obs[i,j]/rates_true[i,j]*100), "%")
	end
      else    									#returns the ratios in fraction form
        ratios[i,j] = string(convert(Int64,rates_obs[i,j]), "/", convert(Int64,rates_true[i,j]))
      end
    end
  end
  return ratios
end

function trim_catalog(list::Array{Float64,1},min::Float64,max::Float64)		#removes all values in a list to stay within bin ranges
  sort!(list)
  while list[1] < min
    shift!(list)
  end
  while list[length(list)] > max
    pop!(list)
  end
  return(list)
end

function test_period_radius_plots()
  sim_param_plots(false, num_param = 3)
end


text:
module Kernel

using StaticArrays, OffsetArrays
using ..ImageFiltering: centered, KernelFactors
import ..ImageFiltering: _reshape

function product2d(kf)
    k1, k2 = kf
    k1[1].*k1[2], k2[1].*k2[2]
end

"""
`Kernel` is a module implementing filtering kernels of full
dimensionality. The following kernels are supported:

  - `sobel`
  - `prewitt`
  - `ando3`, `ando4`, and `ando5`
  - `gaussian`
  - `DoG` (Difference-of-Gaussian)
  - `LoG` (Laplacian-of-Gaussian)
  - `Laplacian`

See also: KernelFactors.
"""
Kernel

"""
    diff1, diff2 = sobel()

Return kernels for two-dimensional gradient compution using the Sobel
operator. `diff1` computes the gradient along the first (y) dimension,
and `diff2` computes the gradient along the second (x) dimension.

See also: KernelFactors.sobel, Kernel.prewitt, Kernel.ando.
"""
sobel() = product2d(KernelFactors.sobel())

sobel(extended, d) = (broadcast(*, KernelFactors.sobel(extended, d)...),)

"""
    diff1, diff2 = prewitt()

Return kernels for two-dimensional gradient compution using the
Prewitt operator.  `diff1` computes the gradient along the first (y)
dimension, and `diff2` computes the gradient along the second (x)
dimension.

See also: KernelFactors.prewitt, Kernel.sobel, Kernel.ando.
"""
prewitt() = product2d(KernelFactors.prewitt())

prewitt(extended, d) = (broadcast(*, KernelFactors.prewitt(extended, d)...),)

"""
    diff1, diff2 = ando3()

Return 3x3 kernels for two-dimensional gradient compution using the
optimal "Ando" filters.  `diff1` computes the gradient along the
y-axis (first dimension), and `diff2` computes the gradient along the
x-axis (second dimension).

# Citation
Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March
2000

See also: KernelFactors.ando3, Kernel.ando4, Kernel.ando5.
"""
ando3() = product2d(KernelFactors.ando3())

ando3(extended, d) = (broadcast(*, KernelFactors.ando3(extended, d)...),)

"""
    diff1, diff2 = ando4()

Return 4x4 kernels for two-dimensional gradient compution using the
optimal "Ando" filters.  `diff1` computes the gradient along the
y-axis (first dimension), and `diff2` computes the gradient along the
x-axis (second dimension).

# Citation
Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March
2000

See also: KernelFactors.ando4, Kernel.ando3, Kernel.ando5.
"""
function ando4()
    f = centered(@SMatrix [ -0.022116 -0.025526  0.025526  0.022116
                            -0.098381 -0.112984  0.112984  0.098381
                            -0.098381 -0.112984  0.112984  0.098381
                            -0.022116 -0.025526  0.025526  0.022116 ])
    return f', f
end

function ando4(extended::Tuple{Bool,Bool}, d)
    all(extended) || error("all dimensions must be extended")
    (ando4()[d],)
end

"""
    diff1, diff2 = ando5()

Return 5x5 kernels for two-dimensional gradient compution using the
optimal "Ando" filters.  `diff1` computes the gradient along the
y-axis (first dimension), and `diff2` computes the gradient along the
x-axis (second dimension).

# Citation
Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March
2000

See also: `KernelFactors.ando5`, `Kernel.ando3`, `Kernel.ando4`.
"""
function ando5()
    f = centered(@SMatrix [ -0.003776 -0.010199  0.0  0.010199  0.003776
                            -0.026786 -0.070844  0.0  0.070844  0.026786
                            -0.046548 -0.122572  0.0  0.122572  0.046548
                            -0.026786 -0.070844  0.0  0.070844  0.026786
                            -0.003776 -0.010199  0.0  0.010199  0.003776 ])
    return f', f
end

function ando5(extended::Tuple{Bool,Bool}, d)
    all(extended) || error("all dimensions must be extended")
    (ando5()[d],)
end

"""
    gaussian((σ1, σ2, ...), [(l1, l2, ...]) -> g
    gaussian(σ)                  -> g

Construct a multidimensional gaussian filter, with standard deviation
`σd` along dimension `d`. Optionally provide the kernel length `l`,
which must be a tuple of the same length.

If `σ` is supplied as a single number, a symmetric 2d kernel is
constructed.

See also: KernelFactors.gaussian.
"""
@inline gaussian{N}(σs::NTuple{N,Real}, ls::NTuple{N,Integer}) =
    broadcast(*, KernelFactors.gaussian(σs, ls)...)
gaussian(σ::Tuple{Real}, l::Tuple{Integer}) = KernelFactors.gaussian(σ[1], l[1])
gaussian(σ::Tuple{}, l::Tuple{}) = reshape([1])  # 0d
gaussian{T<:Real,I<:Integer}(σs::AbstractVector{T}, ls::AbstractVector{I}) =
    gaussian((σs...,), (ls...,))

@inline gaussian{N}(σs::NTuple{N,Real}) = broadcast(*, KernelFactors.gaussian(σs)...)
gaussian{T<:Real}(σs::AbstractVector{T}) = gaussian((σs...,))
gaussian(σ::Tuple{Real}) = KernelFactors.gaussian(σ[1])
gaussian(σ::Tuple{}) = reshape([1])

gaussian(σ::Real) = gaussian((σ, σ))

"""
    DoG((σp1, σp2, ...), (σm1, σm2, ...), [l1, l2, ...]) -> k
    DoG((σ1, σ2, ...))                                   -> k
    DoG(σ::Real)                                         -> k

Construct a multidimensional difference-of-gaussian kernel `k`, equal
to `gaussian(σp, l)-gaussian(σm, l)`.  When only a single `σ` is
supplied, the default is to choose `σp = σ, σm = √2 σ`. Optionally
provide the kernel length `l`; the default is to extend by two
`max(σp,σm)` in each direction from the center. `l` must be odd.

If `σ` is provided as a single number, a symmetric 2d DoG kernel is
returned.

See also: KernelFactors.IIRGaussian.
"""
DoG{N}(σps::NTuple{N,Real}, σms::NTuple{N,Real}, ls::NTuple{N,Integer}) =
    gaussian(σps, ls) - gaussian(σms, ls)
function DoG{N}(σps::NTuple{N,Real})
    σms = map(s->s*√2, σps)
    neg = gaussian(σms)
    l = map(length, indices(neg))
    gaussian(σps, l) - neg
end
DoG(σ::Real) = DoG((σ,σ))

"""
    LoG((σ1, σ2, ...)) -> k
    LoG(σ)             -> k

Construct a Laplacian-of-Gaussian kernel `k`. `σd` is the gaussian width
along dimension `d`.  If `σ` is supplied as a single number, a
symmetric 2d kernel is returned.

See also: KernelFactors.IIRGaussian and Kernel.Laplacian.
"""
function LoG{N}(σs::NTuple{N})
    w = CartesianIndex(map(n->(ceil(Int,8.5*n)>>1), σs))
    R = CartesianRange(-w, w)
    σ = SVector(σs)
    C = 1/(prod(σ)*(2π)^(N/2))
    σ2 = σ.^2
    σ2i = sum(1./σ2)
    function df(I::CartesianIndex, σ2, σ2i)
        x = SVector(I.I)
        xσ = x.^2./σ2
        (sum(xσ./σ2) - σ2i) * exp(-sum(xσ)/2)
    end
    centered([C*df(I, σ2, σ2i) for I in R])
end
LoG(σ::Real) = LoG((σ,σ))

immutable Laplacian{N}
    flags::NTuple{N,Bool}
    offsets::Vector{CartesianIndex{N}}

    function Laplacian(flags::NTuple{N,Bool})
        offsets = Array{CartesianIndex{N}}(0)
        for i = 1:N
            if flags[i]
                push!(offsets,
                      CartesianIndex{N}((ntuple(d->0, i-1)..., 1, ntuple(d->0, N-i)...)))
            end
        end
        new(flags, offsets)
    end
end

"""
    Laplacian((true,true,false,...))
    Laplacian(dims, N)
    Lacplacian()

Laplacian kernel in `N` dimensions, taking derivatives along the
directions marked as `true` in the supplied tuple. Alternatively, one
can pass `dims`, a listing of the dimensions for
differentiation. (However, this variant is not inferrable.)

`Laplacian()` is the 2d laplacian, equivalent to `Laplacian((true,true))`.

The kernel is represented as an opaque type, but you can use
`convert(AbstractArray, L)` to convert it into array format.
"""
Laplacian{N}(flags::NTuple{N,Bool}) = Laplacian{N}(flags)
Laplacian() = Laplacian((true,true))

function Laplacian(dims, N::Int)
    flags = falses(N)
    flags[[dims...]] = true
    Laplacian((flags...,))
end

Base.indices(L::Laplacian) = map(f->f ? (-1:1) : (0:0), L.flags)
Base.isempty(L::Laplacian) = false
function Base.convert{N}(::Type{AbstractArray}, L::Laplacian{N})
    A = fill!(OffsetArray{Int}(indices(L)), 0)
    for I in L.offsets
        A[I] = A[-I] = 1
    end
    A[ntuple(d->0, Val{N})...] = -2*length(L.offsets)
    A
end
_reshape{N}(L::Laplacian{N}, ::Type{Val{N}}) = L

"""
    reflect(kernel) --> reflectedkernel

Compute the pointwise reflection around 0, 0, ... of the kernel
`kernel`.  Using `imfilter` with a `reflectedkernel` performs convolution,
rather than correlation, with respect to the original `kernel`.
"""
function reflect(kernel::AbstractArray)
    inds = map(-, indices(kernel))
    out = similar(kernel, inds)
    for I in CartesianRange(indices(kernel))
        out[-I] = kernel[I]
    end
    out
end

end


text:
using Pkg

Pkg.activate("$(@__DIR__)/docs")
Pkg.add(PackageSpec(path="$(@__DIR__)"))
Pkg.instantiate()
include("$(@__DIR__)/docs/make.jl")


text:
@testset "Random walks" begin
    """is_nonbacktracking: a predicate that tests if a walk is nonbacktracking.
        That is that for no index i walk[i+2]==walk[i]
    """
    function is_nonbacktracking(walk)
        n = length(walk)
        if n < 3
           return true
        end
        for i in 1:(n - 2)
            if walk[i + 2] == walk[i]
                return false
            end
        end
        return true
    end

    """test_nbw: check that a walk is nonbacktracking and print it if is isn't.
        Used only for testing and debugging.
    """
    function test_nbw(g, start, len)
        w = @inferred(non_backtracking_randomwalk(g, start, len))
        if is_nonbacktracking(w)
            return true
        else
            print("walk was:\n  $w")
        end
        return false
    end
    gx = PathDiGraph(10)
    for g in testdigraphs(gx)
      @test @inferred(randomwalk(g, 1, 5)) == [1:5;]
      @test @inferred(randomwalk(g, 2, 100)) == [2:10;]
      @test_throws BoundsError randomwalk(g, 20, 20)
      @test @inferred(non_backtracking_randomwalk(g, 10, 20)) == [10]
      @test @inferred(non_backtracking_randomwalk(g, 1, 20)) == [1:10;]
    end

    gx = PathGraph(10)
    for g in testgraphs(gx)
      @test @inferred(saw(g, 1, 20)) == [1:10;]
      @test_throws BoundsError saw(g, 20, 20)
      @test @inferred(non_backtracking_randomwalk(g, 1, 20)) == [1:10;]
      @test_throws BoundsError non_backtracking_randomwalk(g, 20, 20)
    end

    gx = SimpleDiGraph(PathGraph(10))
    for g in testdigraphs(gx)
      @test @inferred(non_backtracking_randomwalk(g, 1, 20)) == [1:10;]
      @test_throws BoundsError non_backtracking_randomwalk(g, 20, 20)
    end

    gx = CycleGraph(10)
    for g in testgraphs(gx)
      visited = @inferred(non_backtracking_randomwalk(g, 1, 20))
      @test visited == [1:10; 1:10;] || visited == [1; 10:-1:1; 10:-1:2;]
    end

    gx = CycleDiGraph(10)
    for g in testdigraphs(gx)
      @test @inferred(non_backtracking_randomwalk(g, 1, 20)) == [1:10; 1:10;]
    end

    n = 10
    gx = CycleGraph(n)
    for k = 3:(n - 1)
        add_edge!(gx, 1, k)
    end

    for g in testgraphs(gx)
      for len = 1:(3 * n)
        @test test_nbw(g, 1, len)
        @test test_nbw(g, 2, len)
      end
    end
    #test to make sure it works with self loops.
    add_edge!(gx, 1, 1)
    for g in testgraphs(gx)
      for len = 1:(3 * n)
          @test test_nbw(g, 1, len)
          @test test_nbw(g, 2, len)
      end
    end
end


text:
## Convenient tools
# all-zero dynamics (for test)
function ẋ(env::Fym, x, t)
    println("all-zero dynamics for test; see `LazyFym.ẋ`")
    env_names = names(env)
    if env_names == []
        return zero(x)
    else
        zero_values = env_names |> Map(name -> ẋ(getfield(env, name), x[name], t))
        return (; zip(env_names, zero_values)...)  # NamedTuple
    end
end

# default update (for test)
"""
    update(env::Fym, ẋ, x, t, Δt)

Update law of a simulator.
You may have to customise this function to realise high-level agent
such as zero-order-hold (ZOH) input.
# Examples
function update(rlenv::RLEnv, ẋ, x, t, Δt)
    action = command(rlenv.policy, x.env, t)
    _datum = Dict(:x => x, :t => t, :action => action)
    x_next = ∫(rlenv, ẋ, x, t, Δt; action=action)  # `action` will be a kwarg of ẋ
    _datum[:x_next] = x_next
    datum = (; zip(keys(_datum), values(_datum))...)
    return datum, x_next
end
"""
function update(env::Fym, ẋ, x, t, Δt)  # provided
    _datum = Dict(:x => x, :t => t)
    x_next = ∫(env, ẋ, x, t, Δt)
    _datum[:x_next] = x_next
    datum = (; zip(keys(_datum), values(_datum))...)
    return datum, x_next
end

# automatic completion of initial condition
function initial_condition(env::Fym)
    env_names = names(env)
    values = env_names |> Map(name -> initial_condition(getfield(env, name)))
    return (; zip(env_names, values)...)  # NamedTuple
end


text:
using DSGE, HDF5, ModelConstructors, Test, Dates

regen = false

m = Model1002("ss10"; custom_settings = Dict{Symbol, Setting}(:add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true)))
hist_rule = DSGE.taylor_rule()
m <= Setting(:alternative_policies, AltPolicy[hist_rule])
ngdp_altpol = DSGE.ngdp()
m <= Setting(:pgap_type, :ngdp)
m <= Setting(:pgap_value, 12.)
T_ngdp, R_ngdp, C_ngdp = ngdp_altpol.solve(m)
T_hist, R_hist, C_hist = hist_rule.solve(m)
rp = joinpath(dirname(@__FILE__), "../reference/")

@testset "Gensys when alternative policies are not credible" begin
    m <= Setting(:regime_switching, true)
    setup_permanent_altpol!(m, DSGE.ngdp())
    for (i, prob_vec) in enumerate([[1., 0.], [.5, .5], [.999, .001]])
        get_setting(m, :regime_eqcond_info)[2].weights = prob_vec
        T_prs, R_prs, C_prs = DSGE.gensys_uncertain_altpol(m, prob_vec, AltPolicy[hist_rule];
                                                           regime_switching = true, regimes = [2])
        T_prs, R_prs, C_prs = DSGE.augment_states(m, T_prs, R_prs, C_prs)

        m <= Setting(:uncertain_altpolicy, true)
        TTT, RRR, CCC = solve(m; regime_switching = true, regimes = [2], gensys_regimes = [1:2]) # check automatic calculation
        m <= Setting(:uncertain_altpolicy, false) # to make DSGE.gensys_uncertain_altpol above work

        @test TTT ≈ T_prs
        @test RRR ≈ R_prs
        @test CCC ≈ C_prs

        if regen && i == 2
            h5open(joinpath(rp, "gensys_uncertain_altpol.h5"), "w") do file
                write(file, "T", T_prs)
                write(file, "R", R_prs)
                write(file, "C", C_prs)
            end
        end

        if i == 1
            @test T_prs ≈ T_ngdp
            @test R_prs ≈ R_ngdp
            @test C_prs ≈ C_ngdp
        elseif i == 2
            @test !(T_prs ≈ T_ngdp)
            @test !(R_prs ≈ R_ngdp)
            @test C_prs ≈ C_ngdp    # these should be all zeros
            @test !(T_prs ≈ T_hist)
            @test !(R_prs ≈ R_hist)
            @test C_prs ≈ C_hist    # these should be all zeros
            refpath = joinpath(rp, "gensys_uncertain_altpol.h5")
            @test T_prs ≈ h5read(refpath, "T")
            @test R_prs ≈ h5read(refpath, "R")
            @test C_prs ≈ h5read(refpath, "C")
        else
            @test maximum(abs.(T_prs - T_ngdp)) < 2e-2
            @test maximum(abs.(R_prs - R_ngdp)) < 2e-2
            @test C_prs ≈ C_ngdp    # these should be all zeros
        end
    end

    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), 2 => Date(2020, 9, 30), 3 => Date(2020, 12, 31),
                                     4 => Date(2021, 3, 31), 5 => Date(2021, 6, 30)))
    m <= Setting(:date_forecast_start, Date(2020, 9, 30))
    m <= Setting(:regime_switching, true)
    setup_regime_switching_inds!(m)

    m <= Setting(:alternative_policies, AltPolicy[hist_rule])
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:regime_eqcond_info, Dict(i => DSGE.EqcondEntry(DSGE.ngdp(), [.5, .5]) for i in 2:5))
    m <= Setting(:uncertain_altpolicy, true)
    @test !haskey(DSGE.get_settings(m), :gensys2)
    TTTs, RRRs, CCCs = solve(m; regime_switching = true,
                             gensys_regimes = UnitRange{Int}[1:5],
                             regimes = collect(1:5))

    @test !(TTTs[2] ≈ T_ngdp)
    for i in 3:5
        @test TTTs[i] == TTTs[2]
        @test RRRs[i] == RRRs[2]
        @test CCCs[i] == CCCs[2]
    end

    for i in 2:5
        get_setting(m, :regime_eqcond_info)[i].weights = [1., 0.]
    end
    TTTs2, RRRs2, CCCs2 = solve(m; regime_switching = true,
                                gensys_regimes = UnitRange{Int}[1:5],
                                regimes = collect(1:5))

    @test !(TTTs2[2] ≈ TTTs[2])
    for i in 3:5
        @test TTTs2[i] ≈ T_ngdp
        @test RRRs2[i] ≈ R_ngdp
        @test CCCs2[i] ≈ C_ngdp
    end
end


text:
using DFTK: spglib_spacegroup_number, spglib_standardize_cell
using LinearAlgebra
using Test

@testset "spglib" begin
    a = 10.3
    Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
    Ge = ElementPsp(:Ge, psp=load_psp("hgh/lda/Ge-q4"))

    # silicon
    lattice = a / 2 * [[0 1 1.]; [1 0 1.]; [1 1 0.]]
    atoms     = [Si, Si]
    positions = [ones(3)/8, -ones(3)/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 227
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # silicon with Cartesian x coordinates flipped
    lattice = a / 2 * [[0 -1 -1.]; [1 0 1.]; [1 1 0.]]
    atoms     = [Si, Si]
    positions = [ones(3)/8, -ones(3)/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 227
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # silicon with different lattice vectors
    lattice = a / 2 * [[0 1 1.]; [-1 0 1.]; [-1 1 0.]]
    atoms     = [Si, Si]
    positions = [[-1, 1, 1]/8, -[-1, 1, 1]/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 227
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # Zincblende structure
    lattice = a / 2 * [[0 1 1.]; [1 0 1.]; [1 1 0.]]
    atoms     = [Si, Ge]
    positions = [ones(3)/8, -ones(3)/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 216
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # Zincblende structure with Cartesian x coordinates flipped
    lattice = a / 2 * [[0 -1 -1.]; [1 0 1.]; [1 1 0.]]
    atoms     = [Si, Ge]
    positions = [ones(3)/8, -ones(3)/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 216
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # Zincblende structure with different lattice vectors
    lattice = a / 2 * [[0 1 1.]; [-1 0 1.]; [-1 1 0.]]
    atoms     = [Si, Ge]
    positions = [[-1, 1, 1]/8, -[-1, 1, 1]/8]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 216
    @test spglib_standardize_cell(model).lattice ≈ a * I(3)

    # Hexagonal close packed lattice
    lattice = a * [[1. -1/2 0.]; [0. sqrt(3)/2 0.]; [0. 0. sqrt(8/3)]]
    atoms     = [Si, Si]
    positions = [[0, 0, 0], [1/3, 2/3, 1/2]]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 194
    @test spglib_standardize_cell(model).lattice ≈ lattice

    # Hexagonal close packed lattice with x coordinates flipped
    lattice = a * [[-1. 1/2 0.]; [0. sqrt(3)/2 0.]; [0. 0. sqrt(8/3)]]
    atoms     = [Si, Si]
    positions =  [[0, 0, 0], [1/3, 2/3, 1/2]]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 194
    @test !(spglib_standardize_cell(model).lattice ≈ lattice)

    # Hexagonal close packed lattice with different lattice vectors
    lattice = a * [[-1. -1/2 0.]; [0. sqrt(3)/2 0.]; [0. 0. sqrt(8/3)]]
    atoms     = [Si, Si]
    positions =  [[0, 0, 0], [-1/3, 2/3, 1/2]]
    model = model_LDA(lattice, atoms, positions)
    @test spglib_spacegroup_number(model) == 194
    @test (  spglib_standardize_cell(model).lattice
           ≈ a * [[1. -1/2 0.]; [0. sqrt(3)/2 0.]; [0. 0. sqrt(8/3)]])
end


text:
module pgrad

# this is a port of the example found at
#	http://rl-gym-doc.s3-website-us-west-2.amazonaws.com/mlss/lab1.html#starter

using Reinforce
using OpenAIGym
# using Distributions
using Learn

using MLPlots; gr(size=(1400,1400))

# ----------------------------------------------------------------

# """
# Wraps a LearnBase.Transformation which converts an input vector to action values.
# For discrete actions, it chooses the action which produces the highest value.
# For continuous (interval) actions, it squashes actions to [0,1].
# """
# immutable TransformPolicy{T} <: AbstractPolicy
# 	trans::T
# end
#
# Transformations.params(tp::TransformPolicy) = params(tp.trans)
#
# get_action(A::DiscreteSet, ŷ::AbstractVector) = A[indmax(ŷ)]
# function get_action(A::IntervalSet, ŷ)
#     val = 0.5 * (1.0 + clamp(first(ŷ), -1.0, 1.0))
#     val * (A.hi - A.lo) + A.lo
# end
#
# # # discrete: our action is the action which maximizes the affine transform
# # function Reinforce.action(π::TransformPolicy, r, s, A::DiscreteSet)
# # 	A[indmax(transform!(π.trans, s))]
# # end
#
# # most actionsets pass through to get_action
# Reinforce.action(π::TransformPolicy, r, s, A) = get_action(A)
#
# # get an action where each part of the TupleSet is associated with part of the
# # transformation output
# function Reinforce.action(π::TransformPolicy, r, s, A::TupleSet)
#     ŷ = transform!(π.trans, s)
#     a = []
#     i = 0
#     for Aᵢ in A
#         nᵢ = length(Aᵢ)
#         push!(a, get_action(Aᵢ, view(ŷ, i+1:i+nᵢ)))
#         i += nᵢ
#     end
#     # @show a
#     a
# end
#
# # # continuous: return the transform value, squashed to [0,1]
# # # TODO: remove the squashing when a "Affine + Sigmoid" transformation is available
# # function Reinforce.action(π::TransformPolicy, r, s, A::IntervalSet)
# # 	Transformations.sigmoid(transform(π.trans, s)[1]) * (A.amax-A.amin) + A.amin
# # end

# type Critic
#     return_func::Function
#     baseline_func::Function

# ----------------------------------------------------------------

# initialize a policy, do the learning, then return the policy
function doit(sublearners...; env = GymEnv("BipedalWalker-v2"),
                       maxsteps = 500, # max number of steps in one episode
                       maxiter = 1000, # max learning iterations
					   kw...)
    s = state(env)
    ns = length(s)
    A = actions(env,s)
    nA = length(A)
    @show s ns

    # create a stochastic policy which can sample actions from a multivariate normal dist

    # create a multivariate normal transformation with underlying params μ/σ
    μ = zeros(nA)

    # diagonal
    σ = zeros(nA)
    D = MvNormalTransformation(μ, σ)
    nϕ = 2nA

    # # upper-triangular
    # U = eye(nA,nA)
    # D = MvNormalTransformation(μ, U)
    # nϕ = nA*(nA+1)

    @show D

    # create a neural net mapping: s --> ϕ = vec(μ,U) of the MvNormal
    nh = Int[30,20]
    ϕ = nnet(ns, nϕ, nh, :relu, :identity)
    @show ϕ

    # the critic's value function... mapping state to value
    C = nnet(ns, 1, nh, :relu, :identity)
    @show C

    # our discount rates # TODO: can we learn these too??
    γ = 0.9
    λ = 0.6

    # this is a stochastic policy which follows http://www.breloff.com/DeepRL-OnlineGAE/
    policy = OnlineGAE(A, ϕ, D, C, γ, λ,
                       OnlineGradAvg(400, lr=0.1, pu=Adadelta()),
                       OnlineGradAvg(200, lr=0.1, pu=Adadelta()),
                    #    penalty = ElasticNetPenalty(1e-1,0.5)
                        penalty = L2Penalty(1e-4)
                      )

    # --------------------------------
    # set up the custom visualizations

    # chainplots... put one for each of ϕ/C side-by-side
    cp_ϕ = ChainPlot(ϕ)
    cp_C = ChainPlot(C)

    # this will be called on every timestep of every episode
    function eachiteration(ep,i)
        @show i, ep.total_reward
        update!(cp_ϕ)
        update!(cp_C)
        hm1 = heatmap(reshape(D.dist.μ,nA,1), yflip=true,
                    title=string(maximum(D.dist.μ)),
                    xguide=string(minimum(D.dist.μ)),
                    left_margin=150px)
        # Σ = UpperTriangular(D.dist.Σ.chol.factors)
        # Σ = Diagonal(D.dist.Σ.diag)
        Σ = Diagonal(abs.(output_value(ϕ)[nA+1:end]))
        hm2 = heatmap(Σ, yflip=true,
                    title=string(maximum(Σ)),
                    xguide=string(minimum(Σ)))
        plot(cp_ϕ.plt, cp_C.plt, hm1, hm2, layout = @layout([ϕ; C; hm1{0.2w,0.2h} hm2]))
        # i%2000==0 && gui()
        gui()
    end

    function renderfunc(ep,i)
        if mod1(ep.niter, 10) == 1
            OpenAIGym.render(env, ep.niter, nothing)
        end
    end

    learn!(policy, Episodes(
        env,
        freq = 5,
        episode_strats = [MaxIter(1000)],
        epoch_strats = [MaxIter(5000), IterFunction(renderfunc, every=5)],
        iter_strats = [IterFunction(eachiteration, every=1000)]
    ))

    env, policy
end

end #module


text:
@testset "MOVGP" begin
    N, d = 20, 2
    β = 2.0
    likelihoods = Vector{AGP.AbstractLikelihood}([
        LogisticLikelihood(), LaplaceLikelihood(β)
    ])
    k = SqExponentialKernel() ∘ ScaleTransform(10.0)
    X, f = generate_f(N, d, k)
    X, f2 = generate_f(N, d, k; X=X)
    y_logistic = f .> 0
    y_laplace = f2 + rand(Laplace(β), N)
    ys = [y_logistic, y_laplace]
    floattypes = [Float64]
    model = MOVGP(RowVecs(X), ys, k, likelihoods, AnalyticVI(), 2)

    train!(model, 10)
    predict_y(model, X)
    proba_y(model, X)
end


text:
include("simple_debugger.jl")
@inject function f(x)
    @watch! x
    z = x + 2
    @watch! z
    @watch! z+1
    nested(u) = begin
        @watch! u
        u + x
    end
end
@nodebug!
nested = f(1)
nested(1)

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: thautwarm/MLStyle-Playground

max_forks_repo_name: thautwarm/MLStyle-Playground


text:
import PowerModelsACDC
const _PMACDC = PowerModelsACDC
import PowerModels
const _PM = PowerModels
import Ipopt
import Memento
# import CPLEX
import SCS
import Mosek
import MosekTools
import JuMP
import Gurobi  # needs startvalues for all variables!

file_pf_droop = "./test/data/case5_acdc_droop.m"

file_case5acdc = "./test/data/case5_acdc.m"
file_case5dc ="./test/data/case5_dc.m"
file_case5dcgrid ="./test/data/case5_dcgrid.m"
file_case5dcgridb0 ="./test/data/case5_dcgrid_b0.m"
file_case24 = "./test/data/case24_3zones_acdc.m"
file_case39 = "./test/data/case39_acdc.m"
file_case3120 = "./test/data/case3120sp_acdc.m"
file_case5dcgrid = "./test/data/case5_dcgrid.m"
file_case5_b2bdc = "./test/data/case5_b2bdc.m"
file_lcc = "./test/data/lcc_test.m"
file_588sdet_acdc = "./test/data/pglib_opf_case588_sdet_acdc.m"
file = file_case5acdc

data = _PM.parse_file(file)

_PMACDC.process_additional_data!(data)

ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0)
mosek = JuMP.optimizer_with_attributes(Mosek.Optimizer)
scs = JuMP.optimizer_with_attributes(SCS.Optimizer)
gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer, Presolve=0)
# cplex = JuMP.optimizer_with_attributes(CPLEX.Optimizer)

s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)

resultAC = _PMACDC.run_acdcopf(file, _PM.ACPPowerModel, ipopt; setting = s)

resultLPAC = _PMACDC.run_acdcopf(file, _PM.LPACCPowerModel, ipopt; setting = s)

resultQC = _PMACDC.run_acdcopf(file, _PM.QCRMPowerModel, ipopt; setting = s)

resultSOCBIM = _PMACDC.run_acdcopf(file, _PM.SOCWRPowerModel, ipopt; setting = s)
resultACSOCBIM = _PM.run_opf(file, _PM.SOCWRPowerModel, ipopt; setting = s)
# #
resultSOCBFM = _PMACDC.run_acdcopf_bf(file, _PM.SOCBFPowerModel, ipopt; setting = s)
resultSOCBFMConic = _PMACDC.run_acdcopf_bf(file, _PM.SOCBFConicPowerModel, mosek; setting = s)
resultSOCBFMConicSCS = _PMACDC.run_acdcopf_bf(file, _PM.SOCBFConicPowerModel, scs; setting = s)
# #
resultSDP = _PMACDC.run_acdcopf(file, _PM.SDPWRMPowerModel, mosek; setting = s)
# #
resultDC = _PMACDC.run_acdcopf(file, _PM.DCPPowerModel, gurobi; setting = s)
resultACPF5 = _PMACDC.run_acdcpf(file_case5acdc, _PM.ACPPowerModel, ipopt; setting = s)


text:
# Much of this is copied/adapted from RecipesBase.jl. Cred to everyone who has
# worked on that package!

const _debug_recipes = Bool[false]
function debug(v::Bool = true)
    _debug_recipes[1] = v
end

# check for flags as part of the `-->` expression
function _is_arrow_tuple(expr::Expr)
    expr.head == :tuple && !isempty(expr.args) &&
        isa(expr.args[1], Expr) &&
        expr.args[1].head == :(-->)
end

function _equals_symbol(arg::Symbol, sym::Symbol)
    arg == sym
end
function _equals_symbol(arg::Expr, sym::Symbol) #not sure this method is necessary anymore on 0.7
    arg.head == :quote && arg.args[1] == sym
end
function _equals_symbol(arg::QuoteNode, sym::Symbol)
    arg.value == sym
end
_equals_symbol(x, sym::Symbol) = false

function create_kw_body(func_signature::Expr)
    # get the arg list, stripping out any keyword parameters into a
    # bunch of get!(kw, key, value) lines
    func_signature.head == :where && return create_kw_body(func_signature.args[1])
    args = func_signature.args[2:end]
    kw_body = Expr(:block)
    kw_dict = Dict{Symbol, Any}()
    if isa(args[1], Expr) && args[1].head == :parameters
        for kwpair in args[1].args
            k, v = kwpair.args
            if isa(k, Expr) && k.head == :(::)
                k = k.args[1]
                @warn("Type annotations on keyword arguments not currently supported in recipes. Type information has been discarded")
            end
            push!(kw_body.args, :($k = kwargs[$(Meta.quot(k))]))
            kw_dict[k] = v isa QuoteNode ? v.value : v
        end
        args = args[2:end]
    end
    args, kw_body, kw_dict
end

# build an apply_recipe function header from the recipe function header
function get_function_def(func_signature::Expr, args::Vector)
    front = func_signature.args[1]
    kwarg_expr = Expr(:parameters, Expr(:..., esc(:kwargs)))
    if func_signature.head == :where
        Expr(:where, get_function_def(front, args), esc.(func_signature.args[2:end])...)
    elseif func_signature.head == :call
        #= func = Expr(:call, :(Latexify.apply_recipe), esc.(args)..., Expr(:parameters, :kwargs)) =#
        func = Expr(:call, :(Latexify.apply_recipe), kwarg_expr, esc.(args)...)
        if isa(front, Expr) && front.head == :curly
            Expr(:where, func, esc.(front.args[2:end])...)
        else
            func
        end
    else
        error("Expected `func_signature = ...` with func_signature as a call or where Expr... got: $func_signature")
    end
end

# process the body of the recipe recursively.
# when we see the series macro, we split that block off:
    # let
    #   d2 = copy(d)
    #   <process_recipe_body on d2>
    #   RecipeData(d2, block_return)
    # end
# and we push this block onto the series_blocks list.
# then at the end we push the main body onto the series list
function process_recipe_body!(expr::Expr)
    for (i,e) in enumerate(expr.args)
        if isa(e,Expr)

            # process trailing flags, like:
            #   a --> b, :force
            force = false
            if _is_arrow_tuple(e)
                for flag in e.args
                    if _equals_symbol(flag, :force)
                        force = true
                    end
                end
                e = e.args[1]
            end

            if e.head == :(:=)
                force = true
                e.head = :(-->)
            end

            # we are going to recursively swap out `a --> b, flags...` commands
            # note: this means "x may become 5"
            if e.head == :(-->)
                k, v = e.args
                if isa(k, Symbol)
                    k = QuoteNode(k)
                end

                set_expr = if force
                    :(kwargs[$k] = $v)
                else
                    :(haskey(kwargs, $k) || (kwargs[$k] = $v))
                end

                quiet = false
                expr.args[i] = set_expr

            elseif e.head != :call
                # we want to recursively replace the arrows, but not inside function calls
                # as this might include things like Dict(1=>2)
                process_recipe_body!(e)
            end

            if  e.head == :return
                if e.args[1] isa Expr
                    if e.args[1] isa Tuple
                        e.args[1] = :(($(e.args[1]), kwargs))
                    else
                        e.args[1] = :((($(e.args[1]),), kwargs))
                    end
                else
                    e.args[1] = :((($(e.args[1]),), kwargs))
                end
            end
        end
    end
end

macro latexrecipe(funcexpr)
    func_signature, func_body = funcexpr.args

    if !(funcexpr.head in (:(=), :function))
        error("Must wrap a valid function call!")
    end

    if !(isa(func_signature, Expr) && func_signature.head in (:call, :where))
        error("Expected `func_signature = ...` with func_signature as a call or where Expr... got: $func_signature")
    end

    if length(func_signature.args) < 2
        error("Missing function arguments... need something to dispatch on!")
    end

    args, kw_body, kw_dict = create_kw_body(func_signature)
    func = get_function_def(func_signature, args)
    process_recipe_body!(func_body)

    # now build a function definition for apply_recipe
    funcdef = Expr(:function, func, esc(quote
        if Latexify._debug_recipes[1]
            println("apply_recipe args: ", $args)
        end
        kwargs = merge($kw_dict, kwargs)

        $kw_body
        $func_body
    end))

    return funcdef
end


text:
using MbedTLS: digest, MD_SHA256, MD_MD5

@testset "ASync" begin
    s3region="ap-southeast-2"
    s3url="https://s3.$s3region.amazonaws.com"
    stop_pool_dump=false
    configs = [
        [:verbose => 0],
        [:verbose => 0, :reuse_limit => 200],
        [:verbose => 0, :reuse_limit => 50]
    ]
    protocols = ["http", "https"]

    function dump_async_exception(e, st)
        buf = IOBuffer()
        write(buf, "==========\n@async exception:\n==========\n")
        show(buf, "text/plain", e)
        show(buf, "text/plain", st)
        write(buf, "==========\n\n")
        print(String(take!(buf)))
    end

    function startASyncHTTP()
        @async HTTP.listen() do http
            @show HTTP.Sockets.getsockname(http)
            startwrite(http)
            write(http, """

                    <title>HTTP.jl Connection Pool</title>
                    <meta http-equiv="refresh" content="1">

            """)
            write(http, "
")
            buf = IOBuffer()
            HTTP.ConnectionPool.showpoolhtml(buf)
            write(http, take!(buf))
            write(http, "
")
        end

        @async begin
            sleep(5)
            try
                run(`open http://localhost:8081`)
            catch e
                while !stop_pool_dump
                    HTTP.ConnectionPool.showpool(stdout)
                    sleep(1)
                end
            end
        end
    end

    startASyncHTTP()

    @testset "AWS" begin
        if haskey(ENV, "AWS_ACCESS_KEY_ID") || haskey(ENV, "AWS_DEFAULT_PROFILE")
            @testset "async s3 dup$dup, count$count, sz$sz, pipe$pipe, $protocol, $mode" for
                count in [10, 100, 1000],
                dup in [0, 7],
                protocol in protocols,
                sz in [100, 10000],
                mode in [:request, :open],
                pipe in [0, 32]

                if (dup == 0 || pipe == 0) && count > 100
                    continue
                end # if

                println("running async s3 dup$dup, count$count, sz$sz, pipe$pipe, $protocol, $mode")
                put_data_sums = Dict()
                ch = 100
                conf = [:reuse_limit => 90,
                        :verbose => 0,
                        :pipeline_limit => pipe,
                        :connection_limit => dup + 1,
                        :readtimeout => 120]

                @sync for i = 1:count
                    data = rand(UInt8(65):UInt8(75), sz)
                    md5 = bytes2hex(digest(MD_MD5, data))
                    put_data_sums[i] = md5
                    @async try
                        url = "$s3url/http.jl.test/file$i"
                        r = nothing
                        if mode == :open
                            r = HTTP.open("PUT", url, ["Content-Length" => sz];
                                    body_sha256=digest(MD_SHA256, data),
                                    body_md5=digest(MD_MD5, data),
                                    aws_authorization=true,
                                    conf...) do http
                                for n = 1:ch:sz
                                    write(http, data[n:n+(ch-1)])
                                    sleep(rand(1:10)/1000)
                                end  # for
                            end  # http open
                        end # if
                        if mode == :request
                            r = HTTP.request("PUT", url, [], data;
                                    aws_authorization=true, conf...)
                        end # if
                        @test strip(HTTP.header(r, "ETag"), '"') == md5
                    catch e
                        dump_async_exception(e, stacktrace(catch_backtrace()))
                        rethrow(e)
                    end # try-catch
                end # sync

                get_data_sums = Dict()
                @sync begin
                    for i = 1:count
                        @async try
                            url = "$s3url/http.jl.test/file$i"
                            buf = Base.BufferStream()
                            r = nothing
                            if mode == :open
                                r = HTTP.open("GET", url, ["Content-Length" => 0];
                                              aws_authorization=true,
                                              conf...) do http
                                    buf = Base.BufferStream() # in case of retry!
                                    while !eof(http)
                                        write(buf, readavailable(http))
                                        sleep(rand(1:10)/1000)
                                    end # while
                                    close(buf)
                                end # http open
                            end # if
                            if mode == :request
                                r = HTTP.request("GET", url; response_stream=buf,
                                                             aws_authorization=true, conf...)
                            end # if
                            bytes = read(buf)
                            md5 = bytes2hex(digest(MD_MD5, bytes))
                            get_data_sums[i] = (md5, strip(HTTP.header(r, "ETag"), '"'))
                        catch e
                            dump_async_exception(e, stacktrace(catch_backtrace()))
                            rethrow(e)
                        end # try-catch
                    end # for
                end # sync

                for i = 1:count
                    a, b = get_data_sums[i]
                    @test a == b
                    @test a == put_data_sums[i]
                end # for

                if !haskey(ENV, "HTTP_JL_TEST_LONG_ASYNC")
                    break
                end # if
            end # testset
        end # if haskey(ENV, "AWS_ACCESS_KEY_ID")
    end

    @testset "HTTP.request - Headers - $config - $protocol" for config in configs, protocol in protocols
        result = []

        @sync begin
            for i = 1:100
                @async try
                    response = HTTP.request("GET", "$protocol://httpbin.org/headers", ["i" => i]; config...)
                    response = JSON.parse(String(response.body))
                    push!(result, response["headers"]["I"] => string(i))
                catch e
                    dump_async_exception(e, stacktrace(catch_backtrace()))
                    rethrow(e)
                end
            end
        end

        for(a, b) in result
            @test a == b
        end

        HTTP.ConnectionPool.closeall()
    end

    @testset "HTTP.request - Body - $config - $protocol" for config in configs, protocol in protocols
        result = []

        @sync begin
            for i=1:100
                @async try
                    response = HTTP.request("GET", "$protocol://httpbin.org/stream/$i"; config...)
                    response = String(response.body)
                    response = split(strip(response), "\n")
                    push!(result, length(response) => i)
                catch e
                    dump_async_exception(e, stacktrace(catch_backtrace()))
                    rethrow(e)
                end
            end
        end

        for (a,b) in result
            @test a == b
        end

        HTTP.ConnectionPool.closeall()
    end

    @testset "HTTP.open - $config - $protocol" for config in configs, protocol in protocols
        result = []

        @sync begin
            for i=1:100
                @async try
                    open_response = nothing
                    url = "$protocol://httpbin.org/stream/$i"

                    response = HTTP.open("GET", url; config...) do http
                        open_response = String(read(http))
                    end

                    @test response.status == 200

                    open_response = split(strip(open_response), "\n")
                    push!(result, length(open_response) => i)
                catch e
                    dump_async_exception(e, stacktrace(catch_backtrace()))
                    rethrow(e)
                end
            end
        end

        for (a,b) in result
            @test a == b
        end

        HTTP.ConnectionPool.closeall()
    end

    @testset "HTTP.request - Response Stream - $config - $protocol" for config in configs, protocol in protocols
        result = []

        @sync begin
            for i=1:100
                @async try
                    stream_response = nothing

                    # Note: url for $i will give back $i responses split on "\n"
                    url = "$protocol://httpbin.org/stream/$i"

                    try
                        stream = Base.BufferStream()
                        response = HTTP.request("GET", url; response_stream=stream, config...)

                        @test response.status == 200
                        stream_response = String(read(stream))
                    catch e
                        if !HTTP.RetryRequest.isrecoverable(e)
                            rethrow(e)
                        end
                    end

                    stream_response = split(strip(stream_response), "\n")
                    push!(result, length(stream_response) => i)
                catch e
                    dump_async_exception(e, stacktrace(catch_backtrace()))
                    rethrow(e)
                end
            end
        end

        for (a,b) in result
            @test a == b
        end

        HTTP.ConnectionPool.closeall()
    end

    stop_pool_dump = true  # Kill the HTTP server
end


text:
mutable struct PlateUpdateState{U,V,W}
    score::Float64
    is_empty::Bool
    args::U
    nodes::V
    changes::W
    discard::DynamicAssignment
end

function update_existing_trace!(gen::Plate{T,U}, key::Int, state::PlateUpdateState, subtraces, retvals) where {T,U}
    node = haskey(state.nodes, key) ? state.nodes[key] : EmptyAssignment()
    if haskey(state.changes, key)
        kernel_change = state.changes[key]
    else
        kernel_change = nothing
    end
    subtrace::U = subtraces[key]
    prev_call = get_call_record(subtrace)
    kernel_args = get_args_for_key(state.args, key)
    (subtrace, _, kernel_discard) = update(
        gen.kernel, kernel_args, kernel_change, subtrace, node)
    set_internal_node!(state.discard, key, kernel_discard)
    call::CallRecord = get_call_record(subtrace)
    state.score += (call.score - prev_call.score)
    subtraces = assoc(subtraces, key, subtrace)
    retvals = assoc(retvals, key, call.retval::T)
    state.is_empty = state.is_empty && !has_choices(subtraces)
    (subtraces, retvals)
end

function assess_new_trace!(gen::Plate{T,U}, key::Int, state::PlateUpdateState, subtraces, retvals) where {T,U}
    node = haskey(state.nodes, key) ? state.nodes[key] : EmptyAssignment()
    kernel_args = get_args_for_key(state.args, key)
    subtrace::U = assess(gen.kernel, kernel_args, node)
    subtraces = push(subtraces, subtrace)
    call::CallRecord = get_call_record(subtrace)
    state.score += call.score
    retvals = push(retvals, call.retval::T)
    (subtraces, retvals)
end

function update_process_constraints!(nodes, to_visit, key::Int, node, new_length)
    if key <= new_length
        push!(to_visit, key)
        nodes[key] = node
    else
        error("Update did not consume constraints at key $key")
    end
end

function update_process_constraints!(nodes, to_visit, key, node, new_length)
    # key is something other than an Int
    error("Update did not consume constraints at key $key")
end

function _update(gen::Plate{T,U}, args, prev_trace::VectorTrace{T,U}, constraints, to_visit, changes) where {T,U}

    (new_length, prev_length) = get_prev_and_new_lengths(args, prev_trace)

    # collect constraints, indexed by key
    nodes = Dict{Int, Any}()
    for (key, node) in get_internal_nodes(constraints)
        update_process_constraints!(nodes, to_visit, key, node, new_length)
    end

    # discard subtraces of deleted applications
    discard = DynamicAssignment()
    for key=new_length+1:prev_length
        set_internal_node!(discard, key, get_assignment(prev_trace.subtraces[key]))
    end

    # collect initial state
    state = PlateUpdateState(prev_trace.call.score, prev_trace.is_empty, args, nodes, changes, discard)
    subtraces = prev_trace.subtraces
    retvals = prev_trace.call.retval

    # visit existing applications
    for key in to_visit
        (subtraces, retvals) = update_existing_trace!(gen, key, state, subtraces, retvals)
    end

    # visit each new application and assess the new trace
    for key=prev_length+1:new_length
        (subtraces, retvals) = assess_new_trace!(gen, key, state, subtraces, retvals)
    end

    call = CallRecord(state.score, retvals, args)
    trace = VectorTrace(subtraces, call, state.is_empty)

    # compute the weight
    prev_score = prev_trace.call.score
    weight = state.score - prev_score
    retchange = nothing
    (trace, weight, state.discard, retchange)
end

# No change to arguments
function update(gen::Plate{T,U}, args, change::NoChange, trace::VectorTrace{T,U}, constraints) where {T,U}
    to_visit = Set{Int}()
    changes = Dict{Int, T}()
    _update(gen, args, trace, constraints, to_visit, changes)
end

# Known change to arguments
function update(gen::Plate{T,U}, args, change::PlateChange, trace::VectorTrace{T,U}, constraints) where {T,U}

    (new_length, prev_length) = get_prev_and_new_lengths(args, trace)

    # calculate which existing applications to visit
    to_visit = Set{Int}()
    changes = Dict{Int, T}()
    for (i, kernel_change) in zip(change.changed_args, change.sub_changes)
        if i <= new_length
            push!(to_visit, i)
            changes[i] = kernel_change
        end
    end

    _update(gen, args, trace, constraints, to_visit, changes)
end

# Unknown change to arguments
function update(gen::Plate{T,U}, args, change::Nothing, trace::VectorTrace{T,U}, constraints) where {T,U}

    (new_length, prev_length) = get_prev_and_new_lengths(args, trace)

    # visit all existing applications that are preserved
    to_visit = Set(1:min(prev_length, new_length))
    changes = Dict{Int, Any}() # not used

    _update(gen, args, trace, constraints, to_visit, changes)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
function probability_of_failure(
    models::Union{Array{<:UQModel},UQModel},
    performance::Function,
    inputs::Union{Array{<:UQInput},UQInput},
    sim::AbstractMonteCarlo,
)
    samples = sample(inputs, sim)
    evaluate!(models, samples)

    # Probability of failure
    pf = sum(performance(samples) .< 0) / sim.n

    return pf, samples
end

function probability_of_failure(
    models::Union{Array{<:UQModel},UQModel},
    performance::Function,
    inputs::Union{Array{<:UQInput},UQInput},
    sim::LineSampling,
)
    if isempty(sim.direction)
        sim.direction = gradient_in_standard_normal_space(
            [models..., Model(x -> -1 * performance(x), :performance)],
            inputs,
            mean(inputs),
            :performance,
        )
    end

    samples = sample(inputs, sim)
    evaluate!(models, samples)

    p = reshape(performance(samples), length(sim.points), sim.lines)

    ϕ = Normal()
    pf = 0
    roots_found = 0
    x = median(sim.points)
    for i in 1:(sim.lines)
        if all(p[:, i] .< 0)
            @warn "All samples for line $i are inside the failure domain"
            continue
        elseif all(p[:, i] .> 0)
            @warn "All samples for line $i are outside the failure domain"
            continue
        end
        spl = Spline1D(sim.points, p[:, i])
        try
            root = Dierckx.roots(spl)[1]
            pf += cdf.(ϕ, -1 * root)
            roots_found += 1
        catch e
            @warn "Intersection with failure domain not found for line $i ($e)"
        end
    end

    pf /= roots_found

    return pf, samples
end

# Allow to calculate the pf using only a performance function but no model
function probability_of_failure(
    performance::Function, inputs::Union{Array{<:UQInput},UQInput}, sim::Any
)
    return probability_of_failure(UQModel[], performance, inputs, sim)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
"""
Calculate Penalized LogP by using RDKit package
"""

using PyCall
using Printf

struct RDKit
    molecules::PyObject
end

function init()
    push!(pyimport("sys")["path"], @__DIR__)
    molecules = pyimport("molecules")
    return RDKit(molecules)
end

function main()
    println("hello")
end

"""
    penalized_logp()

`smilesstr` is Smiles string.
return penalized LogP value of `smilesstr` molecule.
"""
function penalized_logp(rdkit::RDKit, smilesstr::String)
    mol = rdkit.molecules.smiles_to_mol(smilesstr)
    return rdkit.molecules.penalized_logp(mol)
end

max_stars_repo_name: Soyukke/MolProperties

max_issues_repo_name: Soyukke/MolProperties

max_forks_repo_name: Soyukke/MolProperties


text:
#!plot # ## Quick Tutorial [![Open Notebook](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)](https://nbviewer.jupyter.org/github/Juice-jl/LogicCircuits.jl/blob/gh-pages/dev/generated/usage.ipynb)

#plot # ### Quick Tutorial
#plot # [![Open Notebook](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)](https://nbviewer.jupyter.org/github/Juice-jl/LogicCircuits.jl/blob/gh-pages/dev/generated/usage.ipynb)

# Assuming that the LogicCircuits Julia package has been installed with `julia -e 'using Pkg; Pkg.add("LogicCircuits")'`, we can start using it as follows.

using LogicCircuits

# ### Reasoning with manually constructed circuits

# We begin by creating three positive literals (logical variables) and manually constructing a simple circuit using logical connectives & (and), | (or), and - (not).
sun, rain, rainbow = pos_literals(LogicCircuit, 3)
circuit = (rainbow & sun & rain) | (-rainbow); # rainbow implies sun and rain

# Just like any logical circuit or Boolean function, we can evaluate ours on various inputs.
circuit(false, true, true) # sun is false, rain is true, rainbow is true
@test false == circuit(false, true, true) #jl
#-
circuit(true, true, true) # sun is true, rain is true, rainbow is true
@test true == circuit(true, true, true) #jl

#src TODO: specific for ipython, visualize circuit using plot(circuit)?

# The purpose of this package, however, is to enable more interesting inference scenarios. This is possible by ensuring that the circuit has certain [properties](https://juice-jl.github.io/LogicCircuits.jl/dev/manual/properties/), such as *decomposability*, *determinism*, and more.
# Our current circuit happens to already be decomposable and deterministic by construction:
isdecomposable(circuit) && isdeterministic(circuit)
@test true == isdecomposable(circuit) && isdeterministic(circuit) #jl

# The decomposability property ensures that we can ask whether the circuit is satisfiable (the classical SAT problem) and, surprisingly, still get our answer efficiently. Of course, from the input `true, true, true` tried above, we know the answer to be true.
issatisfiable(circuit) # does there exist an input that outputs true?
@test true == issatisfiable(circuit) #jl
# In addition, the determinism property allows us to efficiently decide whether the circuit is a tautology (always true), or compute its model count, that is, the number of satisfying assignments.
istautology(circuit) # do all inputs give the circuit output true?
@test false == istautology(circuit) #jl
#-
model_count(circuit) # how many possible inputs give the output true?
@test 5 == model_count(circuit) #jl

# ### Reasoning with compiled circuits

# As logical sentences become more complicated, it becomes infeasible to manually write down circuits that have the requisite properties that guarantee tractable inference.

# A process called *compilation* can solve this problem. Concretely, `LogicCircuits` supports compilation into a particular type of circuit called SDD. We construct an SDD manager with four additional variables, and then ask to compile our running example circuit into an SDD:
manager = SddMgr(7, :balanced)
circuit = compile(manager, circuit);

# Now we are able to incorporate many more logical sentences into the same circuit.
sun, rain, rainbow, cloud, snow, los_angeles, belgium = pos_literals(Sdd, manager, 7)
circuit &= (-los_angeles | -belgium) # cannot be in LA and Belgium at the same time
circuit &= (los_angeles ⇒ sun) ∧ (belgium ⇒ cloud) # unicode logical syntax
circuit &= (¬(rain ∨ snow) ⇐ ¬cloud); # no rain or snow without clouds

# Incorporating these constraints has increased the size of our circuit.

#src when plotting is enabled, just show circuit
#plot plot(circuit; simplify=true)

#src when plotting is not enabled, generate image manually and pretend
#!plot using TikzPictures #hide
#!plot file = ENV["JUICE_MAKE_DOC_SRC_GENERATED"]*"/example-circuit.svg" #hide
#!plot rm(file, force=true) #hide
#!plot save(SVG(file), plot(circuit; simplify=true)) #hide
#!plot # ```julia
#!plot # plot(circuit; simplify=true)
#!plot # ```

#!plot # <img src="https://juice-jl.github.io/LogicCircuits.jl/dev/generated/example-circuit.svg" alt="Example Logic Circuit">

# Crucially, the circuit is still decomposable and deterministic.
isdecomposable(circuit) && isdeterministic(circuit)
@test true == isdecomposable(circuit) && isdeterministic(circuit) #jl
# This means that we can still decide satisfiability, count models, and solve various inference tasks efficiently. For example, we can compute the fraction of inputs that gives the output true:
sat_prob(circuit)
@test 29//128 == sat_prob(circuit) #jl

# Moreover, compiled SDD circuits allow for efficiently checking whether one circuit logically entails another circuit, and whether two circuits are logically equivalent.
entails(circuit, (rainbow ⇒ cloud))
@test true == entails(circuit, (rainbow ⇒ cloud)) #jl
#-
entails(circuit, (rainbow ⇒ belgium))
@test false == entails(circuit, (rainbow ⇒ belgium)) #jl
#-
equivalent((rainbow ⇒ belgium), (¬belgium ⇒ ¬rainbow))
@test true == equivalent((rainbow ⇒ belgium), (¬belgium ⇒ ¬rainbow)) #jl

# Logical constraints are often written in conjunctive normal form (CNF). These can be loaded from file and compiled into circuits, using an SDD manager whose decomposition structure is specified by a *vtree* file.
manager = SddMgr(zoo_vtree("iscas89/s208.1.scan.min.vtree"))
circuit = compile(manager, zoo_cnf("iscas89/s208.1.scan.cnf")) # CNF has 285 clauses
@test 262144 == model_count(circuit) #jl
@test 3115 == num_nodes(circuit) #jl
@test 5826 == num_edges(circuit) #jl
"This CNF has $(model_count(circuit)) satisfying assignments. Its circuit has $(num_nodes(circuit)) nodes and $(num_edges(circuit)) edges."

# ### Advanced functionality

# `LogicCircuits` further provides
#  * CPU (SIMD) and GPU (CUDA) kernels to efficiently compute satisfiability, model counts, etc., for large numbers of inputs, parallelizing over both circuit nodes and data inputs.
#  * Algorithms that transform circuits in non-trivial ways (split, clone, smooth, condition, etc.), verify and enforce structural properties.
#  * Functionality to load and save circuits in various file formats

#src TODO: example of downward pass conditional probability or conditional SAT?
#src issmooth(circuit)
#src example of smoothing a compiled SDD
#src smooth_circuit = smooth(circuit)
#src isdecomposable(smooth_circuit) && isdeterministic(smooth_circuit) && issmooth(smooth_circuit)
#src model_var_prob(smooth_circuit)

#src TODO example of vectorized inference?

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
include("prop_logic_theory.jl")
include("prover.jl")

using Test

Metatheory.options.verbose = true

ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
@profview prove(t, ex, 2, 7)

using Metatheory
using Metatheory.Classic
@metatheory_init ()
using Metatheory.EGraphs
using Metatheory.Library
using Metatheory.Util
using Metatheory.EGraphs.Schedulers

Metatheory.options.printiter = true
Metatheory.options.verbose = true

function rep(x, op, n::Int)
    foldl((x, y) -> :(($op)($x, $y)), repeat([x], n))
end

rep(:a, :*, 3)

Mid = @theory begin
    a * :ε => a
    :ε * a => a
end

Massoc = @theory begin
    a * (b * c) => (a * b) * c
    (a * b) * c => a * (b * c)
end

T = [
    @rule :b*:B => :ε
    @rule :a*:a => :ε
    @rule :b*:b*:b => :ε
    @rule :B * :B => :B
    RewriteRule(Pattern(rep(:(:a*:b), :*, 7)), Pattern(:(:ε)))
    RewriteRule(Pattern(rep(:(:a*:b*:a*:B), :*, 5)), Pattern(:(:ε)))
]

G = Mid∪Massoc∪T
expr = :(a*b* a*a*a * b*b*b * a * B*B*B*B * a)

ex = expr
g = EGraph(expr)
params = SaturationParams(timeout=8, scheduler=BackoffScheduler)#, schedulerparams=(128,4))#, scheduler=SimpleScheduler)
@timev saturate!(g, G, params)
ex = extract!(g, astsize)
println(ex)
rewrite(ex, Mid) |> println

another_expr = :(b*B)
g = EGraph(another_expr)
@timev saturate!(g, G, params)

# another_expr = :(a*a*a*a)
some_eclass = addexpr!(g, another_expr)
# g.root = some_eclass.id
ex = extract!(g, astsize)