text:
# TODO: Add these to the .gitignore of the created project
#    .vscode
#    .vscode/*

# ==================================================================================
#                Custom Templates
# ==================================================================================

# Template: Tests.yml
# This is the template templates\github\workflows\CI.yml of PkgTemplates v0.7.18 with,
# - Remove the job "docs" that build/deploy the doc for Documenter.jl
# - Changed the triggers, we keept the condition [push, pull_request], but only on branch "main"

# Template: CompatHelper.yml
# This is the template templates\github\workflows\CompatHelper.yml of PkgTemplates v0.7.18 with,
# - Remove the authentication line in env: COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}-
#   If the package deploy doc to another repo, it will be a DOCUMENTER_KEY. If a COMPATHELPER_PRIV
#   is present with a non-empty key, then CompatHelper will use it to deploy in preference to
#   the GITHUB_TOKEN. However, the DOCUMENTER_KEY may NOT give access to the repo itself but
#   only to the repo to which the doc is deployed.

# Template: README.md
# This is the template templates\README.md of PkgTemplates v0.7.18 with,
# - Very verbose extension...

# Template: docs_make.jlt
# This is the template templates\docs\make.jl of PkgTemplates v0.7.18 with,
# - Changed extension from .jl to .jlt to avoid spurious IDE errors
# - Changed quite a lot to allow deployment to another repo

# Template: docs_index.md
# This is the template templates\docs\src\index.md of PkgTemplates v0.7.18 with,
# - Currently, nothing changed.

# ==================================================================================
#                Type definitions
# ==================================================================================

using PkgTemplates: @plugin, @with_kw_noshow, Plugin, FilePlugin

@plugin struct GitPagesDoc <: FilePlugin
    privaterepo::Bool = false              # true if this package will upstream-ed to a private GitHub repo
    deploy_gituser::StrOpt = nothing       # IF nothing: deploy doc in the package repo and in gh-pages branch
    deploy_gitrepo::StrOpt = nothing       # deploy doc into repo "github/ylvain/PkDoc"
    deploy_branch::StrOpt = nothing        # deploy doc into branch main of the above repo
    deploy_subdir::StrOpt = nothing        # If nothing deploys at the root ("/") of the target repo filsystem
                                           # If "<PKG>" deploys at "/PkgRepoName" where PkgRepoName is the name of the packge repo
end

# ==================================================================================
#                Basics accessors
# ==================================================================================

function getplugin(::Type{T}, t::Template) where {T<:Plugin}
    idx = findfirst(x -> x isa T, t.plugins)
    isnothing(idx) && error("A $T plugin must be declared in the template")
    return t.plugins[idx] :: T
end

function isprivaterepo(t::Template)
    getplugin(GitPagesDoc, t).privaterepo
end

function package_reponame(t::Template, pkg::AbstractString)
    getplugin(Git, t).jl ? pkg * ".jl" : pkg
end

function userrepo_code(t::Template, pkg::AbstractString)
    t.user * "/" * package_reponame(t, pkg)
end

function hostuserrepo_code(t::Template, pkg::AbstractString)
    t.host * "/" * userrepo_code(t, pkg)
end

function userrepo_docs(t::Template, pkg::AbstractString)
    c = getplugin(GitPagesDoc, t)
    p = package_reponame(t, pkg)

    if isnothing(c.deploy_gituser)
        "$(t.user)/$p"
    else # Deploymant to another repo
        "$(c.deploy_gituser)/$(c.deploy_gitrepo)"
    end
end

function hostuserrepo_docs(t::Template, pkg::AbstractString)
    t.host * "/" * userrepo_docs(t, pkg)
end

function deploy_samerepo(t::Template, pkg::AbstractString)
    userrepo_code(t, pkg) == userrepo_docs(t, pkg)
end

function deploy_docdir(t::Template, pkg::AbstractString)
    c = getplugin(GitPagesDoc, t)
    if isnothing(c.deploy_subdir)
        return nothing
    end

    if c.deploy_subdir == "<PKG>"
        package_reponame(t, pkg)
    else
        c.deploy_subdir
    end
end

function deploy_branch(t::Template, pkg::AbstractString)
    c = getplugin(GitPagesDoc, t)
    if isnothing(c.deploy_branch)
        return deploy_samerepo(t, pkg) ? "gh-pages" : "main"
    end
    return c.deploy_branch
end

function pages_url_base(t::Template, pkg::AbstractString)
    c = getplugin(GitPagesDoc, t)
    p = package_reponame(t, pkg)
    dir_opt = deploy_docdir(t, pkg)
    d = ""
    if !isnothing(dir_opt)
        d = "/" * dir_opt
    end

    if isnothing(c.deploy_gituser)
        "https://$(t.user).github.io/$p" * d
    else
        "https://$(c.deploy_gituser).github.io/$(c.deploy_gitrepo)" * d
    end
end

pages_url_stable(t::Template, pkg::AbstractString) = pages_url_base(t, pkg) * "/stable"
pages_url_dev(t::Template, pkg::AbstractString)    = pages_url_base(t, pkg) * "/dev"

# ==================================================================================
#                Extends sub-environments
# ==================================================================================

const DEP_DOCUMENTER = PackageSpec(;name="Documenter", uuid="e30172f5-a6a5-5a46-863b-614d45cd2de4")
const DEP_TEST       = PackageSpec(;name="Test"      , uuid="8dfed614-e22c-5e08-85e1-65c5234f0b40")

# We cannot extends the way PkgTemplates generates the Test environment without
# creating another full Plugin (that would delete most of the work to the existing
# one). See: PkgTemplates \src\plugins\tests.jl.
# We also avoid to overload "make_test_project" in the above file.
# So we simply add our changes after package creation.
#
# NOTE: pkgconfig_env_test ==> need to commit these

function pkgconfig_env_test(pkg_dir::String)
    PkgTemplates.with_project(joinpath(pkg_dir, "test")) do
        # Pkg.add(DEP_TEST) (done by PkgTemples)
        Pkg.add(DEP_DOCUMENTER)
    end
end

# ==================================================================================
#                PkgTemplates
# ==================================================================================

import PkgTemplates: Badge

PkgTemplates.source(::GitPagesDoc) = template("Documenter.yml")
PkgTemplates.destination(::GitPagesDoc) = joinpath(".github","workflows","Documenter.yml")
PkgTemplates.tags(::GitPagesDoc) = "<<", ">>"
PkgTemplates.view(::GitPagesDoc, t::Template, pkg::AbstractString) = Dict(
    "PKG" => pkg,
    "LOCAL_DOC" => deploy_samerepo(t, pkg) ? "yes" : nothing
)

PkgTemplates.make_canonical(::Type{GitPagesDoc}) = pages_url_stable

#
#  This was giving `WARNING: Method definition badges(PkgTemplates.GitHubActions)`
#  during precompilation. So we simply handle the Badge ouself in the README plugin.
#
# Type piracy. All URL generated are of form github.com/{{{USER}}}/{{{PKG}}}.jl
# badges(::GitHubActions : The .jl suffix is always present, event if the Git plugin has jl=false
# badges(::Codecov) : the branch name 'master' is hardcoded! We correct both here.
# badges(::Documenter{GitPagesDoc} : need custom doc URL
# function PkgTemplates.badges(::GitHubActions)
#     Badge("Build Status",
#     "https://github.com/{{{USER}}}/{{{PKG_REPO}}}/workflows/Tests/badge.svg",
#     "https://github.com/{{{USER}}}/{{{PKG_REPO}}}/actions")
# end
# function PkgTemplates.badges(::Codecov)
#     Badge("Coverage",
#     "https://codecov.io/gh/{{{USER}}}/{{{PKG_REPO}}}/branch/main/graph/badge.svg{{{BADGE_TOKEN}}}",
#     "https://codecov.io/gh/{{{USER}}}/{{{PKG_REPO}}}")
# end
# function PkgTemplates.badges(::Documenter{GitPagesDoc})
#     [Badge("Stable", "https://img.shields.io/badge/docs-stable-blue.svg", "{{{DOC_URL}}}/stable"),
#      Badge("Dev",    "https://img.shields.io/badge/docs-dev-blue.svg",    "{{{DOC_URL}}}/dev")   ]
# end

function PkgTemplates.user_view(::Documenter{GitPagesDoc}, t::Template, pkg::AbstractString)
    Dict(
        "HAS_DEPLOY" => true,
        # The repo into which the doc is deployed
        "DEPLOY_REPO" => userrepo_docs(t, pkg),
        # A sub-dir of root into which the doc is deployed (maybe nothing)
        "DEPLOY_DIR" => deploy_docdir(t, pkg),
        # The branch of the doc repo under which we deploy
        "DEPLOY_BRC" => deploy_branch(t, pkg),
        # The source code repo from which the doc is built
        # Used to build (checkout the repo)
        # Used in the doc to link back to the code repo
        "REPO" => hostuserrepo_code(t,pkg),

    )
end

function PkgTemplates.user_view(::Readme, t::Template, pkg::AbstractString)
    reponame = package_reponame(t, pkg)
    Dict(
        "TODAY" => today(),
        "USER"  => t.user,
        "PKG_REPO" => reponame,
        "DOC_URL" => pages_url_base(t, pkg),
        "BADGE_TOKEN" => isprivaterepo(t) ? "?token=SOME_10_CHARS_STRING" : "",
        "REMOTE_DOC_URL" => deploy_samerepo(t, pkg) ? nothing : "https://" * hostuserrepo_docs(t, pkg),
        "REMOTE_DOC" => deploy_samerepo(t, pkg) ? nothing : userrepo_docs(t, pkg),
        "PRIVATE_REPO" => isprivaterepo(t) ? "yes" : nothing,
        "CODECOV_HOME" => "https://app.codecov.io/gh/$(t.user)",
        "CODECOV_PKG_SETTING" => "https://app.codecov.io/gh/$(t.user)/$reponame/settings",
        "GITHUB_PKG_SETTING_SECRETS" => "https://github.com/$(t.user)/$reponame/settings/secrets/actions",
    )
end


text:
using CSV, Distributions, HypothesisTests

data = CSV.read("../data/machine1.csv", header=false)[:,1]
xBar, n = mean(data), length(data)
sigma = 1.2
mu0A, mu0B = 52.2, 53

testStatistic = ( xBar - mu0A ) / ( sigma/sqrt(n) )
pVal = 2*ccdf(Normal(), abs(testStatistic))

testA = OneSampleZTest(xBar, sigma, n, mu0A)
testB = OneSampleZTest(xBar, sigma, n, mu0B)

println("Results for mu0 = ", mu0A,":")
println("Manually calculated test statistic: ", testStatistic)
println("Manually calculated p-value: ", pVal,"\n")
println(testA)

println("\n In case of  mu0 = ", mu0B, " then p-value = ", pvalue(testB))

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="first", value = 1),
    dcc_input(id="second", value = 1),
    dcc_input(id="third", value = 1)
end

callback!(
    ClientsideFunction("clientside", "add1_no_update_at_11"),
    app, callid"{second.value, third.value} first.value=>second.value, third.value"
    )

run_server(app)


text:
# Author: Lukasz Bratos

# import Pkg
# Pkg.add("Polynomials")
using Polynomials

# wspolczynniki wielomianu
p = [1, -210.0, 20615.0,-1256850.0, 53327946.0,-1672280820.0, 40171771630.0,
    -756111184500.0, 11310276995381.0, -135585182899530.0, 1307535010540395.0,
    -10142299865511450.0, 63030812099294896.0, -311333643161390640.0,
    1206647803780373360.0, -3599979517947607200.0, 8037811822645051776.0,
    -12870931245150988800.0, 13803759753640704000.0, -8752948036761600000.0,
    2432902008176640000.0]

# wspolczynniki postaci iloczynowej
p_roots = [x for x in 1.0:20.0]
println(p_roots)

# wspolczynniki zmodyfikowanego wielomianu
p_mod = [1, -210.0 - 2 ^ -23, 20615.0,-1256850.0, 53327946.0,-1672280820.0, 40171771630.0,
-756111184500.0, 11310276995381.0, -135585182899530.0, 1307535010540395.0,
-10142299865511450.0, 63030812099294896.0, -311333643161390640.0,
1206647803780373360.0, -3599979517947607200.0, 8037811822645051776.0,
-12870931245150988800.0, 13803759753640704000.0, -8752948036761600000.0,
2432902008176640000.0]

# wielomianu w postaci normalnej
A = Poly(reverse(p))
# wielomian w postaci iloczynowej
B = poly(p_roots)
# zmodyfikowany wielomian w postaci normalnej
C = Poly(reverse(p_mod))

# pierwiastki wielomianow
A_roots = reverse(roots(A))
C_roots = reverse(roots(C))

for i in 1:20
    println("$i & $(A_roots[i]) & $(abs(polyval(A, A_roots[i]))) & $(abs(polyval(B, A_roots[i]))) & ", abs(A_roots[i] - i), " \\\\ \\hline")
end

for i in 1:20
    println("$i & $(C_roots[i]) & $(abs(polyval(C, C_roots[i]))) & ", abs(C_roots[i] - i), " \\\\ \\hline")
end


text:
using HTTP

HTTP.listen("0.0.0.0", 8000; max_connections=500) do http
    HTTP.setstatus(http, 200)
    write(http, "OK")
    return
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
module SymbolicCodegen

using SymbolicUtils
using SymbolicUtils: Sym, Term, Symbolic, operation, arguments, symtype
const RW = SymbolicUtils.Rewriters
const MaybeSym{T} = Union{T, Symbolic{T}}
using MacroTools
using MacroTools: @q
using RuntimeGeneratedFunctions

RuntimeGeneratedFunctions.init(@__MODULE__)

include("utils.jl")
include("base.jl")
include("sums.jl")
include("arrays.jl")
include("cse.jl")
include("foldconstants.jl")

export codegen

function codegen(s::Symbolic{T}) where {T}
    assignments = cse(s)

    q = @q begin end

    for a in assignments
        x = a[1]
        rhs = _codegen(a[2])
        push!(q.args, @q begin $x = $rhs end)
    end

    MacroTools.flatten(q)
end

ACRULES = [
    @acrule (~a + ~b)*(~c) => (~a) * (~c) + (~b) * (~c)
]

# export rewrite

# function rewrite(s)
#     simplify(s; polynorm=true) |> RW.Fixpoint(RW.Prewalk(RW.Chain(RULES))) |> simplify
# end

end


text:
using NonstationaryProcesses
using Plots

lorenz = lorenzSim(
    X0 = [0.0, -0.01, 9.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (10.0, 28.0, 8/3), # Sprott's recomendation
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.01,
    tmax = 8000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))

𝔩𝔬𝔯𝔢𝔫𝔷 = corruptphase(lorenz, rampInterval, (0.0, 0.2, 1000.0, 8000.0), nwindows = 160)()

a = animatespectrum(𝔩𝔬𝔯𝔢𝔫𝔷, downsample=500, trail=5000, colorgradient=cgrad([:black, :crimson]), phasogram=true, nwindows=160, dpi=100)
gif(a, "./thresholdcorruptlorenzsym.mp4", fps=48)


text:
export angular_distance, along_line_angular_distance, angular_length
# export angular_distance, distance, along_line_angular_distance, angular_length
# """
#     Radius Earth [m]
#
# Source: en.wikipedia.org/wiki/Earth_radius
# """
# const Rₑ_m = 6371008.7714

"""
    angular_distance(point₁::Point, point₂::Point)

Return the `angular_distance` [deg] of the great circle line between the
positions `point₁` [deg] and `point₂` [deg] on a unit sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function angular_distance(point₁::Point, point₂::Point)
    Δpos = point₂ - point₁
    a = sind(Δpos.ϕ/2.0)^2 + cosd(point₁.ϕ)*cosd(point₂.ϕ)*sind(Δpos.λ/2.0)^2
    c = 2.0 * atand(√a, √(1.0 - a))
    return c
end

"""
    angular_length(arc::Arc)

Return the `angular_length` [deg] of the great circle line section on a
unit sphere.
"""
angular_length(arc::Arc) =
angular_distance(arc.point₁, arc.point₂)

"""
    angular_length(points::Vector{Points{Float64}})

Return the `angular_length` [deg] along the subsequent points on a unit sphere.
"""
function angular_length(points::Vector{Point{Float64}})
    point₁ = points[1]
    dist_deg = 0.0
    for point₂ in points[2:end]
        dist_deg += angular_distance(point₁, point₂)
        point₁ = point₂
    end
    return dist_deg
end

"""
    angular_length(arcs::Arcs)

Return the `angular_length` [deg] along the arcs on a unit sphere.
"""
angular_length(arcs::Arcs) = angular_length(arcs.points)

"""
    angular_length(polygon::Polygon)

Return the `angular_length` [deg] (e.g. circumference) along the polygon on a unit
sphere.
"""
angular_length(polygon::Polygon) = angular_length(polygon.points)

# """
#     distance(angular_distance::Float64, radius::Float64=Rₑ_m)
#
# Convert the `angular_distance` [deg] into a distance on a sphere of given radius
#  [m]. The default radius is the radius of the Earth.
# """
# function distance(angular_distance::Float64, radius::Float64=Rₑ_m)
#     return deg2rad(angular_distance) * radius
# end

"""
    angular_distance(point₃::Point, point₁::Point, azimuth₁₂::Float64)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
great circle line starting in `point₁` [deg] with `azimuth₁₂` [deg] from
`point₁` to `point₂`. It is assumed that the great circle line does not stop in
`point₂`, but continuous around the unit sphere. A positive value indicates
being right of the line, and negative being left of the line.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function angular_distance(point₃::Point, point₁::Point, azimuth₁₂::Float64)
    angular_distance₁₃ = angular_distance(point₁, point₃)
    azimuth₁₃ = azimuth(point₁, point₃)
    return angular_distance(angular_distance₁₃, azimuth₁₃, azimuth₁₂)
end

"""
    angular_distance(point₃::Point, line₁::Line)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
great circle line starting in `point₁` [deg] with `azimuth₁₂` [deg]. A positive
value indicates being right of the line, and negative being left of the line.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
angular_distance(point₃::Point, line₁::Line) = angular_distance(point₃,
line₁.point, line₁.azimuth)

"""
    angular_distance(angular_distance₁₃::Float64, azimuth₁₃::Float64,
    azimuth₁₂::Float64)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
great circle line starting in `point₁` [deg] with `azimuth₁₂` [deg]. For the
calculation we need the `angular_distance₁₃` and `azimuth₁₃` between `point₁`
and `point₃` [deg], and the `azimuth₁₂` from `point₁` to `point₂` along the
great circle line. It is assumed that the great circle line does not stop in
`point₂`, but continuous around the unit sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function angular_distance(angular_distance₁₃::Float64, azimuth₁₃::Float64,
azimuth₁₂::Float64)
    return asind(sind(angular_distance₁₃) * sind(azimuth₁₃ - azimuth₁₂))
end

"""
    angular_distance(point₃::Point, point₁::Point, point₂::Point)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
great circle line section starting in `point₁` [deg] and ending in `point₂` [deg].
The great circle line section does not continue around the unit sphere.

The `angular_distance` does not change sign when being left or right of the arc.
"""
function angular_distance(point₃::Point, point₁::Point, point₂::Point, tolerance::Float64=tolerance_deg)
    angular_distance₁₂ = angular_distance(point₁, point₂)
    angular_distance₁₃ = angular_distance(point₁, point₃)
    angular_distance₂₃ = angular_distance(point₂, point₃)
    azimuth₁₂ = azimuth(point₁, point₂)
    angular_distance_line3 = angular_distance(point₃, point₁, azimuth₁₂)
    # along_line_angular_dist = along_line_angular_distance(angular_distance₁₃,
    # angular_distance_line3)
    along_line_pnt_3a = along_line_point(point₃, Line(point₁, azimuth₁₂))
    if abs(angular_distance(along_line_pnt_3a, point₁) +
        angular_distance(along_line_pnt_3a, point₂) -
        angular_distance₁₂) < tolerance
        return abs(angular_distance_line3)
    else
        return min(angular_distance₁₃, angular_distance₂₃)
    end
end

"""
    angular_distance(point₃::Point, points::Vector{Points{T}}) where T<:Float64

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
an set of points representing an set of arcs.

The `angular_distance` does not change sign when being left or right of the arc.
"""
function angular_distance(point₃::Point, points::Vector{Point{T}}) where
    T<:Float64
    point₁ = points[1]
    point₂ = points[2]
    min_dist = angular_distance(point₃, point₁, point₂)
    point₁ = point₂
    for point₂ in points[2:end]
        min_dist = min(min_dist, angular_distance(point₃, point₁, point₂))
        point₁ = point₂
    end
    return min_dist
end

"""
    angular_distance(point₃::Point, arcs::Arcs)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
an `arcs`.

The `angular_distance` does not change sign when being left or right of the arc.
"""
angular_distance(point₃::Point, arcs::Arcs) = angular_distance(point₃, arcs.points)

"""
    angular_distance(point₃::Point, polygon::Polygon)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on
the border of a `polygon`.

The `angular_distance` does not change sign when being left or right of the arc.
"""
angular_distance(point₃::Point, polygon::Polygon) = angular_distance(point₃,
polygon.points)

"""
    angular_distance(point₃::Point, arc₁₂::Arc)

Return the `angular_distance` [deg] from point₃ [deg] to the closest point on a
great circle line section starting in `point₁` [deg] and ending in `point₂` [deg].
The great circle line section does not continue around the unit sphere.
"""
angular_distance(point₃::Point, arc₁₂::Arc) = angular_distance(point₃, arc₁₂.point₁, arc₁₂.point₂)

"""
    along_line_angular_distance(angular_distance₁₃::Float64,
    angular_distance_line3::Float64)

The `along_line_angular_distance` from the start `point₁` [deg] to the closest
point on the great circle line to `point₃`. As input we need the
`angular_distance₁₃` [deg] between `point₁` and `point₃` [deg] and the closest
distance `angular_distance_line3` [deg] between the great circle line and
`point3` [deg]. The distance is non-directional.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function along_line_angular_distance(angular_distance₁₃::Float64,
angular_distance_line3::Float64)
    return acosd(cosd(angular_distance₁₃) / cosd(angular_distance_line3))
end

"""
    along_line_angular_distance(point₃::Point, line₁::Line)

The `along_line_angular_distance` from the start `point₁` [deg] to the closest
point on the great circle line `line₁` to `point₃`. The distance is
non-directional.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function along_line_angular_distance(point₃::Point, line₁::Line)
    angular_distance₁₃ = angular_distance(line₁.point, point₃)
    angular_distance_line3 = angular_distance(point₃::Point, line₁::Line)
    return along_line_angular_distance(angular_distance₁₃,
    angular_distance_line3)
end


text:
# If we list all the natural numbers below 10 that are multiples of 3 or 5, we
# get 3, 5, 6 and 9. The sum of these multiples is 23.
#
# Find the sum of all the multiples of 3 or 5 below 1000.

using ProjectEulerSolutions

# Brute force approach, loop through checking for multiples.
function p001solution_loop(n::Integer=9)::Integer
    total_sum = 0
    @simd for i in 1:n
        if i % 3 == 0 || i % 5 == 0
            total_sum += i
        end
    end
    return total_sum
end

# Faster approach, sum of real numbers in a sequence is n*(n+1)/2, just scale
# and remove double counts.
function p001solution_closedform(n::Integer=9)::Integer
    n3  = div(n, 3)
    n5  = div(n, 5)
    n15 = div(n, 15)
    return  3 * div(n3 * (n3 + 1), 2) +
            5 * div(n5 * (n5 + 1), 2) -
           15 * div(n15 * (n15 + 1), 2)
end

p001 = Problems.Problem(Dict("Closed form" => p001solution_closedform,
                             "Loop" => p001solution_loop))

Problems.benchmark(p001, 999)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using Documenter
using Plots
using TimeSeriesClustering

makedocs(sitename="TimeSeriesClustering.jl",
    authors = "Holger Teichgraeber, and Elias Kuepper",
    pages = [
        "Introduction" => "index.md",
        "Quick Start Guide" => "quickstart.md",
        "Load Data" => "load_data.md",
        "Representative Periods" => "repr_per.md",
        "Optimization" => "opt.md"
        ],
    format = Documenter.HTML(assets=["assets/clust_for_opt_text.svg"])
    )

deploydocs(repo = "github.com/holgerteichgraeber/TimeSeriesClustering.jl.git", devbranch="dev")


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_generators_variability(setup::Dict, path::AbstractString, inputs_genvar::Dict)

Function for reading input parameters related to hourly maximum capacity factors for all generators (plus storage and flexible demand resources)
"""
function load_generators_variability(setup::Dict, path::AbstractString, inputs_genvar::Dict)
	# Hourly capacity factors
	#data_directory = chop(replace(path, pwd() => ""), head = 1, tail = 0)
	data_directory = joinpath(path, setup["TimeDomainReductionFolder"])
	if setup["TimeDomainReduction"] == 1  && isfile(joinpath(data_directory,"Load_data.csv")) && isfile(joinpath(data_directory,"Generators_variability.csv")) && isfile(joinpath(data_directory,"Fuels_data.csv")) # Use Time Domain Reduced data for GenX
		gen_var = DataFrame(CSV.File(joinpath(data_directory,"Generators_variability.csv"), header=true), copycols=true)
	else # Run without Time Domain Reduction OR Getting original input data for Time Domain Reduction
		gen_var = DataFrame(CSV.File(joinpath(path,"Generators_variability.csv"), header=true), copycols=true)
	end

	# Reorder DataFrame to R_ID order (order provided in Generators_data.csv)
	select!(gen_var, [:Time_Index; Symbol.(inputs_genvar["RESOURCES"]) ])

	# Maximum power output and variability of each energy resource
	inputs_genvar["pP_Max"] = transpose(Matrix{Float64}(gen_var[1:inputs_genvar["T"],2:(inputs_genvar["G"]+1)]))

	println("Generators_variability.csv Successfully Read!")

	return inputs_genvar
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
import Dance

include("./utils/main.jl")

Dance.start_project("demo")
project_settings()
routes("html")

## Test all routes with pending slash, to ensure is removed ##
@testset "HTTP.listen" begin
    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(1)

    make_and_test_request_get("/", 200, Dict("Content-Type" => "text/html; charset=UTF-8"), 217, false, "Hello World")

    # Test decimal url param
    make_and_test_request_get("/dict/12.3/", 200, Dict("Content-Type" => "text/html; charset=UTF-8"), 216, true, Dict("a" => 12.3))

    # Test int and string url params
    make_and_test_request_get("/dict/abc/123/", 200, Dict("Content-Type" => "text/html; charset=UTF-8"), 217, true, Dict("abc" => 123))

    make_and_test_request_get(
        "/dataframe/",
        200,
        Dict("Content-Type" => "text/html; charset=UTF-8"),
        249,
        true,
        [
            ["A", "B"],
            [1, "M"],
            [2, "F"],
            [3, "F"],
            [4, "M"]
        ]
    )
end

delete_project()


text:
using Test

include("Grids.jl")
include("StateVecs.jl")
include("StateVecFuncs.jl")

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
using DataFrames
using Test

m = page_model()
include("../src/components/MarketDamages.jl")

marketdamages = addmarketdamages(m)

set_param!(m, :MarketDamages, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/rtl_realizedtemperature.csv"))
set_param!(m, :MarketDamages, :rcons_per_cap_SLRRemainConsumption, readpagedata(m,"test/validationdata/rcons_per_cap_SLRRemainConsumption.csv"))
set_param!(m, :MarketDamages, :rgdp_per_cap_SLRRemainGDP, readpagedata(m,"test/validationdata/rgdp_per_cap_SLRRemainGDP.csv"))
set_param!(m, :MarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m,"test/validationdata/atl_adjustedtolerableleveloftemprise_market.csv"))
set_param!(m, :MarketDamages, :imp_actualreduction, readpagedata(m,"test/validationdata/imp_actualreduction_market.csv"))
set_param!(m, :MarketDamages, :isatg_impactfxnsaturation, 28.333333333333336)

p = load_parameters(m)
p["y_year_0"] = 2008.
p["y_year"] = Mimi.dim_keys(m.md, :time)
set_leftover_params!(m, p)

run(m)

rcons_per_cap = m[:MarketDamages, :rcons_per_cap_MarketRemainConsumption]
rcons_per_cap_compare = readpagedata(m, "test/validationdata/rcons_per_cap_MarketRemainConsumption.csv")
@test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-1

rgdp_per_cap = m[:MarketDamages, :rgdp_per_cap_MarketRemainGDP]
rgdp_per_cap_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_MarketRemainGDP.csv")
@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2


text:
# test ScatterAlignPose2

using Test
using IncrementalInference
using Distributions

##
@testset "Test HeatmapGridDensity" begin
##

x = -10:0.1:10;
y = -10:0.1:10;

img = zeros(length(x),length(y))

g = (x,y)->pdf(MvNormal([1.0;1.0]),[x;y])

for (i,x_) in enumerate(x), (j,y_) in enumerate(y)
  img[i,j] = g(x_,y_)
end

hgd = IIF.HeatmapGridDensity(img,(x,y),nothing, 0.07, N=1000)

# check conversions to packed types
phgd = convert(PackedSamplableBelief, hgd)
hgd_ = convert(SamplableBelief, phgd)

@test isapprox( hgd, hgd_ )

# use in graph

##
end


text:
#=
    Top level program containing the function main which runs the c_sw kernel.
=#

module program

include("./SwCore.jl")
include("./states.jl")

using TOML, Printf
using NCDatasets, OffsetArrays
using .SwCoreModule

export main

    function main(input::Dict)

    # Input configuration variables
    datasize = input["name"]
    nc_input_file = input["input_file"]
    print_output_file = input["test_out_file"]
    nc_output_file = input["output_file"]
    interpFactor = input["interpFactor"]

    # Open the statistics log file
    io = open(print_output_file,"w")

    # Define the dataset using NCDatasets pkg for each datafile
    ds = NCDataset(nc_input_file, "a")

    # Assign Variables from the NetCDF file to a ::State Julia Struct
    current_state = State(ds)

    # Print the input state statistics to the log file
    print_state("Input State - Original", current_state, io)

    # Interpolate the data according to the interpolation factor.
    if (interpFactor > 0)
      interpolate_state!(current_state, interpFactor)
      # Write the input state statistics to the log
      print_state("Input State - Interpolated", current_state, io)
    elseif (interpFactor == 0)
      # Do nothing.
    else
      println("Error, InterpFactor less than zero.")
    end

    # Run the kernel
    println("Run the kernel :  $datasize  , Num Threads : ", Threads.nthreads())
    @time Threads.@threads for k = 1 : current_state.npz
        c_sw!(current_state, k)
    end

    # Print the output state statistics to the log file
    print_state("Output State", current_state, io)

    # Write a new NetCDF file
    println("Write NetCDF file $datasize : ")
    @time write_state(nc_output_file, current_state)

    # Close the statistics log file
    close(io)

    return true

    end # function main

end # module program


text:
function transform(data_dictionary::Dict{String,Any}, system_flag::String)::DataFrame

    # initialize -
    df = DataFrame(biophysical_parameter_id=String[], parameter_key=String[],
        parameter_value=String[], parameter_units=String[], parameter_source_reference = Union{String,Missing}[],
        system_key = String[])

    try

        # ok, so lets populate the data from the dictionary -
        for (key, dict_value) in data_dictionary

            # build a local array, we'll add stuff to the df row wise -
            local_data_row = Array{Union{String,Missing},1}()

            # build a UUID for this record -
            record_id = string(uuid4())
            push!(local_data_row,record_id)
            push!(local_data_row,key)

            # get the value -
            parameter_value = dict_value["value"]
            push!(local_data_row, parameter_value)

            # get the units -
            parameter_units = dict_value["units"]
            push!(local_data_row, parameter_units)

            # parameter reference -
            parameter_reference = dict_value["bionumbers_reference"]
            push!(local_data_row, parameter_reference)

            # what is the system key -
            push!(local_data_row, system_flag)

            # add to the df -
            push!(df, tuple(local_data_row...))
        end

        # return -
        return df
    catch error
        throw(error)
    end
end


text:
# Utility functions
function trapz(x,y)
   n = length(y)
   if n != length(x)
       error("Input x,y must be of same length")
   end
       r = zero(zero(eltype(x))*zero(eltype(y)))
   for i in 2:n
       r += (x[i] - x[i-1]) * (y[i] + y[i-1])
   end
       return r*0.5
end


text:
a = 1
b = 5
println(a + b)


text:
module RealisticLESbrary

export InterpolatedProfile, InterpolatedProfileTimeSeries, ∂z, ∂t
export validate_sose_dates, dates_to_offset

using Dates

include("interpolated_profiles.jl")

const sose_start_date = Date(2013, 1, 1)
const sose_end_date = Date(2018, 1, 1)

function validate_sose_dates(start_date, end_date)
    @assert start_date < end_date
    @assert start_date >= sose_start_date
    @assert end_date <= sose_end_date
end

function dates_to_offset(start_date, end_date; step=1)
    offset = (start_date - sose_start_date).value / step + 1
    n_dates = (end_date - start_date).value / step
    return floor(Int, offset), ceil(Int, n_dates)
end

end # module


text:
#!/usr/bin/env julia
push!(LOAD_PATH, ".")
using HelloWorld
using Test

@test HelloWorld.greet("John") == "Hello, John"

max_stars_repo_name: vtavernier/jomw

max_issues_repo_name: vtavernier/jomw

max_forks_repo_name: vtavernier/jomw


text:
export hist_eig

#Uses the method of Sturm sequences to compute the histogram of eigenvalues of a matrix
#
# Reference
#   Cy P. Chan, "Sturm Sequences and the Eigenvalue Distribution of
#   the Beta-Hermite Random Matrix Ensemble", M.Sc. thesis (MIT), 2007
#
# For better numerical stability, this algorithm is modified so that only the sign
# of the determinant is stored.
#
#For general matrices, this is slower than hist(eigvals(M), bins) and is NOT recommended
function hist_eig{GridPoint <: Number}(M::AbstractMatrix, bins::Vector{GridPoint})
  n = size(M)[1]
  NumBins = length(bins)
  histogram = zeros(NumBins)
  SturmSequence = zeros(n)
  for BinId=1:NumBins
    K = M - bins[BinId]*eye(n)
    #Progression of determinants of lower-right submatrices
    #SturmSequence = [det(K[end-j+1:end,end-j+1:end]) for j=1:n]
    #SturmRatioSequence = [SturmSequence[j+1]/SturmSequence[j] for j=1:n]

    SturmSequenceSigns = [int(sign(det(K[end-j+1:end,end-j+1:end]))) for j=1:n]
    SturmSequenceSigns = [1; SturmSequenceSigns]
    SturmRatioSequence = [SturmSequenceSigns[j+1]/SturmSequenceSigns[j] for j=1:n]
    histogram[BinId] = sum([r < 0 for r in SturmRatioSequence])
  end
  histogram = int(diff([histogram; n]))
end

#Uses the method of Sturm sequences to compute the histogram of eigenvalues of
#a symmetric tridiagonal mmatrix
#
# Reference
#   Cy P. Chan, "Sturm Sequences and the Eigenvalue Distribution of
#   the Beta-Hermite Random Matrix Ensemble", M.Sc. thesis (MIT), 2007
#
function hist_eig{GridPoint <: Number}(M::SymTridiagonal, bins::Vector{GridPoint})
  n = size(M)[1]
  NumBins = length(bins)
  histogram = zeros(NumBins)
  r = zeros(n)
  for BinId=1:NumBins
    a, b = M.dv - bins[BinId], M.ev
    #Formula (2.3)
    r[1] = a[1]
    for i=2:n
      r[i] = a[i] - b[i-1]^2/r[i-1]
    end
    histogram[BinId] = sum([SturmRatio < 0 for SturmRatio in r])
  end
  histogram = int(diff([histogram; n]))
end


text:
using Test
using Quiqbox

@testset "OneBody.jl tests" begin
    errT1 = 1e-10
    errT2 = 1e-10
    nucs = fill("H", 2)
    cens = [[-0.7, 0.0, 0.0], [ 0.7, 0.0, 0.0]]
    bf1 = genBasisFunc(cens[1], ("STO-3G", "H"))[]
    bf2 = genBasisFunc(cens[2], ("STO-3G", "H"))[]
    bs = [bf1, bf2]

    S = [1.0 0.6593182058508896; 0.6593182058508896 1.0]
    @test isapprox(overlap(bf1, bf1), ones(1,1), atol=errT1)
    @test isapprox(overlap(bf1, bf2), S[2]*ones(1,1), atol=errT1)
    @test isapprox(overlaps(bs), S)

    V = [-1.8804408905227634 -1.1948346220535715; -1.1948346220535715 -1.8804408905227632]
    @test isapprox(nucAttraction(bf1, bf2, nucs, cens), V[2]*ones(1,1), atol=errT1)
    @test isapprox(nucAttractions(bs, nucs, cens), V, atol=errT1)

    T = [0.7600318799755844 0.23645465829079276; 0.23645465829079276 0.7600318799755844]
    @test isapprox(elecKinetic(bf1, bf2), T[2]*ones(1,1,1), atol=errT1)
    @test isapprox(elecKinetics(bs), T)

    Hc = coreH(bs, nucs, cens)
    @test isapprox(Hc, T+V, atol=errT1)
    @test isapprox(coreHij(bf1, bf1, nucs, cens)[], Hc[1], atol=errT2)
    @test isapprox(coreHij(bf2, bf1, nucs, cens)[],
                   coreHij(bf1, bf2, nucs, cens)[], atol=errT2)
    @test isapprox(coreHij(bf2, bf1, nucs, cens)[], Hc[3], atol=errT2)
    @test isapprox(coreHij(bf2, bf2, nucs, cens)[], Hc[4], atol=errT2)
end


text:
struct Constraint{T}
    evaluate::Any
    jacobian::Any
    hessian::Any
    num_state::Int
    num_action::Int
    num_parameter::Int
    num_constraint::Int
    num_jacobian::Int
    num_hessian::Int
    jacobian_sparsity::Vector{Vector{Int}}
    hessian_sparsity::Vector{Vector{Int}}
    evaluate_cache::Vector{T}
    jacobian_cache::Vector{T}
    hessian_cache::Vector{T}
    indices_inequality::Vector{Int}
end

Constraints{T} = Vector{Constraint{T}} where T

function Constraint(f::Function, num_state::Int, num_action::Int;
    num_parameter::Int=0,
    indices_inequality=collect(1:0),
    evaluate_hessian=false)

    #TODO: option to load/save methods
    @variables x[1:num_state], u[1:num_action], w[1:num_parameter]
    evaluate = f(x, u, w)
    jac = Symbolics.sparsejacobian(evaluate, [x; u])
    evaluate_func = eval(Symbolics.build_function(evaluate, x, u, w)[2])
    jac_func = eval(Symbolics.build_function(jac.nzval, x, u, w)[2])
    num_constraint = length(evaluate)
    num_jacobian = length(jac.nzval)
    jacobian_sparsity = [findnz(jac)[1:2]...]
    if evaluate_hessian
        @variables λ[1:num_constraint]
        lag_con = dot(λ, evaluate)
        hessian = Symbolics.sparsehessian(lag_con, [x; u])
        hess_func = eval(Symbolics.build_function(hessian.nzval, x, u, w, λ)[2])
        hessian_sparsity = [findnz(hessian)[1:2]...]
        num_hessian = length(hessian.nzval)
    else
        hess_func = Expr(:null)
        hessian_sparsity = [Int[]]
        num_hessian = 0
    end

    return Constraint(
        evaluate_func,
        jac_func,
        hess_func,
        num_state,
        num_action,
        num_parameter,
        num_constraint,
        num_jacobian,
        num_hessian,
        jacobian_sparsity,
        hessian_sparsity,
        zeros(num_constraint),
        zeros(num_jacobian),
        zeros(num_hessian),
        indices_inequality)
end

function Constraint()
    return Constraint(
        (constraint, state, action, parameter) -> nothing,
        (jacobian,   state, action, parameter) -> nothing,
        (hessian,    state, action, parameter) -> nothing,
        0, 0, 0, 0, 0, 0,
        [Int[], Int[]],
        [Int[], Int[]],
        Float64[],
        Float64[],
        Float64[],
        collect(1:0))
end

function constraints!(violations, indices, constraints::Constraints{T}, states, actions, parameters) where T
    for (t, con) in enumerate(constraints)
        con.evaluate(con.evaluate_cache, states[t], actions[t], parameters[t])
        @views violations[indices[t]] .= con.evaluate_cache
        fill!(con.evaluate_cache, 0.0) # TODO: confirm this is necessary
    end
end

function jacobian!(jacobians, indices, constraints::Constraints{T}, states, actions, parameters) where T
    for (t, con) in enumerate(constraints)
        con.jacobian(con.jacobian_cache, states[t], actions[t], parameters[t])
        @views jacobians[indices[t]] .= con.jacobian_cache
        fill!(con.jacobian_cache, 0.0) # TODO: confirm this is necessary
    end
end

function hessian_lagrangian!(hessians, indices, constraints::Constraints{T}, states, actions, parameters, duals) where T
    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_cache)
            con.hessian(con.hessian_cache, states[t], actions[t], parameters[t], duals[t])
            @views hessians[indices[t]] .+= con.hessian_cache
            fill!(con.hessian_cache, 0.0) # TODO: confirm this is necessary
        end
    end
end

function sparsity_jacobian(constraints::Constraints{T}, num_state::Vector{Int}, num_action::Vector{Int};
    row_shift=0) where T

    row = Int[]
    col = Int[]

    for (t, con) in enumerate(constraints)
        col_shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
        push!(row, (con.jacobian_sparsity[1] .+ row_shift)...)
        push!(col, (con.jacobian_sparsity[2] .+ col_shift)...)
        row_shift += con.num_constraint
    end

    return collect(zip(row, col))
end

function sparsity_hessian(constraints::Constraints{T}, num_state::Vector{Int}, num_action::Vector{Int}) where T
    row = Int[]
    col = Int[]

    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_sparsity[1])
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
            push!(row, (con.hessian_sparsity[1] .+ shift)...)
            push!(col, (con.hessian_sparsity[2] .+ shift)...)
        end
    end

    return collect(zip(row, col))
end

num_constraint(constraints::Constraints{T}) where T = sum([con.num_constraint for con in constraints])
num_jacobian(constraints::Constraints{T}) where T = sum([con.num_jacobian for con in constraints])
# num_hessian(constraints::Constraints{T}) where T = sum([con.num_hessian for con in constraints])

function constraint_indices(constraints::Constraints{T};
    shift=0) where T

    indices = Vector{Int}[]

    for (t, con) in enumerate(constraints)
        indices = [indices..., collect(shift .+ (1:con.num_constraint)),]
        shift += con.num_constraint
    end

    return indices
end

function jacobian_indices(constraints::Constraints{T};
    shift=0) where T

    indices = Vector{Int}[]

    for (t, con) in enumerate(constraints)
        push!(indices, collect(shift .+ (1:con.num_jacobian)))
        shift += con.num_jacobian
    end

    return indices
end

function hessian_indices(constraints::Constraints{T}, key::Vector{Tuple{Int,Int}}, num_state::Vector{Int}, num_action::Vector{Int}) where T
    indices = Vector{Int}[]
    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_sparsity[1])
            row = Int[]
            col = Int[]
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
            push!(row, (con.hessian_sparsity[1] .+ shift)...)
            push!(col, (con.hessian_sparsity[2] .+ shift)...)
            rc = collect(zip(row, col))
            push!(indices, [findfirst(x -> x == i, key) for i in rc])
        else
            push!(indices, Int[])
        end
    end
    return indices
end


text:
"""
Original publication:
L. Yin, H. Wang, and W. Fan. Active learning based
support vector data description method for robust
novelty detection. Knowledge-Based Systems, pages
40–52, Aug. 2018.
"""
struct BoundaryNeighborCombinationPQs{M <: SVDD.OCClassifier} <: HybridPQs
    occ::M
    nn_dist
    η::Float64
    p::Float64
end

BoundaryNeighborCombinationPQs{M}(args...;kwargs...) where M = BoundaryNeighborCombinationPQs(args...; kwargs...)

function BoundaryNeighborCombinationPQs(occ::M, x::Array{T, 2}; η=0.7, p=0.15) where {T <: Real, M <:SVDD.OCClassifier}
    BoundaryNeighborCombinationPQs(occ, knn_mean_dist(x; k=1), η, p)
end

function BoundaryNeighborCombinationPQs(occ::M, x::Array{T, 2}; η=0.7, p=0.15) where {T <: Real, M <: SVDD.SubOCClassifier}
    nn_dist = [knn_mean_dist(x[s,:]; k=1) for s in occ.subspaces]
    BoundaryNeighborCombinationPQs(occ, nn_dist, η, p)
end

function bnc(nn_dist, η, p, labels, predictions)
    pred = abs.(predictions)
    d_hyp = (pred .- minimum(pred[labels[:U]])) ./ maximum(pred[labels[:U]])
    d_nn = (nn_dist .- minimum(nn_dist[labels[:U]])) ./ maximum(nn_dist[labels[:U]])
    return η * (-d_hyp) + (1 - η) * (-d_nn)
end

function qs_score(qs::BoundaryNeighborCombinationPQs{M},
                  x::Array{T, 2},
                  labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where {T <: Real, M <: SVDD.OCClassifier}
    @assert length(qs.nn_dist) == size(x, 2)
    if rand() < qs.p
        return qs_score(RandomPQs(), x, labels)
    end
    return bnc(qs.nn_dist, qs.η, qs.p, labels, SVDD.predict(qs.occ, x))
end

function qs_score(qs::BoundaryNeighborCombinationPQs{M},
                  x::Array{T, 2},
                  labels::Dict{Symbol, Array{Int, 1}},
                  subspaces::Vector{Vector{Int}}) where {T <: Real, M <: SVDD.SubOCClassifier}
    length(qs.nn_dist) == length(subspaces) || throw(DimensionMismatch("Number of nn_dist arrays must match the number of subspaces."))
    all(length.(qs.nn_dist) .== size(x, 2)) || throw(DimensionMismatch("Number of neighbors in nn_dist must match number of observations in x."))
    size(qs.occ.data, 1) == size(x, 1) || throw(DimensionMismatch("Number of dimensions in x must match the number of dimensions of the classifier."))
    # return random for all observations (and not per subspace)
    if rand() < qs.p
        return qs_score(RandomPQs(), x, labels, subspaces)
    end
    map(i -> bnc(qs.nn_dist[i], qs.η, qs.p, labels, SVDD.predict(qs.occ, x[subspaces[i], :], i)), eachindex(subspaces))
end


text:
# Case I example: Bad diamond II
# Claudia January 2015

# types in "types.jl"
if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
end

tree = "((((8,10))#H1,7),6,(4,#H1));" # Case I Bad diamond II
#f = open("prueba_tree.txt","w")
#write(f,tree)
#close(f)
net = readTopologyLevel1(tree)

#printEdges(net)
#printNodes(net)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1ObjectMeta <: SwaggerModel
    annotations::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: annotations
    clusterName::Any # spec type: Union{ Nothing, String } # spec name: clusterName
    creationTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: creationTimestamp
    deletionGracePeriodSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: deletionGracePeriodSeconds
    deletionTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: deletionTimestamp
    finalizers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: finalizers
    generateName::Any # spec type: Union{ Nothing, String } # spec name: generateName
    generation::Any # spec type: Union{ Nothing, Int64 } # spec name: generation
    initializers::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Initializers } # spec name: initializers
    labels::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: labels
    managedFields::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry} } # spec name: managedFields
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace
    ownerReferences::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1OwnerReference} } # spec name: ownerReferences
    resourceVersion::Any # spec type: Union{ Nothing, String } # spec name: resourceVersion
    selfLink::Any # spec type: Union{ Nothing, String } # spec name: selfLink
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid

    function IoK8sApimachineryPkgApisMetaV1ObjectMeta(;annotations=nothing, clusterName=nothing, creationTimestamp=nothing, deletionGracePeriodSeconds=nothing, deletionTimestamp=nothing, finalizers=nothing, generateName=nothing, generation=nothing, initializers=nothing, labels=nothing, managedFields=nothing, name=nothing, namespace=nothing, ownerReferences=nothing, resourceVersion=nothing, selfLink=nothing, uid=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("annotations"), annotations)
        setfield!(o, Symbol("annotations"), annotations)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("clusterName"), clusterName)
        setfield!(o, Symbol("clusterName"), clusterName)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("creationTimestamp"), creationTimestamp)
        setfield!(o, Symbol("creationTimestamp"), creationTimestamp)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("deletionGracePeriodSeconds"), deletionGracePeriodSeconds)
        setfield!(o, Symbol("deletionGracePeriodSeconds"), deletionGracePeriodSeconds)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("deletionTimestamp"), deletionTimestamp)
        setfield!(o, Symbol("deletionTimestamp"), deletionTimestamp)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("finalizers"), finalizers)
        setfield!(o, Symbol("finalizers"), finalizers)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("generateName"), generateName)
        setfield!(o, Symbol("generateName"), generateName)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("generation"), generation)
        setfield!(o, Symbol("generation"), generation)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("initializers"), initializers)
        setfield!(o, Symbol("initializers"), initializers)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("labels"), labels)
        setfield!(o, Symbol("labels"), labels)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("managedFields"), managedFields)
        setfield!(o, Symbol("managedFields"), managedFields)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("ownerReferences"), ownerReferences)
        setfield!(o, Symbol("ownerReferences"), ownerReferences)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("resourceVersion"), resourceVersion)
        setfield!(o, Symbol("resourceVersion"), resourceVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("selfLink"), selfLink)
        setfield!(o, Symbol("selfLink"), selfLink)
        validate_property(IoK8sApimachineryPkgApisMetaV1ObjectMeta, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1ObjectMeta

const _property_map_IoK8sApimachineryPkgApisMetaV1ObjectMeta = Dict{Symbol,Symbol}(Symbol("annotations")=>Symbol("annotations"), Symbol("clusterName")=>Symbol("clusterName"), Symbol("creationTimestamp")=>Symbol("creationTimestamp"), Symbol("deletionGracePeriodSeconds")=>Symbol("deletionGracePeriodSeconds"), Symbol("deletionTimestamp")=>Symbol("deletionTimestamp"), Symbol("finalizers")=>Symbol("finalizers"), Symbol("generateName")=>Symbol("generateName"), Symbol("generation")=>Symbol("generation"), Symbol("initializers")=>Symbol("initializers"), Symbol("labels")=>Symbol("labels"), Symbol("managedFields")=>Symbol("managedFields"), Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"), Symbol("ownerReferences")=>Symbol("ownerReferences"), Symbol("resourceVersion")=>Symbol("resourceVersion"), Symbol("selfLink")=>Symbol("selfLink"), Symbol("uid")=>Symbol("uid"))
const _property_types_IoK8sApimachineryPkgApisMetaV1ObjectMeta = Dict{Symbol,String}(Symbol("annotations")=>"Dict{String, String}", Symbol("clusterName")=>"String", Symbol("creationTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("deletionGracePeriodSeconds")=>"Int64", Symbol("deletionTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("finalizers")=>"Vector{String}", Symbol("generateName")=>"String", Symbol("generation")=>"Int64", Symbol("initializers")=>"IoK8sApimachineryPkgApisMetaV1Initializers", Symbol("labels")=>"Dict{String, String}", Symbol("managedFields")=>"Vector{IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry}", Symbol("name")=>"String", Symbol("namespace")=>"String", Symbol("ownerReferences")=>"Vector{IoK8sApimachineryPkgApisMetaV1OwnerReference}", Symbol("resourceVersion")=>"String", Symbol("selfLink")=>"String", Symbol("uid")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1ObjectMeta }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1ObjectMeta))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1ObjectMeta }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1ObjectMeta[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1ObjectMeta }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1ObjectMeta[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1ObjectMeta)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1ObjectMeta }, name::Symbol, val)
end


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for PMCMC constructor.

# Fields
$(TYPEDFIELDS)
"""
struct ParticleMetropolis{P<:ParticleFilter,M<:MCMC} <: PMCMCKernel
    "Particle Filter kernel to estimate latent trajectory."
    pf::P
    "MCMC kernel to sample continuous model parameter."
    mcmc::M
    function ParticleMetropolis(pf::P, mcmc::M) where {P<:ParticleFilter,M<:MCMC}
        @argcheck isa(pf.tune.referencing, Marginal) "Cannot use Conditional Filter in Particle Metropolis setting, use Marginal referencing"
        #!NOTE: Only create warning, but let user decide if they want proceed.
        #@argcheck isa(mcmc.kernel, Metropolis) "Only gradient free mcmc kernels allowed as tuning is based on log-target approximation from particle filter instead of objective(val). For gradient based mcmc kernels, use ParticleGibbs instead.")
        if isa(mcmc.tune.stepsize.adaption, BaytesCore.UpdateTrue)
            println("If no objective(val) specified, it is recommended to use a fixed stepsize for the mcmc kernel, as otherwise the kernel will be tuned based on acceptance rate of mcmc kernel.")
        end
        return new{P,M}(pf, mcmc)
    end
end
function ParticleMetropolis(
    filter::F, mcmc::M; kwargs...
) where {F<:ParticleFilterConstructor,M<:MCMCConstructor}
    return PMCMC(ParticleMetropolis, filter, mcmc; kwargs...)
end

############################################################################################
"""
$(SIGNATURES)
Propose new parameter with pmcmc sampler. If update=true, objective function will be updated with input model and data.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    pmcmc::ParticleMetropolis,
    model::ModelWrapper,
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0),
    update::U=BaytesCore.UpdateTrue(),
) where {D,F<:AbstractFloat,U<:BaytesCore.UpdateBool}
    ## Compute initial logposterior and save initial model value
    ℓpostₜ =
        pmcmc.pf.particles.ℓobjective.cumulative +
        ModelWrappers.log_prior_with_transform(model, pmcmc.mcmc.tune.tagged)
    val = deepcopy(model.val)
    ## Update Objective with new model parameter from other MCMC samplers and/or new/latent data
    objective = Objective(model, data, pmcmc.mcmc.tune.tagged, temperature)
    ## Update Kernel with current objective/configs
    update!(pmcmc.mcmc.kernel, objective, update)
    ## Make MCMC Proposal step
    resultᵖ, divergent, accept, sampler_statistic = propagate(
        _rng, pmcmc.mcmc.kernel, pmcmc.mcmc.tune, objective
    )
    ## Make PF proposal step with proposed parameter
    if !divergent
        ModelWrappers.unflatten_constrain!(
            objective.model, pmcmc.mcmc.tune.tagged, resultᵖ.θᵤ
        )
    end
    _, pf_diagnostics = propose!(_rng, pmcmc.pf, objective.model, objective.data, temperature, update)
    ## Compute acceptance rate
    ℓpostₜᵖ =
        pmcmc.pf.particles.ℓobjective.cumulative +
        ModelWrappers.log_prior_with_transform(objective.model, pmcmc.mcmc.tune.tagged)
    acceptᵖ = BaytesCore.AcceptStatistic(_rng, model.info.reconstruct.default.output(ℓpostₜᵖ - ℓpostₜ))
    if acceptᵖ.accepted
        ## If accepted, update MCMC kernel
        pmcmc.mcmc.kernel.result = resultᵖ
    else
        ## Else set model back to initial value
        model.val = val
    end
    #!NOTE: Would need to acceptᵖ here, but makes sampler very noisy. Better to print to either use fixed stepsize or assign objective(val).
    update!(pmcmc.mcmc.tune, pmcmc.mcmc.kernel.result, accept.rate)
    mcmc_diagnostics = MCMCDiagnostics(
        BaytesCore.BaseDiagnostics(
            pmcmc.mcmc.kernel.result.ℓθᵤ,
            objective.temperature,
            predict(_rng, objective),
            pmcmc.mcmc.tune.iter.current
        ),
        sampler_statistic,
        divergent,
        acceptᵖ,
        generate(_rng, objective, pmcmc.mcmc.tune.generated),
    )
    ## Assign base diagnostics - ℓobjective and predictions are taken from particle filter
    diagnostics = BaytesCore.BaseDiagnostics(
        pf_diagnostics.base.ℓobjective, pf_diagnostics.base.temperature, pf_diagnostics.base.prediction,
        mcmc_diagnostics.base.iter
    )
    ## Return pmcmc output
    return model.val, PMCMCDiagnostics(diagnostics, pf_diagnostics, mcmc_diagnostics)
end

############################################################################################
#export
export ParticleMetropolis, propose!


text:
using MonteCarloX
using Random
using StatsBase
using Printf
using ProgressMeter
using LinearAlgebra
using DelimitedFiles
using StaticArrays
"""
# I just want to see if this works. OK it does :-)
"""
# generate a stochastic trajectory
#SIR dynamics is goverened by differential equation
#``
#    \\frac{dI}{dt} = \\sum_{i=1}^{I}\\lambda_i\\frac{S}{N} - \\mu I
#``

# Arguments
#* distribution for the heterogeneous rates P(lambda)
"""
function trajectory!(rng, list_T, list_S, list_I, list_R, system)
    rates = current_rates(system)
    pass_update!(rates, index) = update!(rates, index, system, rng)

    list_S[1] = system.measure_S
    list_I[1] = system.measure_I
    list_R[1] = system.measure_R
    time_simulation = Float64(list_T[1])
    for i in 2:length(list_T)
        if time_simulation < list_T[i]
            dT = list_T[i] - time_simulation
            dT_sim = advance!(KineticMonteCarlo(), rng, rates, pass_update!, dT)
            time_simulation += dT_sim
        end
        list_S[i] = system.measure_S
        list_I[i] = system.measure_I
        list_R[i] = system.measure_R
    end
end

"""
#Ingreedients:
# - the rate of each agent to become infected or to recover
# -> the rate of infected to make a connection with any other agent (well-mixed),
#    which adds to the connection matrix, and only infects with a certain
#    probability
#    !! Connection matrix needs to be deleted after 2 weeks -> LightGraphs plus a
#    queue where we store the edges and the time when stored. (clean the graph whenever it gets updated)
#
# - the rate to be tested for those that will develop symptoms (probability)
#     -> as soon as S->I in update we draw a random number, and if the id is symptom based then add this rate to the system
# - the rate to be tested randomly
# -> upon positive test: add contacts to a list of tracing
# - the rate to test people from tracing list
"""
mutable struct SIR{D}
    epsilon::Float64
    mu::Float64
    P_lambda::D
    current_lambda::Vector{Float64}
    sum_current_lambda::Float64
    update_current_lambda::Int
    N::Int
    S::Int
    I::Int
    R::Int
    measure_S::Int
    measure_I::Int
    measure_R::Int

    function SIR{D}(rng::AbstractRNG, P_lambda::D, epsilon::Number, mu::Number, S0::Int, I0::Int, R0::Int) where D
        N = S0 + I0 + R0
        current_lambda = rand(rng, P_lambda, I0)
        sum_current_lambda = sum(current_lambda)
        new(Float64(epsilon), Float64(mu), P_lambda, current_lambda, sum_current_lambda, 0, N, S0, I0, R0, S0, I0, R0)
    end
end

"""
# update rates and system according to the last event that happend:
# index -> event type and agent id
# e.g.: index = 1-10000 would mean update agent i out of 10.000 = N
#       index = 10001 -> would mean to random test
#       index = 10002++ -> would mean do a test from some list that we need to define
"""
function update!(rates::AbstractVector, index::Int, system::SIR, rng::AbstractRNG)
    system.measure_S = system.S
    system.measure_I = system.I
    system.measure_R = system.R
    if index == 1 # recovery
        random_I = rand(rng, 1:system.I)
        delete_from_system!(system, random_I)
        system.I -= 1
        system.R += 1
    elseif index == 2 # infection
        add_to_system!(system, rand(rng, system.P_lambda))
        system.S -= 1
        system.I += 1
        @assert length(system.current_lambda) == system.I
    elseif index == 0 # absorbing state of zero infected
        system.measure_S = system.S
        system.measure_I = system.I
        system.measure_R = system.R
    else
        throw(UndefVarError(:index))
    end

    rates .= current_rates(system)
end

function delete_from_system!(system, index)
    update_sum!(system, -system.current_lambda[index])
    deleteat!(system.current_lambda, index)
end

function add_to_system!(system, lambda)
    push!(system.current_lambda, lambda)
    update_sum!(system, lambda)
end

function update_sum!(system, change)
    system.sum_current_lambda += change
    system.update_current_lambda += 1
    if system.update_current_lambda > 1e5
        system.sum_current_lambda = sum(system.current_lambda)
        system.update_current_lambda = 0
    end
end

"""
#evaluate current rates of SIR system
"""

function current_rates(system::SIR)
    rate_recovery = system.mu * system.I
    rate_infection = system.sum_current_lambda* system.S/system.N + system.epsilon
    return MVector{2,Float64}(rate_recovery, rate_infection)
end
"""


text:
using GenerativeAD
using FileIO, BSON
using ValueHistories, DistributionsAD
using Flux
using ConditionalDists
using GenerativeModels
using EvalMetrics
using Plots
using Statistics
using DrWatson

dataset = "arrhythmia"
dataset = "wall-following-robot"
dataset = "yeast"
dataset = "letter-recognition"
dataset = "kdd99_small"
seed=1
outpath = datadir("normal_validation/tst_auc_vs_val_score")
mkpath(outpath)
masterpath = datadir("experiments/tabular/")
models = [
	"aae_full", "adVAE", "GANomaly", "vae_full", "wae_full", "abod", "hbos", "if", "knn", "loda", "lof",
	"ocsvm_rbf", "ocsvm", "pidforest", "MAF", "RealNVP", "sptn", "fmgan", "gan", "MO_GAAL", "DeepSVDD",
	"vae_knn", "vae_ocsvm"
	]

"""
Returns AUC and mean anomaly score of normal validation samples.
"""
function auc_val_score(f)
	data = load(f)
	try
		auc = EvalMetrics.auc_trapezoidal(EvalMetrics.roccurve(data[:tst_labels], vec(data[:tst_scores]))...)
		val_score = Flux.mean(data[:val_scores][data[:val_labels].==0])
		return auc, val_score
	catch e
		return NaN, NaN
	end
end

# loop over all models, create plots and compute correlation
function explore_scores(model, dataset, seed)
	sfs = readdir(joinpath(masterpath, model, dataset, "seed=$seed"), join=true)
	# filter out modelfiles
	sfs = filter(x->!(occursin("model", x)),sfs)
	# filter only sampled reconstruction scores
	if any(occursin.("score=reconstruction-sampled", sfs))
		sfs = filter(x->occursin("reconstruction-sampled", x),sfs)
	end
	res = map(auc_val_score, sfs)
	aucs = [x[1] for x in res]
	val_scores = [x[2] for x in res]
	valid_inds = .!isnan.(val_scores)
	aucs, val_scores = aucs[valid_inds], val_scores[valid_inds]
	cor = round(Statistics.cor(val_scores, aucs), digits=3)
	scatter(val_scores, aucs, xlabel="validation anomaly score - normal data", ylabel="test AUC", title="$model - $dataset, correlation=$(cor)")
	sf = joinpath(outpath, "$(dataset)_$(model)_val_score_tst_auc_cor=$(cor).png")
	savefig(sf)
	@info "Saved $sf."
	return aucs, val_scores
end

for model in models
	@info "Processing $model"
	try
		explore_scores(model, dataset, seed)
	catch e
		@warn "Model $model not processed"
	end
end


text:
struct Broyden
end

macro broydentrace(stepnorm)
    esc(quote
        if tracing
            dt = Dict()
            if extended_trace
                dt["x"] = copy(x)
                dt["f(x)"] = copy(value(df))
                dt["g(x)"] = copy(jacobian(df))
            end
            update!(tr,
                    it,
                    maximum(abs, value(df)),
                    $stepnorm,
                    dt,
                    store_trace,
                    show_trace)
        end
    end)
end

function broyden_(df::Union{NonDifferentiable, OnceDifferentiable},
                    initial_x::AbstractArray{T},
                    xtol::T,
                    ftol::T,
                    iterations::Integer,
                    store_trace::Bool,
                    show_trace::Bool,
                    extended_trace::Bool,
                    linesearch) where T
    # setup
    n = length(initial_x)
    x = vec(copy(initial_x))

    value!!(df, x)

    vecvalue = vec(value(df))
    fold, xold = similar(vecvalue), similar(x)
    fold .= T(0)
    xold .= T(0)

    p = similar(x)
    g = similar(x)
    Jinv = Matrix{T}(I, n, n)
    check_isfinite(value(df))
    it = 0
    x_converged, f_converged, converged = assess_convergence(value(df), ftol)

    # FIXME: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
  #  lsr = LineSearches.LineSearchResults(T)

    tr = SolverTrace()
    tracing = store_trace || show_trace || extended_trace
    @broydentrace T(NaN)

    maybe_stuck = false
    max_resets = 3
    resets = 0
    while !converged && it < iterations

        it += 1

        copyto!(xold, x)
        copyto!(fold, value(df))
        p = - Jinv*fold

        ρ = T(0.9)
        σ₂ = T(0.001)
        x = xold + p
        value!(df, x)
        if norm(value(df), 2) <= ρ*norm(fold, 2) - σ₂*norm(p, 2)^2 # condition 2.7
            α = T(1.0)
        else
            α = approximate_norm_descent(x->value!(df, x), x, p)
            x = xold + α*p
        end

        value!(df, x)

        Δx = x-xold
        Δf = vecvalue .- fold

        maybe_stuck = all(abs.(Δx) .<= 1e-12) || all(abs.(Δf) .<= 1e-12)

        if maybe_stuck
            Jinv = Matrix{T}(I, n, n)
            resets += 1
            if resets > max_resets
                maybe_stuck = false # to allow convergence check
            end
        else
            Jinv = Jinv + ((Δx - Jinv *Δf)/(Δx'Jinv*Δf))*Δx'Jinv
        end

        if !maybe_stuck
            x_converged, f_converged, converged = assess_convergence(x, xold, value(df), xtol, ftol)
        end

        maybe_stuck = false
        @broydentrace sqeuclidean(x, xold)
    end
    return SolverResults("broyden without line-search",
                         initial_x, copyto!(similar(initial_x), x), norm(value(df), Inf),
                         it, x_converged, xtol, f_converged, ftol, tr,
                         first(df.f_calls), 0)
end

function broyden(df::Union{NonDifferentiable, OnceDifferentiable},
                   initial_x::AbstractArray{T},
                   xtol::Real,
                   ftol::Real,
                   iterations::Integer,
                   store_trace::Bool,
                   show_trace::Bool,
                   extended_trace::Bool,
                   linesearch) where T
    broyden_(df, initial_x, convert(T, xtol), convert(T, ftol), iterations, store_trace, show_trace, extended_trace, linesearch)
end

# A derivative-free line search and global convergence
# of Broyden-like method for nonlinear equations
# by Dong-Hui Li & Masao Fukushima
# https://doi.org/10.1080/10556780008805782
function approximate_norm_descent(F, x::AbstractArray{T}, p;
                                  lambda_0=T(1.0), beta=T(0.5),
                                  sigma_1=T(0.001), eta=T(0.1),
                                  nan_max_iter=5, max_iter=50) where T
    β, σ₁, η = beta, sigma_1, eta
    λ₂ = lambda_0
    λ₁ = λ₂

    # Calculate norms
    Fx_norm = norm(F(x), 2)
    if isnan(Fx_norm)
        throw(ErrorException("Exiting line search: 2-norm of current residuals is not a number. Value is: ", Fx_norm))
    elseif !isfinite(Fx_norm)
        throw(ErrorException("Exiting line search: 2-norm of current residuals is not finite. Value is: ", Fx_norm))
    end

    Fxλp_norm = norm(F(x - λ₂*p), 2)

    # If the step produces something that is not finite or a number, backtrack
    # a maximum number of nan_max_iter number of times to try to find a suitable
    # range of values for which F(x + λp) is well-behaved
    if isnan(Fxλp_norm) || !isfinite(Fxλp_norm)
         λ₁, λ₂ = nantrack(F, x, λ₂, p; β=β, nan_max_iter=nan_max_iter)
    end

    j = 0
    converged = norm(F(x + λ₂*p), 2) <= Fx_norm - σ₁*norm(λ₂*p, 2)^2 + η*norm(F(x), 2)

    while j < max_iter && !converged
        j += 1

        λ₁, λ₂ = λ₂, β*λ₂
        Fxλp_norm = norm(F(x - λ₂*p), 2)

        converged = norm(F(x + λ₂*p), 2) <= Fx_norm - σ₁*norm(λ₂*p, 2)^2 + η*norm(F(x), 2)
    end
    if j >= max_iter && !converged
        throw(ErrorException("Exiting line search: failed to satisfy condition (2.4)."))
    end
    λ₂
end

function nantrack(F, x, λ₂, p; β=0.5, nan_max_iter=5)
    i = 0
    nan_converged = false
    while i  <= nan_max_iter && !nan_converged
        i +=1

        λ₁, λ₂ = λ₂, β*λ₂
        Fxλp_norm = norm(F(x - λ₂*p), 2)

        converged = !isnan(Fxλp_norm)
    end
    if i >= nan_max_iter
        throw(ErrorException("Exiting line search: failed to find finite 2-norm of residuals in $nan_max_iter trials."))
    end
    return λ₁, λ₂
end


text:
using Test
using AbstractAlgebra
using Groups

using LinearAlgebra

@testset "Groups" begin

   @testset "generate balls" begin
      M = MatrixAlgebra(zz, 3)
      w = one(M); w[1,2] = 1;
      r = one(M); r[2,3] = -3;
      s = one(M); s[1,3] = 2; s[3,2] = -1;

      S = [w,r,s]; S = unique([S; inv.(S)]);
      _, sizes = Groups.generate_balls(S, radius=4);
      @test sizes == [7, 33, 141, 561]
   end

   include("FreeGroup-tests.jl")
   include("AutGroup-tests.jl")
   include("DirectPower-tests.jl")
   include("WreathProd-tests.jl")
end


text:
using Multigraphs, Graphs, SparseArrays
try
    m2 = spzeros(Int, 2, 3)
    dg = Multigraph(m2)
catch e
    @test e !== nothing
end
try
    m2 = spzeros(Int, 2, 2)
    m2[1, 2] = 2
    dg = Multigraph(m2)
catch e
    @test e !== nothing
end
try
    m2 = spzeros(Int, 2, 2)
    m2[1, 2] = -1
    m2[2, 1] = -1
    dg = Multigraph(m2)
catch e
    @test e !== nothing
end

m = spzeros(Int, 4, 4)
m
m[1,2] = 2
m[2,1] = 2
m[2,3] = 2
m[3,2] = 2
m[3,4] += 1
m[3,4] = 0
m[4,3] += 1
m[4,3] = 0
g = Multigraph(m)
g = Multigraph(Matrix(m))
g0 = Multigraph(2)
@test !add_edge!(g0, 2, 3) && !rem_edge!(g0, 1, 2)
g1 = Multigraph(path_graph(3))
@test adjacency_matrix(g) == m

@test !is_directed(g)
@test edgetype(g) == MultipleEdge{Int, Int}
@test size(adjacency_matrix(g), 1) == 4

@test nv(g) == 4 && ne(g, count_mul = true) == 4 && ne(g) == 2

add_vertices!(g, 3)
@test nv(g) == 7

@test has_edge(g, 1, 2, 2)
@test rem_vertices!(g, [7, 5, 4, 6])
add_edge!(g, [2, 3, 2])
rem_edge!(g, [2, 3, 2])
add_edge!(g, 2, 3)
rem_edge!(g, 2, 3)
add_edge!(g, 2, 3, 2)
rem_edge!(g, 2, 3, 1)

@test has_edge(g, 2, 3) && has_edge(g, [2, 3])
@test has_edge(g, 2, 3, 2) && has_edge(g, (2, 3, 2))
@test !has_edge(g, 2, 2) && !has_edge(g, 2, 5)
@test has_vertex(g, 1) && !has_vertex(g, 5)
for v in vertices(g)
    @test inneighbors(g, v) == outneighbors(g, v)
    @test degree(g, v) == indegree(g, v) && indegree(g, v) == outdegree(g, v)
end
add_vertex!(g)
@test indegree(g) == outdegree(g)

mg0 = Multigraph(0)
@test nv(mg0) == ne(mg0) == 0


text:
# Demonstrate the simulateMS function
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import RSDeltaSigmaPort: fft
import RSDeltaSigmaPort: BoundingBox
import Printf: @sprintf
j=im

#==Demo5Module: Module to define types, etc
===============================================================================#
module Demo5Module

struct SimConfig
	A::Float64
	f::Float64
	mtf
	dither::Float64
#	dw
	id::String
end

struct SimResult
	fin::Int #Input frequency
	sv::Array
	Svv::Array
	Sdd::Array
	leg::String
end

end #Demo5Module

import .Demo5Module: SimConfig, SimResult

#==
===============================================================================#
show_usage = false
OSR = 25
M = 16
N = 2^14
sigma_d = 0.01 #1% mismatch

dsm2b = RealDSM(order=2, OSR=OSR, M=M, opt=1, Hinf=2) #Second-order shaping
dsm4 = RealDSM(order=4, OSR=round(Int, OSR*0.9), M=M, opt=1, Hinf=1.3) #Fourth-order shaping
dsm6 = RealDSM(order=6, OSR=OSR, M=M, opt=1, Hinf=4)

mtf1 = _zpk(1,0,1,1)         #First-order shaping
mtf2 = _zpk([ 1 1 ], [ 0.3 0.3 ], 1, 1) #Second-order shaping
mtf2b = synthesizeNTF(dsm2b) #Second-order shaping
mtf4 = synthesizeNTF(dsm4)   #Fourth-order shaping

cases = [
	SimConfig(undbv(-3),  0.01,    [],   0, "Thermometer"),
	SimConfig(undbv(-3),  0.01,  mtf1,   0, "Rotation"),
	SimConfig(undbv(-3),  0.01,  mtf2,   0, "2^{nd}-order"),
	SimConfig(undbv(-30), 0.01,  mtf1,   0, "Rotation"),
	SimConfig(undbv(-30), 0.01,  mtf1, 0.5, "Rot + dither"),
	SimConfig(undbv(-3),  0.01, mtf2b,   0, "2^{nd}-order with zero"),
	SimConfig(undbv(-3),  0.01,  mtf4,   0, "4^{th}-order"),
]
comparisons = [ [1,2], [2,3], [4,5], [3,6], [6,7] ]

if false #Debug
	cases = cases[1:2]
	comparisons = comparisons[1:1]
end

nconfig = length(cases)
println("\nNumber of configurations to simulate: $nconfig"); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------

NTF = synthesizeNTF(dsm6)
window = ds_hann(N)'/(M*N/8)
windowM = repeat(window, M)

println()
@info("Running simulations..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
simresults = SimResult[] #Simulation results
for cfg in cases
	local u #Avoid complaints wrt "global u"
	local sv
	println("\t", cfg.id); flush(stdout); flush(stderr)
	fin = round(Int, cfg.f*N)
	inband = setdiff(2:1 .+ ceil(Int, 0.5*N/OSR), 1 .+ [0 1 fin .+ [-1 0 1]])
	w = (2π/N)*fin
	u = M*cfg.A*sin.(w*(0:N-1))
	v = simulateDSM(u, NTF, nlev=M+1).v #M unit elements requires an (M+1)-level quant.
	data = applywnd(v, window)
	Svv = abs.(fft(applywnd(v, window))) .^ 2
	if isa(cfg.mtf, Array) && isempty(cfg.mtf)
		sv = ds_therm(v, M)
	else
		sv = simulateMS(v, M=M, mtf=cfg.mtf, d=cfg.dither).sv
	end
	Sdd = sigma_d^2 * sum( abs.(fft(applywnd(sv, windowM),2)) .^ 2 , dims=1)
	mnp = sum(Sdd[inband])/1.5
	leg = @sprintf("%s (MNP= %.0f dBFS)", cfg.id, dbp(mnp))
	push!(simresults, SimResult(fin, sv, Svv, Sdd, leg))
end
println("\tdone.")

println()
@info("Plotting results..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
for case_nums in comparisons
	local plot #Avoid complaints wrt "global plot"
	nc = length(case_nums)
	#_cases = cases[case_nums]
	_results = simresults[case_nums]
	id1 = _results[1].leg; id2 = _results[2].leg #Only print out first 2
	println("Mismatch-Shaping Unit-Element DAC")
	println("\tComparing $id1 vs. $id2")
	if show_usage
		pcoll = cons(:plot_collection, ncolumns=1, title="Element Usage")
		T = 25
		for resi in _results
			plot = plotUsage(resi.sv[:,1:T])
			push!(pcoll, plot)
		end
		displaygui(pcoll)
	end

	colors = [:blue, :magenta, :red, :green, :orange]
	plot = plotSpectrum()
#		plot.title = "Error Spectra"
	set(plot,
		xyaxes=set(xmin=1e-3, xmax=0.5, ymin=-140, ymax=-50),
		labels=set(yaxis="Error PSD"),
	)
	for (i, resi) in enumerate(_results)
		plotSpectrum!(plot, sqrt.(resi.Sdd), resi.fin, lw=2, color=colors[i], id=resi.leg, n=4)
	end
	resi = _results[end]
	A = dbp(resi.Svv[resi.fin])
	plotSpectrum!(plot, sqrt.(resi.Svv), resi.fin, color=:green, n=5)
	lwfrm = waveform([1e-3, 0.5/OSR], -140*[1, 1])
	nbw_str = @sprintf("NBW = %.1e", 1.5/N) #orig. coords: (0.5, -140)
	push!(plot,
		cons(:wfrm, lwfrm, line=set(style=:solid, color=:black, width=8)),
		cons(:atext, nbw_str, y=-140, reloffset=set(x=0.95), align=:br),
	)
	plot.title = @sprintf("A = %.0fdBFS", A)
	displaygui(plot)
end

if true #Quadrature example
	println()
	@warn("TODO: Quadrature example")
end

:END_OF_DEMO


text:
using FractalLattices
using Test

@testset "FractalLattices.jl" begin
    # Write your tests here.
end


text:
@doc doc"""
`Q` state clock model with energy $E = -\sum_{ij} J_{ij} \cos(\theta_i - \theta_j)$,
where $\theta_i = 2\pi \sigma_i/Q$ and $\sigma_i$ takes an integer value from $1$ to $Q$.
"""
mutable struct Clock <: Model
    lat :: Lattice
    Q :: Int
    spins :: Matrix{Int}
    cosines :: Vector{Float64}
    sines :: Vector{Float64}
    sines_sw :: Vector{Float64}
    rng :: Random.MersenneTwister

    function Clock(lat::Lattice, Q::Integer, rng::Random.AbstractRNG)
        if Q < 2
            error("Q should be 2 or more: Q = $Q")
        end
        spins = rand(rng, 1:Q, 1, numsites(lat))
        cosines = [cospi(2s/Q) for s in 1:Q]
        sines = [sinpi(2s/Q) for s in 1:Q]
        sines_sw = [sinpi(2(s-0.5)/Q) for s in 1:Q]
        return new(lat, Q, spins, cosines, sines, sines_sw, rng)
    end
end
Clock(lat::Lattice, Q::Integer) = Clock(lat, Q, Random.seed!(Random.MersenneTwister(0)))
Clock(lat::Lattice, Q::Integer, seed) = Clock(lat, Q, Random.seed!(Random.MersenneTwister(0), seed...))

@doc doc"""
    Clock(param)

Generates `Clock` using `param["Lattice"]`, `param["Q"]`,  and `param["Seed"]` (if defined).
Each spin $\sigma_i$ will be initialized randomly and independently.
"""
function Clock(param::Parameter)
    lat = generatelattice(param)
    Q = param["Q"]
    if "Seed" in keys(param)
        return Clock(lat, Q, param["Seed"])
    else
        return Clock(lat, Q)
    end
end

include("update.jl")
include("estimator.jl")
include("postproc.jl")


text:
#=

```shell
julia make.jl
jupyter nbconvert --to notebook --execute --inplace example.ipynb
```

The first line converts src.jl to a notebook.
The second line runs the notebook.

If you want to save that in another format do the following

jupyter nbconvert --to html     example.ipynb
jupyter nbconvert --to webpdf   example.ipynb
jupyter nbconvert --to markdown example.ipynb
```

=#

using Literate

config = Dict(
    "name"          => "example",
    "documenter"    => false,
    "keep_comments" => false,
    "execute"       => false,
    "credit"        => false,
)

Literate.notebook("src.jl"; config=config)


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0, 1f0,360))
end
description = """
You can simply visualize a range with the particle system. The positions
will be placed on a 1D grid.
"""

# last argument can be used to control the granularity of the resulting mesh
sphere = GLNormalMesh(Sphere(Point3f0(0.5), 0.5f0), 24)
c = collect(linspace(0.1f0,1.0f0,10f0))
rotation = map(rotationmatrix_z, const_lift(*, timesignal, 2f0*pi))
# create a visualisation
vis = visualize((sphere, c))
_view(vis, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
function query_expression_translation_phase_A(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if i==1 && clause.head==:macrocall && clause.args[1]==Symbol("@from")
			if !(isa(clause.args[2].args[3], Expr) && clause.args[2].args[3].head==:macrocall && isa(clause.args[2].args[3].args[1],Expr) && clause.args[2].args[3].args[1].head==:. && clause.args[2].args[3].args[1].args[1]==:Query)
				clause.args[2].args[3] = :(Query.query($(esc(clause.args[2].args[3]))))
			end
		elseif clause.head==:macrocall && clause.args[1]==Symbol("@from")
			clause.args[2].args[3] = :(Query.query($(clause.args[2].args[3])))
		elseif clause.head==:macrocall && clause.args[1]==Symbol("@join")
			clause.args[2].args[3] = :(Query.query($(esc(clause.args[2].args[3]))))
		elseif clause.head==:macrocall && clause.args[1]==Symbol("@select")
			if isa(clause.args[2], Expr) && clause.args[2].head==:cell1d
				clause.args[2] = Expr(:macrocall, Symbol("@NT"), clause.args[2].args...)
				for (j,field_in_NT) in enumerate(clause.args[2].args[2:end])
					if isa(field_in_NT, Expr) && field_in_NT.head==:(=)
						clause.args[2].args[j+1] = Expr(:(=>), field_in_NT.args...)
					elseif isa(field_in_NT, Expr) && field_in_NT.head==:.
						name_to_use = field_in_NT.args[2].args[1]
						clause.args[2].args[j+1] = Expr(:(=>), name_to_use, field_in_NT)
					elseif isa(field_in_NT, Symbol)
						clause.args[2].args[j+1] = Expr(:(=>), field_in_NT, field_in_NT)
					end
				end
			end
		end
		i+=1
	end
end

function query_expression_translation_phase_1(qe)
	done = false
	while !done
		group_into_index = findfirst(i->i.head==:macrocall && i.args[1]==Symbol("@group") && length(i.args)==6 && i.args[5]==:into,qe)
		select_into_index = findfirst(i->i.head==:macrocall && i.args[1]==Symbol("@select") && length(i.args)==4 && i.args[3]==:into,qe)
		if length(qe)>=2 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && group_into_index>0
			x = qe[group_into_index].args[6]

			sub_query = Expr(:block, qe[1:group_into_index]...)
			deleteat!(sub_query.args[end].args,6)
			deleteat!(sub_query.args[end].args,5)

			translate_query(sub_query)

			if length(sub_query.args)>1
				error("Subquery too long")
			end

			qe[1] = :( @from $x in $(sub_query.args[1]) )
			for i=group_into_index:-1:2
				deleteat!(qe,i)
			end
		elseif length(qe)>=2 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && select_into_index>0
			x = qe[select_into_index].args[4]

			sub_query = Expr(:block, qe[1:select_into_index]...)
			deleteat!(sub_query.args[end].args,4)
			deleteat!(sub_query.args[end].args,3)

			translate_query(sub_query)

			if length(sub_query.args)>1
				error("Subquery too long")
			end

			qe[1] = :( @from $x in $(sub_query.args[1]) )
			for i=select_into_index:-1:2
				deleteat!(qe,i)
			end
		else
			done = true
		end
	end
end

function query_expression_translation_phase_3(qe)
	done = false
	while !done
		if length(qe)>=2 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@select") && qe[1].args[2].args[2]==qe[2].args[2]
			x = qe[1].args[2].args[2]
			e = qe[1].args[2].args[3]

			qe[1] = :( Query.@select_internal($e,x->x) )
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

function query_expression_translation_phase_4(qe)
	done = false
	while !done
		if length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@from") && qe[3].head==:macrocall && qe[3].args[1]==Symbol("@select")
			x1 = qe[1].args[2].args[2]
			x2 = qe[2].args[2].args[2]
			e1 = qe[1].args[2].args[3]
			e2 = qe[2].args[2].args[3]
			v = qe[3].args[2]

			f_collection_selector = Expr(:->, x1, e2)
			f_result_selector = Expr(:->, Expr(:tuple,x1,x2), v)

			qe[1] = :( Query.@select_many_internal($e1, $(esc(f_collection_selector)), $(esc(f_result_selector))) )
			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@from")
			x1 = qe[1].args[2].args[2]
			x2 = qe[2].args[2].args[2]
			e1 = qe[1].args[2].args[3]
			e2 = qe[2].args[2].args[3]

			f_collection_selector = Expr(:->, x1, e2)
			f_result_selector = Expr(:->, Expr(:tuple,x1,x2), :(@NT($x1=>$x1,$x2=>$x2)))

			qe[1].args[2].args[2] = Expr(:transparentidentifier, gensym(:t), x1, x2)
			qe[1].args[2].args[3] = :( Query.@select_many_internal($e1, $(esc(f_collection_selector)), $(esc(f_result_selector))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@let")
			x = qe[1].args[2].args[2]
			e = qe[1].args[2].args[3]
			y = qe[2].args[2].args[1]
			f = qe[2].args[2].args[2]

			f_selector = Expr(:->, x, :(@NT($x=>$x,$y=>$f)))

			qe[1].args[2].args[2] = Expr(:transparentidentifier, gensym(:t), x, y)
			qe[1].args[2].args[3] = :( Query.@select_internal($e,$(esc(f_selector))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@where")
			x = qe[1].args[2].args[2]
			e = qe[1].args[2].args[3]
			f = qe[2].args[2]

			f_condition = Expr(:->, x, f)

			qe[1].args[2].args[3] = :( Query.@where_internal($e,$(esc(f_condition))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@join") && length(qe[2].args)==6 && qe[3].head==:macrocall && qe[3].args[1]==Symbol("@select")
			outer = qe[1].args[2].args[3]
			inner = qe[2].args[2].args[3]
			outer_range_var = qe[1].args[2].args[2]
			inner_range_var = qe[2].args[2].args[2]
			f_outer_key = Expr(:->, outer_range_var, qe[2].args[4])
			f_inner_key = Expr(:->, inner_range_var, qe[2].args[6])
			f_result = Expr(:->, Expr(:tuple,outer_range_var,inner_range_var), qe[3].args[2])
			qe[1] = :(
				Query.@join_internal($outer, $inner, $(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result)))
				)

			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@join") && length(qe[2].args)==6
			e1 = qe[1].args[2].args[3]
			e2 = qe[2].args[2].args[3]
			x1 = qe[1].args[2].args[2]
			x2 = qe[2].args[2].args[2]
			k1 = qe[2].args[4]
			k2 = qe[2].args[6]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,x2), :(@NT($x1=>$x1,$x2=>$x2)) )

			qe[1].args[2].args[2] = Expr(:transparentidentifier, gensym(:t), x1, x2)
			qe[1].args[2].args[3] = :( Query.@join_internal($e1,$e2,$(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@join") && length(qe[2].args)==8 && qe[3].head==:macrocall && qe[3].args[1]==Symbol("@select")
			e1 = qe[1].args[2].args[3]
			e2 = qe[2].args[2].args[3]
			x1 = qe[1].args[2].args[2]
			x2 = qe[2].args[2].args[2]
			k1 = qe[2].args[4]
			k2 = qe[2].args[6]
			g = qe[2].args[8]
			v = qe[3].args[2]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,g), v)
			qe[1] = :( Query.@group_join_internal($e1, $e2, $(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )

			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@join") && length(qe[2].args)==8
			e1 = qe[1].args[2].args[3]
			e2 = qe[2].args[2].args[3]
			x1 = qe[1].args[2].args[2]
			x2 = qe[2].args[2].args[2]
			k1 = qe[2].args[4]
			k2 = qe[2].args[6]
			g = qe[2].args[8]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,g), :(@NT($x1=>$x1,$g=>$g)) )

			qe[1].args[2].args[2] = Expr(:transparentidentifier, gensym(:t), x1, g)
			qe[1].args[2].args[3] = :( Query.@group_join_internal($e1,$e2,$(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@orderby")
			e = qe[1].args[2].args[3]
			x = qe[1].args[2].args[2]
			ks = []
			if isa(qe[2].args[2], Expr) && qe[2].args[2].head==:tuple
				for sort_clause in qe[2].args[2].args
					if isa(sort_clause, Expr) && sort_clause.head==:call && sort_clause.args[1]==:descending
						k = sort_clause.args[2]
						direction = :descending
					elseif isa(sort_clause, Expr) && sort_clause.head==:call && sort_clause.args[1]==:ascending
						k = sort_clause.args[2]
						direction = :ascending
					else
						k = sort_clause
						direction = :ascending
					end
					push!(ks, (k, direction))
				end
			else
				if isa(qe[2].args[2], Expr) && qe[2].args[2].head==:call && qe[2].args[2].args[1]==:descending
					k = qe[2].args[2].args[2]
					direction = :descending
				elseif isa(qe[2].args[2], Expr) && qe[2].args[2].head==:call && qe[2].args[2].args[1]==:ascending
					k = qe[2].args[2].args[2]
					direction = :ascending
				else
					k = qe[2].args[2]
					direction = :ascending
				end
				push!(ks, (k, direction))
			end

			for (i,sort_clause) in enumerate(ks)
				f_condition = Expr(:->, x, sort_clause[1])

				if sort_clause[2]==:ascending
					if i==1
						qe[1].args[2].args[3] = :( Query.@orderby_internal($e,$(esc(f_condition))) )
					else
						qe[1].args[2].args[3] = :( Query.@thenby_internal($(qe[1].args[2].args[3]),$(esc(f_condition))) )
					end
				elseif sort_clause[2]==:descending
					if i==1
						qe[1].args[2].args[3] = :( Query.@orderby_descending_internal($e,$(esc(f_condition))) )
					else
						qe[1].args[2].args[3] = :( Query.@thenby_descending_internal($(qe[1].args[2].args[3]),$(esc(f_condition))) )
					end
				end
			end
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

function query_expression_translation_phase_5(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if clause.head==:macrocall && clause.args[1]==Symbol("@select")
			from_clause = qe[i-1]
			if from_clause.head!=:macrocall || from_clause.args[1]!=Symbol("@from")
				error("Error in phase 5")
			end
			range_var = from_clause.args[2].args[2]
			source = from_clause.args[2].args[3]
			if clause.args[2]==range_var
				qe[i-1] = source
			else
				func_call = Expr(:->, range_var, clause.args[2])
				qe[i-1] = :( Query.@select_internal($source, $(esc(func_call))) )
			end
			deleteat!(qe,i)
		else
			i+=1
		end
	end
end

function query_expression_translation_phase_6(qe)
	done = false
	while !done
		if length(qe)>=2 && qe[1].head==:macrocall && qe[1].args[1]==Symbol("@from") && qe[2].head==:macrocall && qe[2].args[1]==Symbol("@group")
			e = qe[1].args[2].args[3]
			x = qe[1].args[2].args[2]
			v = qe[2].args[2]
			k = qe[2].args[4]

			f_elementSelector = Expr(:->, x, k)
			f_resultSelector = Expr(:->, x, v)

			if v==x
				qe[1] = :( Query.@group_by_internal_simple($e, $(esc(f_elementSelector))) )
			else
				qe[1] = :( Query.@group_by_internal($e, $(esc(f_elementSelector)), $(esc(f_resultSelector))) )
			end
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

# Phase 7

function replace_transparent_identifier_in_anonym_func(ex::Expr, names_to_put_in_scope)
	for (i,child_ex) in enumerate(ex.args)
		if isa(child_ex, Expr)
			replace_transparent_identifier_in_anonym_func(child_ex, names_to_put_in_scope)
		elseif isa(child_ex, Symbol)
			index_of_name = findfirst(j->child_ex==j[2], names_to_put_in_scope)
			if index_of_name>0 && !(ex.head==Symbol("=>") && i==1)
				ex.args[i] = Expr(:., names_to_put_in_scope[index_of_name][1], QuoteNode(child_ex))
			end
		end
	end
end

function find_names_to_put_in_scope(ex::Expr)
	names = []
	for child_ex in ex.args[2:end]
		if isa(child_ex,Expr) && child_ex.head==:transparentidentifier
			child_names = find_names_to_put_in_scope(child_ex)
			for child_name in child_names
				push!(names, (Expr(:., ex.args[1], QuoteNode(child_name[1])), child_name[2]))
			end
		elseif isa(child_ex, Symbol)
			push!(names,(ex.args[1],child_ex))
		elseif isa(child_ex, Expr) && child_ex.head==:.
			push!(names,(ex.args[1],child_ex.args[2].value))
		else
			error()
		end
	end
	return names
end

function find_and_translate_transparent_identifier(ex::Expr)
	# First expand any transparent identifiers in lambdas
	if ex.head==:-> && isa(ex.args[1], Expr) && ex.args[1].head==:transparentidentifier
		names_to_put_in_scope = find_names_to_put_in_scope(ex.args[1])
		ex.args[1] = ex.args[1].args[1]
		replace_transparent_identifier_in_anonym_func(ex, names_to_put_in_scope)
	end

	for (i,child_ex) in enumerate(ex.args)
		if isa(child_ex, Expr) && child_ex.head==:transparentidentifier
			ex.args[i] = child_ex.args[1]
		elseif isa(child_ex, Expr)
			find_and_translate_transparent_identifier(child_ex)
		end
	end
end

function query_expression_translation_phase_7(qe)
	for clause in qe
		if isa(clause, Expr)
			find_and_translate_transparent_identifier(clause)
		end
	end
end

function query_expression_translation_phase_B(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if clause.head==:macrocall && clause.args[1]==Symbol("@collect")
			previous_clause = qe[i-1]
			if length(clause.args)==1
				qe[i-1] = :( collect($previous_clause) )
			else
				qe[i-1] = :( collect($previous_clause, $(esc(clause.args[2]))) )
			end
			deleteat!(qe,i)
		else
			i+=1
		end
	end
end

function translate_query(body)
	debug_output = false

	debug_output && println("AT START")
	debug_output && println(body)

	query_expression_translation_phase_1(body.args)
	debug_output && println("AFTER 1")
	debug_output && println(body)

	query_expression_translation_phase_A(body.args)
	debug_output && println("AFTER A")
	debug_output && println(body)

	query_expression_translation_phase_3(body.args)
	debug_output && println("AFTER 3")
	debug_output && println(body)

	query_expression_translation_phase_4(body.args)
	debug_output && println("AFTER 4")
	debug_output && println(body)

	query_expression_translation_phase_5(body.args)
	debug_output && println("AFTER 5")
	debug_output && println(body)

	query_expression_translation_phase_6(body.args)
	debug_output && println("AFTER 6")
	debug_output && println(body)

	query_expression_translation_phase_7(body.args)
	debug_output && println("AFTER 7")
	debug_output && println(body)

	query_expression_translation_phase_B(body.args)
	debug_output && println("AFTER B")
	debug_output && println(body)
end


text:
#####
### This is the specialized code I wrote to simplify the similarity
### data. There is otherwise too much variation to find patterns.
### We use binning to reduce values to categorical values.
### For most of the analysis we simply reduce it to whole
### percents.
#####

using DataFrames, DataArrays

"""
    find_bin(x::Float64, limits::Vector{Float64})
This is a pretty simple function. It associates a value with largest value
in the limits vector.
# Examples
```julia-repl
julia> find_bin(0.357, [0.0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
3
```
"""
function find_bin(x::Float64, limits::Vector{Float64})
    bin = length(limits) + 1
    for i in 1:length(limits)
        if x < limits[i]
            bin = i
            break
        end
    end
    bin
end

"""
    find_bin(x::Float64, limits::StepRangeLen{Float64})
This is a pretty simple function. It associates a value with largest value
in the limits range.
# Examples
```julia-repl
julia> find_bin(0.357, range(0.0,0.1,11))
3
julia> find_bin(0.357, 0.0:0.1:11))
3
```
"""
function find_bin(x::Float64, limits::StepRangeLen{Float64})
    bin = length(limits) + 1
    for i in 1:length(limits)
        if x < limits[i]
            bin = i
            break
        end
    end
    bin
end

"""
    find_bins(dv::DataFrames.DataVector{Float64}, limits::Vector{Float64})
Takes a DataVector, like a column in a DataFrame, and computes the bins based
on the given range.
"""
function find_bins(dv::DataFrames.DataVector{Float64}, limits::Vector{Float64})
    vbins = zeros(Int64, length(dv))
    for i in 1:length(dv)
        vbins[i] = find_bin(dv[i], limits)
    end
    vbins
end

"""
    find_bins(dv::DataFrames.DataVector{Float64}, limits::StepRangeLen{Float64})
Takes a DataVector, like a column in a DataFrame, and computes the bins based
on the given range.
"""
function find_bins(dv::DataFrames.DataVector{Float64}, limits::StepRangeLen{Float64})
    vbins = zeros(Int64, length(dv))
    for i in 1:length(dv)
        vbins[i] = find_bin(dv[i], limits)
    end
    vbins
end

"""
    group_bins{T}(df::AbstractDataFrame, col::T, limits::Vector{Float64})
This is currently not working. It's an attempt to do a groupby with the
bins, but its probably just better to incorportate the bins as a new
column and do the groupby on the updated data frame. My kungfu wasn't
strong enough to fix this yet.
"""
function group_bins{T}(df::AbstractDataFrame, col::T, limits::Vector{Float64})
    n_bins = length(limits) + 1

    vbins = find_bins(df[col], limits)
    (idx, starts) = DataArrays.groupsort_indexer(vbins, n_bins)

    # Remove zero-length groupings
    starts = _uniqueofsorted(starts)
    ends = [starts[2:end] - 1]
    GroupedDataFrame(df, [col], idx, starts[1:end-1], ends)
end

"""
    group_bins{T}(df::AbstractDataFrame, col::T, limits::StepRangeLen{Float64})
This is currently not working. It's an attempt to do a groupby with the
bins, but its probably just better to incorportate the bins as a new
column and do the groupby on the updated data frame. My kungfu wasn't
strong enough to fix this yet.
"""
function group_bins{T}(df::AbstractDataFrame, col::T, limits::StepRangeLen{Float64})
    n_bins = length(limits) + 1

    df[:idx] = find_bins(df[col], limits)
    return groupby(df, :idx)
    # (idx, starts) = DataArrays.groupsort_indexer(vbins, n_bins)
    #
    # # Remove NA groupings
    # starts = unique(starts)
    # ends = [starts[2:end] - 1]
    # #GroupedDataFrame(df, [col], idx, starts[1:end-1], ends)
    # groupby(df, [col], idx, starts[1:end-1], ends)
end

#####
### Here are some tests, they are useful in getting a feel for how binning
### operates.
# x_df = DataFrame(X=randn(30), Y=rand(30))
#
# find_bins(x_df[:X], Array(range(0.0, 0.5, 21)))
#####

#####
# the next file you should run is gapanalyisis_stay_change.jl
#####


text:
# Temporary fix... just to make tests pass

abstract type AbstractPriceData{Tprice,Tvol} end

const Price = Float64 #Nullable{Float64}
const Volume = Float64 #Nullable{Float64}

null_volume = NaN