text:
using Gadfly, StatsBase, DSP
import Gadfly.plot, StatsBase.pacf

colorscale = Scale.color_discrete_hue()

function buildlayers(f::Function, colnames, colors)
  [f(coln, colr) for (coln, colr) in zip(colnames, colors)]
end #buildlayers

# Raw data visualization

function plot(ta::TimeArray; subplots=false, title="")
  subplots ? error("Subplotting not yet supported") :
  linecolors = colorscale.f(length(ta.colnames))
  layers = buildlayers(ta.colnames, linecolors) do coln, colr
    layer(x=ta.timestamp, y=ta[coln].values, Geom.line, Theme(default_color=colr))
  end
  plot(layers..., Guide.xlabel(""), Guide.ylabel(""), Guide.title(title), Guide.manual_color_key("", ta.colnames, linecolors))
end #plot

function histogram(ta::TimeArray; subplots=false, title="")
  subplots ? error("Subplotting not yet supported") :
  linecolors = colorscale.f(length(ta.colnames))
  layers = buildlayers(ta.colnames, linecolors) do coln, colr
    layer(x=ta[coln].values, Geom.histogram, Theme(default_color=colr))
  end
  plot(layers..., Guide.xlabel(""), Guide.ylabel("Counts"), Guide.title(title), Guide.manual_color_key("", ta.colnames, linecolors))
end #hist

function density(ta::TimeArray; subplots=false, title="")
  subplots ? error("Subplotting not yet supported") :
  linecolors = colorscale.f(length(ta.colnames))
  layers = buildlayers(ta.colnames, linecolors) do coln, colr
    layer(x=ta[coln].values, Geom.density, Theme(default_color=colr))
  end
  plot(layers..., Guide.xlabel(""), Guide.ylabel("Probability Density"), Guide.title(title), Guide.manual_color_key("", ta.colnames, linecolors))
end # density

# Serial correlation & cross-correlation visualization

function directedscatter(ta::TimeArray; naive=false, subplots=false, nlags=1, title="")
  naive || require_clean(ta)
  subplots ? error("Subplotting not yet supported") :
  linecolors = colorscale.f(length(ta.colnames))
  layers = buildlayers(ta.colnames, linecolors) do coln, colr
    layer(x=lag(ta[coln], nlags).values, y=ta[coln].values[nlags+1:end], Geom.line(preserve_order=true), Theme(default_color=colr))
  end
  plot(layers..., Guide.xlabel("t-"*string(nlags)), Guide.ylabel("t"), Guide.title(title), Guide.manual_color_key("", ta.colnames, linecolors))
end # directedscatter

function laghistogram(ta::TimeArray; naive=false, nlags=1, title="")
  naive || require_clean(ta)
  length(ta.colnames) > 1 ? error("Subplotting not yet supported. Use a single-column TimeArray.") :
  plot(
    x=lag(ta, nlags).values, y=ta.values[nlags+1:end], Geom.histogram2d,
    Guide.xlabel("t-"*string(nlags)), Guide.ylabel("t"), Guide.title(title)
  )
end # laghistogram

function acf(ta::TimeArray; naive=false, lags=1:30)
  naive || require_clean(ta)
  length(ta.colnames) > 1 ? error("Subplotting not yet supported. Use a single-column TimeArray.") :
  conf_int = 1.96 / sqrt(length(ta))
  plot(x=lags, y=autocor(ta.values, lags), xmin=collect(lags)-.1, xmax=lags, yintercept=[-conf_int, +conf_int], Geom.bar, Geom.hline(color=colorant"orange"), Scale.y_continuous(minvalue=-1, maxvalue=1), Guide.xlabel("lags"), Guide.ylabel(""), Guide.title(ta.colnames[1]*" ACF"))
end # acf

function pacf(ta::TimeArray; naive=false, lags=1:30)
  naive || require_clean(ta)
  length(ta.colnames) > 1 ? error("Subplotting not yet supported. Use a single-column TimeArray.") :
  conf_int = 1.96 / sqrt(length(ta))
  plot(x=lags, y=pacf(ta.values, lags), xmin=collect(lags)-.1, xmax=lags, yintercept=[-conf_int, +conf_int], Geom.bar, Geom.hline(color=colorant"orange"), Scale.y_continuous(minvalue=-1, maxvalue=1), Guide.xlabel("lags"), Guide.ylabel(""), Guide.title(ta.colnames[1]*" PACF"))
end # pacf

function ccf(ta::TimeArray; naive=false, lags=-15:15)
  naive || require_clean(ta)
  length(ta.colnames) != 2 ? error("`ccf` requires a 2-column TimeArray") :
  conf_int = 1.96 / sqrt(length(ta))
  plot(x=lags, y=crosscor(ta[ta.colnames[1]].values, ta[ta.colnames[2]].values, lags), xmin=collect(lags)-.1, xmax=lags, yintercept=[-conf_int, +conf_int], Geom.bar, Geom.hline(color=colorant"orange"), Scale.y_continuous(minvalue=-1, maxvalue=1), Guide.xlabel("lags"), Guide.ylabel(""), Guide.title(ta.colnames[1]*", "*ta.colnames[2]*" CCF"))
end # ccf

# Spectral analysis

function cumulative_periodogram(ta::TimeArray; naive=false)
  naive || require_clean(ta)
  length(ta.colnames) > 1 ? error("Subplotting not yet supported. Use a single-column TimeArray.") :
  p = periodogram(ta.values)
  plot(x=freq(p), y=cumsum(power(p)), Geom.line, Guide.xlabel("Frequency"), Guide.ylabel("Cumulative Power"), Guide.title(ta.colnames[1]*" Cumulative Periodogram"))
end # cumulatve_periodogram


text:
#=

Example 1: Govt spending is AR(1) and state is (g, 1).

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-21

References
----------

Simple port of the file examples/lqramsey_ar1.py

http://quant-econ.net/lqramsey.html

=#
include("lqramsey.jl")

# == Parameters == #
bet = 1 / 1.05
rho, mg = .7, .35
A = eye(2)
A = [rho mg*(1 - rho); 0.0 1.0]
C = [sqrt(1 - rho^2)*mg/10 0.0]'
Sg = [1.0 0.0]
Sd = [0.0 0.0]
Sb = [0 2.135]
Ss = [0.0 0.0]
discrete = false
proc = ContStochProcess(A, C)

econ = Economy(bet, Sg, Sd, Sb, Ss, discrete, proc)
T = 50

path = compute_paths(econ, T)

gen_fig_1(path)

ext: jl

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
# It is suggested that test code for MLJ.jl include files be placed in
# a file of the same name under "test/" (and included below) and that
# this test code be wrapped in a module. Any new module name will do -
# eg, `module TestDatasets` for code testing `datasets.jl`.

# using Revise
using MLJBase
using Test

@testset "traits" begin
  @test include("traits.jl")
end

@testset "distributions" begin
  @test include("distributions.jl")
end

@testset "data" begin
  @test include("data.jl")
end


text:
using DataStructures
using Pipe

function parse_template(l)::Pair{Tuple{Char, Char}, Char}
    lhs, rhs = split(l, " -> ")
    Tuple(lhs) => rhs[1]
end

const word = readuntil(Base.stdin, "\n\n")
const templates = Dict(parse_template.(readlines()))

pairwise(xs) = zip(xs, last(Iterators.peel(xs)))

function step(polymer)
    new_polymer = []
    for (c1, c2) in pairwise(polymer)
        push!(new_polymer, c1)
        if haskey(templates, (c1, c2))
            push!(new_polymer, templates[(c1, c2)])
        end
    end
    push!(new_polymer, polymer[end])
    join(new_polymer)
end

function step(polymer, n::Int)
    for _ in 1:n
        polymer = step(polymer)
    end
    polymer
end

# Part One
@pipe step(word, 10) |>
      collect |>
      counter |>
      values |>
      maximum(_) - minimum(_) |>
      println

const all_chars = unique(values(templates))
const n_chars = length(all_chars)
const n_pairs = n_chars^2
const pairs = [(c1, c2) for c1 in all_chars for c2 in all_chars]

# Building a dp table n_pairs × n_pairs
# Each line provides what the current pair will produce after one step
# It is either itself (1 on the diagonal)
#       or 2 new pairs formed by inserting a char inbetween
const dp_mem = begin
    dp = zeros(Int, n_pairs, n_pairs)
    for i in 1:n_pairs
        if haskey(templates, pairs[i])
            c_inter = templates[pairs[i]]
            for t in ((pairs[i][1], c_inter), (c_inter, pairs[i][2]))
                j = findfirst(==(t), pairs)
                dp[i, j] = 1
            end
        else
            dp[i, i] = 1
        end
    end
    dp
end

function count_occurences(polymer, n::Int)
    # Exponentiations on `dp_mem` computes what a given pair produces after `n` steps
    # Julia automatically uses fast exponentiations <3 so we could do the same computation
    # for very large numbers (using BigInt instead of Int64)
    occurence_table = dp_mem^n

    char_count = zeros(Int, n_chars)
    pair_count = zeros(Int, n_pairs)

    # Counting pairs first
    for p in pairwise(polymer)
        i = findfirst(==(p), pairs)
        pair_count .+= occurence_table[i, :]
    end
    # Now counting chars
    for (i, n_p) in enumerate(pair_count)
        for c in pairs[i]
            char_count[findfirst(==(c), all_chars)] += n_p
        end
    end
    # Chars are counted twice (they all belong to 2 pairs)
    # Except for first and last char of polymer sequence
    char_count[findfirst(==(polymer[1]), all_chars)] += 1
    char_count[findfirst(==(polymer[end]), all_chars)] += 1
    char_count .÷ 2
end

# Part Two
@pipe count_occurences(word, 40) |>
      maximum(_) - minimum(_) |>
      println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using KernelDensity, Distributions, Interpolations
import KernelDensity: kernel_dist
import Distributions: pdf

kernel_dist(::Type{Epanechnikov}, w::Float64) = Epanechnikov(0.0, w)
kernel(data) = kde(data; kernel=Epanechnikov)

struct Prob <: ContinuousUnivariateDistribution
    θ::Float64
end

pdf(dist::Prob, data::Real) = dist.θ

function generate_data(;set_sizes, n_trials, kwargs...)
    experiments = map(x->Experiment(n_trials=n_trials,
        set_size=x, populate_visicon=conjunctive_set), set_sizes)
    map(x->run_condition!(x; kwargs...), experiments)
    return experiments
end

function condition_LL(condition; kwargs...)
    raw_preds = deepcopy(condition)
    raw_preds.n_trials = 10^3
    run_condition!(raw_preds; kwargs...)
    preds = generate_kde(raw_preds)
    LL = 0.0
    for d in condition.data
        dist = preds[d.trial_type]
        LL += max(log(pdf(dist, d.rt)), -1000.0)
    end
    return LL
end

function compute_LL(conditions; kwargs...)
    LL = 0.0
    for condition in conditions
        LL += condition_LL(condition; kwargs...)
    end
    return LL
end

generate_kde(experiment::Experiment) = generate_kde(experiment.data)

function generate_kde(data)
    f(present, response, x) = x.target_present==present && x.response==response
    hits = filter(x->f(:present, :present, x), data)
    misses = filter(x->f(:present, :absent, x), data)
    hits_rt = map(x->x.rt, hits)
    p_hit = length(hits)/(length(hits) + length(misses))
    hit_dens = kernel(hits_rt)
    hit_dens.density *= p_hit
    miss_dens = Prob(1-p_hit)
    # misses = filter(x->f(:present, :absent, x), data)
    # misses_rt = map(x->x.rt, misses)
    # if isempty(misses_rt)
    #     dens = kernel([0.0])
    #     dens.density *= (1-p_hit)
    #     preds[:miss] = dens
    # else
    #     dens = kernel(misses_rt)
    #     dens.density *= (1-p_hit)
    #     preds[:miss] = dens
    # end
    cr = filter(x->f(:absent, :absent, x), data)
    cr_rt = map(x->x.rt, cr)
    cr_dens = kernel(cr_rt)
    return (hit=hit_dens, miss=miss_dens, cr=cr_dens)
end

# function generate_kde(data)
#     T1 = typeof(kernel(rand(2)))
#     preds = Dict{Symbol,Union{Prob,T1}}()
#     f(present, response, x) = x.target_present==present && x.response==response
#     hits = filter(x->f(:present, :present, x), data)
#     misses = filter(x->f(:present, :absent, x), data)
#     hits_rt = map(x->x.rt, hits)
#     p_hit = length(hits)/(length(hits) + length(misses))
#     dens = kernel(hits_rt)
#     dens.density *= p_hit
#     preds[:hit] = dens
#     preds[:miss] = Prob(1-p_hit)
#     miss_dens = Prob(1-p_hit)
#     cr = filter(x->f(:absent, :absent, x), data)
#     cr_rt = map(x->x.rt, cr)
#     preds[:cr] = kernel(cr_rt)
#     cr_dens = kernel(cr_rt)
#     return preds
# end

max_stars_repo_path: Bayesian PDA/KDE.jl

max_issues_repo_path: Bayesian PDA/KDE.jl

max_forks_repo_path: Bayesian PDA/KDE.jl


text:
## General HTML rendering - static pages and IJulia –––––––––––––––––––––––––––––––––––––

function save(file::AbstractString, mime::MIME"text/html", doc::Metadata, config::Config)
    config.include_internal ||
            throw(ArgumentError("`config` option `include_internal` must be true for html"))
    ents = EntriesHtml()
    # Write the main file.
    isfile(file) || mkpath(dirname(file))
    open(file, "w") do f
        info("writing documentation to $(file)")
        headerhtml(f, doc, config)
        ents = writehtml(f, doc, ents, config)
        footerhtml(f, doc, config)
    end

    # copy static files
    src = joinpath(dirname(@__FILE__), "..", "..", "static")
    dst = joinpath(dirname(file), "static")
    isdir(dst) || mkpath(dst)
    for file in readdir(src)
        info("copying $(file) to $(dst)")
        cp(joinpath(src, file), joinpath(dst, file); remove_destination=true)
    end
    return ents
end

type EntriesHtml
    entries::Vector{@compat(Tuple{Module, Any, AbstractEntry})}
end
EntriesHtml() = EntriesHtml(@compat(Tuple{Module, Any, AbstractEntry})[])

function push!(ents::EntriesHtml, modulename::Module, obj, ent::AbstractEntry)
    push!(ents.entries, (modulename, obj, ent))
end

function writehtml(io::IO, doc::Metadata, ents::EntriesHtml, config::Config)
    rootdir = isfile(root(doc)) ? dirname(root(doc)) : root(doc)
    for file in manual(doc)
        writemime(io, MIME("text/html"), Markdown.parse(readall(joinpath(rootdir, file))))
    end

    index = Dict{Symbol, Any}()
    for (obj, entry) in entries(doc)
        addentry!(index, obj, entry)
    end

    if !isempty(index)
        println(io, "
Reference
")
        wrap(io, "ul", "class='index'") do
            for k in config.category_order
                haskey(index, k) || continue
                wrap(io, "li") do
                    println(io, "$(k)s:")
                end
                wrap(io, "li") do
                    wrap(io, "ul") do
                        for (s, obj) in index[k]
                            push!(ents, modulename(doc), obj, entries(doc)[obj])
                            wrap(io, "li") do
                                href = string(k, "_", generate_html_id(s))
                                print(io, "", s, "")
                            end
                        end
                    end
                end
            end
        end
        writehtml(io,ents, config)
    end
    return ents
end

function writehtml(io::IO, ents::EntriesHtml, config::Config)
    wrap(io, "div", "class='entries'") do
        for (modname, obj, ent) in ents.entries
            writehtml(io, modname, obj, ent, config)
        end
    end
end

function writehtml{category}(io::IO, modname, obj, ent::Entry{category}, config::Config)
    wrap(io, "div", "class='entry'") do
        objname = writeobj(obj, ent)
        idname = "$(category)_" * generate_html_id(objname)
        wrap(io, "div", "id='$(idname)' class='entry-name category-$(category)'") do
            print(io, "
[$(category)] —
 ")
            println(io, objname)
        end
        wrap(io, "div", "class='entry-body'") do
            writehtml(io, docs(ent))
            if has_output_metadata(ent, config)
                wrap(io, "div", "class='entry-meta'") do
                    println(io, "Details:")
                    wrap(io, "table", "class='meta-table'") do
                        for m in config.metadata_order
                            if haskey(ent.data, m)
                                wrap(io, "tr") do
                                    print(io, "", m, ":	")
                                    wrap(io, "td") do
                                        writehtml(io, Meta{m}(ent.data[m]))
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

function writehtml(io::IO, md::Meta)
    println(io, "", md.content, "")
end

function writemime(io::IO, mime::MIME"text/html", m::Meta{:parameters})
    wrap(io, "table") do
        for (k, v) in m.content
            wrap(io, "tr") do
                wrap(io, "td") do
                    println(io, "", k, "")
                end
                wrap(io, "td") do
                    writemime(io, MIME("text/html"), Markdown.parse(v))
                end
            end
        end
    end
end

function writehtml(io::IO, m::Meta{:source})
    path = last(split(m.content[2], r"v[\d\.]+(/|\\)"))
    print(io, "$(path):$(m.content[1])")
end

function headerhtml(io::IO, doc::Metadata, config::Config)
    println(io, """
    <!doctype html>

    <meta charset="utf-8">

    <title>$(doc.modname)</title>

    <link rel="stylesheet" type="text/css"
          href="https://cdnjs.cloudflare.com/ajax/libs/codemirror/4.5.0/codemirror.min.css">

    <link rel="stylesheet" type="text/css" href="static/custom.css">

$(doc.modname)

    """)
end

function footerhtml(io::IO, doc::Metadata, config::Config)
    println(io, """

    $(config.mathjax ? "" : "")

    """)
end

function wrap(fn::Function, io::IO, tag::AbstractString, attributes::AbstractString = "")
    println(io, "<", tag, " ", attributes, ">")
    fn()
    println(io, "</", tag, ">")
end

## Docs-specific rendering ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

function writehtml(io::IO, docs::Docile.Interface.Docs{:md})
    writemime(io, MIME("text/html"), parsed(docs))
end

### API-Index ----------------------------------------------------------------------------

function save(file::AbstractString, mime::MIME"text/html", index::Index, c::Config)
    throw(ArgumentError("The html format does currently not support saving of separate API-Index pages.)"))
end


text:
#=
4 target center out

Graphics Library(ies):
Cairo

States:
1) Start
2) In center
3) Target appears
4) Enter target
5) Hold time achieved (success)
6) Reward
7) End
8) Fail

Target locations:
0, 90, 180, 270 degrees

Target and cursor are circles

=#

export Task_COT

#=
Data structure for 4 target center out
=#
mutable struct Task_COT <: Task
    win_t::Gtk.GtkWindowLeaf#task Window handle
    c1::Gtk.GtkCanvasLeaf#task canvas handle
    win_c::Gtk.GtkWindowLeaf#control window handle
    tar::NTuple{2,Float64} #Target location
    success::Int64 #total successes
    total::Int64 #total trials
    state::Int64 #current state
    hold_time::Float64 #hold time
    juice_amount::Float64
    cursor_size::Float64
    target_size::Float64
    reach_radius::Float64
    it_time::Float64 #iter-trial interval
end

#=
Constructors
=#

function Task_COT()

    #Draw Task Screen
    c1=Canvas(800,800)
    @guarded draw(c1) do widget
    ctx = getgc(c1)
    set_source_rgb(ctx,0,0,0)
    paint(ctx)
    end
    show(c1)
    grid1 = Grid()
    push!(grid1,c1)
    win1 = Window(grid1,"Center Out Task")
    showall(win1)

    #Draw Control widget
    grid2 = Grid()
    button_start=Button("Start Task")
    grid2[1,1]=button_start
    win2 = Window(grid2,"Center Out Control")
    showall(win2)

    hold_time=300.0
    juice_amount=1.0
    cursor_size=5.0
    target_size=5.0
    reach_radius=5.0
    it_time=1.0

    Task_COT(win1,c1,win2,(1.0,0.0),0,0,0,hold_time,juice_amount,cursor_size,target_size,reach_radius,it_time)
end

#=
This will create any additional GUIs necessary for experimental control

=#
function init_task(myt::Task_COT,rhd::RHD2000)

end

function do_task(myt::Task_COT,rhd::RHD2000)

    #Task logic

    #Check if total time has exceeded threshold

    if myt.state==1

    elseif myt.state==2
        #Determine if in center
        #If time > x, then progress to state=3, make target appear
        #if time < x, then continue waiting
    elseif myt.state==3
        #check if in target, if yes then progress to state=4
    elseif myt.state==4
        #determine if in target
        #if time > x, then progress to state 5
        #if time < x, then continue waiting

        #if outside target, go back to state 3
    elseif myt.state==5
        #If time > x, deliver reward, go to 6
    elseif myt.state==6
        #if time > it_time go to 7
    elseif myt.state==7

    else

    end
    nothing
end

function save_task(myt::Task_COT,rhd::RHD2000)

end

#=
GUI callbacks
=#


text:
mutable struct ConstraintParams{T}
	ϕ::T  	    # penalty scaling parameter
	μ0::T 	    # initial penalty parameter
	μ_max::T    # max penalty parameter
	λ_max::T    # max Lagrange multiplier
end

function ConstraintParams(ϕ::T1 = 10, μ0::T2 = 1.0, μ_max::T3 = 1e8, λ_max::T4 = 1e8) where {T1,T2,T3,T4}
	T = promote_type(T1,T2,T3,T4)
	ConstraintParams(T(ϕ), T(μ0), T(μ_max), T(λ_max))
end

# Iteration
Base.iterate(conSet::AbstractConstraintSet) =
	isempty(get_convals(conSet)) ? nothing : (get_convals(conSet)[1].con,1)
Base.iterate(conSet::AbstractConstraintSet, state::Int) =
	state >= length(conSet) ? nothing : (get_convals(conSet)[state+1].con, state+1)
@inline Base.length(conSet) = length(get_convals(conSet))
Base.IteratorSize(::AbstractConstraintSet) = Base.HasLength()
Base.IteratorEltype(::AbstractConstraintSet) = Base.HasEltype()
Base.eltype(::AbstractConstraintSet) = AbstractConstraint

# Constraint Evaluation
function evaluate!(conSet::AbstractConstraintSet, Z::AbstractTrajectory)
    for conval in get_convals(conSet)
        evaluate!(conval, Z)
    end
end

function jacobian!(conSet::AbstractConstraintSet, Z::AbstractTrajectory)
    for conval in get_convals(conSet)
        jacobian!(conval, Z)
    end
end

function ∇jacobian!(G::Vector{<:Matrix}, conSet::AbstractConstraintSet, Z::AbstractTrajectory,
		λ::Vector{<:Vector})
	for (i,conval) in enumerate(get_convals(conSet))
		∇jacobian!(G[i], conval, Z, λ[i])
	end
end

function error_expansion!(conSet::AbstractConstraintSet, model::AbstractModel, G)
	@assert get_convals(conSet) == get_errvals(conSet)
	return nothing
end

function error_expansion!(conSet::AbstractConstraintSet, model::LieGroupModel, G)
	convals = get_convals(conSet)
	errvals = get_errvals(conSet)
	for i in eachindex(errvals)
		error_expansion!(errvals[i], convals[i], model, G)
	end
end

# Max values
function max_violation(conSet::AbstractConstraintSet)
	max_violation!(conSet)
	if isempty(conSet)
		return 0
	else
		return maximum(conSet.c_max)
	end
end

function max_violation!(conSet::AbstractConstraintSet)
	convals = get_convals(conSet)
	T = eltype(conSet.c_max)
    for i in eachindex(convals)
        max_violation!(convals[i])
        conSet.c_max[i] = maximum(convals[i].c_max::Vector{T})
    end
	return nothing
end

"""
	norm_violation(conSet, [p=2])
	norm_violation(conVal, [p=2])

Calculate the norm of the violation of a `AbstractConstraintSet`, a `Conval` using the `p`-norm.

	norm_violation(sense::ConstraintSense, v::AbstractVector, [p=2])

Calculate the `p`-norm of constraint violations given the vector of constraint values `v`,
	for a constraint of sense `sense` (either `Inequality()` or `Equality()`).
	Assumes that positive values are violations for inequality constraints.
"""
function norm_violation(conSet::AbstractConstraintSet, p=2)
	norm_violation!(conSet, p)
	norm(conSet.c_max, p)
end

function norm_violation!(conSet::AbstractConstraintSet, p=2)
	convals = get_convals(conSet)
	T = eltype(conSet.c_max)
	for i in eachindex(convals)
		norm_violation!(convals[i], p)
		c_max = convals[i].c_max::Vector{T}
		conSet.c_max[i] = norm(c_max, p)
	end
end

function norm_dgrad(conSet::AbstractConstraintSet, dx::AbstractTrajectory, p=1)
	convals = get_convals(conSet)
	T = eltype(conSet.c_max)
	for i in eachindex(convals)
		norm_dgrad!(convals[i], dx, p)
		c_max = convals[i].c_max::Vector{T}
		conSet.c_max[i] = sum(c_max)
	end
	return sum(conSet.c_max)
end

"""
	findmax_violation(conSet)

Return details on the where the largest violation occurs. Returns a string giving the
constraint type, time step index, and index into the constraint.
"""
function findmax_violation(conSet::AbstractConstraintSet)
	max_violation!(conSet)
	c_max0, j_con = findmax(conSet.c_max) # which constraint
	if c_max0 < eps()
		return "No constraints violated"
	end
	convals = get_convals(conSet)
	conval = convals[j_con]
	i_con = findmax(conval.c_max)[2]  # which index
	k_con = conval.inds[i_con] # time step
	con_sense = sense(conval.con)
	viol = [violation(con_sense, v) for v in conval.vals[i_con]]
	c_max, i_max = findmax(viol)  # index into constraint
	@assert c_max == c_max0
	con_name = string(typeof(conval.con).name)
	return con_name * " at time step $k_con at " * con_label(conval.con, i_max)
end


text:
using Clang

# Customize how functions, constants, and structs are written
# Copied and modified from CUDArt.
const skip_expr = []
const skip_error_check = []
function rewriter(ex::Expr)
    if in(ex, skip_expr)
        return :()
    end
    # Empty types get converted to Void
    if ex.head == :type
        a3 = ex.args[3]
        if isempty(a3.args)
            objname = ex.args[2]
            return :(typealias $objname Void)
        end
    end
    ex.head == :function || return ex
    decl, body = ex.args[1], ex.args[2]
    # omit types from function prototypes
    for i = 2:length(decl.args)
        a = decl.args[i]
        if a.head == :(::)
            decl.args[i] = a.args[1]
        end
    end
    # Error-check functions that return a cudaError_t (with some omissions)
    ccallexpr = body.args[1]
    if ccallexpr.head != :ccall
        error("Unexpected body expression: ", body)
    end
    rettype = ccallexpr.args[2]
    if rettype == :cusparseStatus_t
        fname = decl.args[1]
        if !in(fname, skip_error_check)
            body.args[1] = Expr(:call, :cusparseCheck, deepcopy(ccallexpr))
        end
    end
    ex
end

rewriter(A::Array) = [rewriter(a) for a in A]

rewriter(s::Symbol) = string(s)

rewriter(arg) = arg

run(wrap_c.init(headers = ["/share/apps/cuda/cuda-6.5/include/cusparse.h"],
                # index = None,
                common_file="types.jl",
                output_file="libcusparse.jl",
                # output_dir = "",
                # clang_args = ASCIIString[],
                # clang_includes = ["/usr/lib/clang/3.4.2/include","/usr/local/cuda/include","/usr/include"],
                clang_includes = ["/usr/local/cuda-6.5/include"],
                # clang_diagnostics = true,
                header_wrapped=(header,cursorname)->(contains(cursorname,"cusparse")),
                # header_wrapped=(header,cursorname)->true,
                header_library=x->"libcusparse",
                # header_outputfile = None,
                # cursor_wrapped=(cursorname,cursor)->!isempty(cursorname),
                cursor_wrapped=(cursorname,cursor)->(!isempty(cursorname) && (cursorname!="CUSPARSEAPI")),
                # options = InternalOptions(),
                rewriter = rewriter,
                ))


text:
for version in (v"1.0.0", v"1.1.0")
    v_description = VERSION_DESCRIPTION[version]
    hash = VERSION_HASH[version]

    dataset2description = VERSION_DATASET_DESCRIPTION[version]
    dataset2checksum = VERSION_DATASET_CHECKSUM[version]

    for dataset in ("pudl-eia860-eia923", "pudl-eia860-eia923-epacems", "pudl-ferc1")
        name = dataset * "-v$version"
        description = v_description * "\n" * dataset2description[dataset]
        checksum = get(dataset2checksum, dataset, nothing)

        register(
            DataDep(
                name,
                description,
                "https://zenodo.org/api/files/$(hash)/$(dataset).tgz",
                checksum,
                post_fetch_method=unpack,
            )
        )
    end
end


text:
# AbstractArray interface for JSON.Array

Base.IteratorSize(::Type{JSON.Array{W, T}}) where W where T = Base.SizeUnknown()
Base.IteratorEltype(::Type{JSON.Array{W, T}}) where W where T = Base.EltypeUnknown()

# Access

Base.getindex(a::JSON.Array{W}, i::Integer) where W =
    getvalue(W, a.s, getindex_ic(a, i)...)

# Dimensions

function Base.length(a::JSON.Array)
    @lazywarn "length(::JSON.Array) is slow!" stacktrace()
    collection_length(a)
end

Base.size(a::JSON.Array) = (length(a), )

# Iterate

Base.iterate(j::JSON.Array, i = (j.i, 0x00)) = _iterate(j, i)
function _iterate(j::JSON.Array{W}, (i, c)) where W
    i, c = nextindex(j, i, c)
    if c == ']'
        return nothing
    end
    return getvalue(W, j.s, i, c), (i, c)
end

# IOString Wrappers

Base.length(j::JSON.Array{W, IOString{T}}) where W where T =
    pump(() -> collection_length(j), j.s)

Base.getindex(j::JSON.Array{W, IOString{T}}, i::Integer) where W where T =
    pump(() -> getvalue(W, j.s, getindex_ic(j, i)...), j.s)

Base.iterate(j::JSON.Array{W, IOString{T}}, i = (j.i, 0x00)) where W where T =
    pump(() -> _iterate(j, i), j.s)


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

const SMARTS_BOND_SYMBOL = Dict(
    '-' => QueryFormula(:bondorder, 1),
    '=' => QueryFormula(:bondorder, 2),
    '#' => QueryFormula(:bondorder, 3),
    '@' => QueryFormula(:isringbond, true),
    ':' => QueryFormula(:isaromaticbond, true),
    '/' => QueryFormula(:stereo, :up),
    '\\' => QueryFormula(:stereo, :down)
)

defaultbond(state::SmilesParser) = SmilesBond()
defaultbond(state::SmartsParser) = SmartsBond(QueryFormula(:defaultbond, true))

"""
    bondsymbol!(state::SmartsParserState) -> Union{Pair,Nothing}

BondSymbol <- [-=#@:/\\] / '/?' / '\\?'
"""
function bondsymbol!(state::SmartsParserState)
    sym1 = read(state)
    sym2 = lookahead(state, 1)
    if sym1 == '/' && sym2 == '?'
        forward!(state, 2)
        return QueryFormula(:not, QueryFormula(:stereo, :down))
    elseif sym1 == '\\' && sym2 == '?'
        forward!(state, 2)
        return QueryFormula(:not, QueryFormula(:stereo, :up))
    elseif sym1 in keys(SMARTS_BOND_SYMBOL)
        forward!(state)
        return SMARTS_BOND_SYMBOL[sym1]
    end
    # Implicit single bond returns nothing
end

"""
    bond!(state::SmilesParser) -> Union{SmilesBond,Nothing}

Bond <- BondSymbol?
"""
function bond!(state::SmilesParser)
    fml = bondsymbol!(state)
    if fml === nothing
        return
    elseif fml.key == :bondorder
        return SmilesBond(fml.value)
    elseif fml.key == :isaromaticbond
        return SmilesBond(1, true, :unspecified)
    elseif fml.key == :stereo
        return SmilesBond(1, false, fml.value)
    end
    # return nothing
end

"""
    bond!(state::SmartsParser) -> Union{SmartsBond,Nothing}

Bond <- '~' / (BondSymbol / LogicalCond)?
"""
function bond!(state::SmartsParser)
    if read(state) == '~'
        forward!(state)
        return SmartsBond()
    end
    fml = lglowand!(state, bondsymbol!)
    if fml !== nothing
        fml = tidyformula(fml)
        return SmartsBond(fml)
    end
    # return nothing: Invalid bond token or implicit single bond
end


text:
# https://adventofcode.com/2021/day/5
using AdventOfCode
input = parse.("(".*replace.(readlines("2021/data/day_5.txt"), "->"=>")=>(").*")")
rangestep(x,y,step)=range(x,y,step=step)
function part_1(input)
    lines=filter(x->any(first(x).==last(x)),input)
    limit=maximum(maximum.(vcat(first.(lines),last.(lines))))
    field=zeros(Int,limit,limit)
    ranges=(x->range.(first(x),last(x),step=sign(sum(last(x).-first(x))))).(lines)
    for r in ranges
        field[r[1],r[2]].+=1
    end
    sum(>(1),field)
end

function part_2(input)
    lines=input
    limit=maximum(maximum.(vcat(first.(lines),last.(lines))))
    field=zeros(Int,limit,limit)
    ranges=(x->rangestep.(first(x),last(x),sign.(last(x).-first(x).+eps()).|>Int)).(lines)
    for r in ranges
        field[CartesianIndex.(r[1],r[2])].+=1
    end
    sum(>(1),field)
end


text:
include("Input.jl")
include("Output.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
"""
    validate(instance, solution)

Check if the given solution fits the instance

"""
function validate(instance::Instance, L, T, J)
end

# Checking the feasibility of a solution equates Checking
# - presence of all jobs, with starting times after the earliest arrival times
# - respect for precedence constraints
# - no intersections (nor invali)

function check_jobdata(solution)
    J = jobs(solution)
    PG = precedence_graph(J)
    ts = zeros(Int, length(J))
    td = zeros(Int, length(J))
    ispresent = zeros(Int, length(J))
    for T in trajectories(solution)
        S = starttimes(T)
        for (j, t_s) in pairs(S)
            @test t_arrival(J[j]) <= t_s
            ts[j] = t_s
            ispresent[j] += 1
        end
        D = durations(T)
        for (j, t_d) in pairs(D)
            @test t_processing(J[j]) == t_d
            td[j] = t_d
        end
    end
    tc = ts .+ td
    for (s,d) in edges(PG)
        tc[s] < ts[d]
    end

end


text:
# I want to try a pure MC simulation that doesn't rely on a neural network and
# see how things
using Statistics: mean, quantile
using StatsBase: countmap
using Game2048
using Game2048: bitboard_to_array, DIRS, transform_4corners_to_4cols
using Serialization: serialize, deserialize
using Base.Threads
using BSON
using BenchmarkTools
using Flux.Data: DataLoader
using Flux: throttle
using Flux: softmax, onecold
# make a nn to study these games
using Flux, CUDA
using Flux: onehot, logitbinarycrossentropy, onehotbatch
using Chain
CUDA.allowscalar(false)

model_group_cnn = Flux.Chain(
    Conv((4,1), 1=>128, relu),
    x->reshape(x, 2, 2, 128, :),
    Conv((2,2), 128=>64),
    Flux.flatten,
    SkipConnection(Dense(64, 64, relu), (mx,x) -> mx.+x),
    SkipConnection(Dense(64, 64), (mx,x) -> mx.+x),
    Dense(64, 4)
) |> gpu

model_group_cnn_cpu = model_group_cnn |> cpu

function group_cnn_policy(bitboard)
    @chain bitboard begin
        bitboard_to_array
        transform_4corners_to_4cols
        reshape(_, 4, 4, 1, :) ./ Float32(11)
        model_group_cnn_cpu
        softmax
        reshape(4)
    end
end

function get_move_probs(one_game::GameRecorder)
    # which move did I regret the most?
    game_matrix = mapreduce((a,b)->cat(a,b, dims=4), one_game.states) do state
        @chain state begin
            bitboard_to_array
            _ ./ 11
            transform_4corners_to_4cols
            reshape(4, 4, 1)
        end
    end |> cu

    softmax(model_group_cnn(game_matrix)) |> cpu
end

if false && isfile("group_pcpu")
    pcpu = deserialize("group_pcpu")
    Flux.loadparams!(model_group_cpu, pcpu)
end

# includet("src/neural-network-utils.jl")
# includet("src/policy-player.jl")

@time y = deserialize("y-data.jls");

if false
    @time x = deserialize("x-data.jls");
    @time x = mapslices(transform_4corners_to_4cols, x, dims=[1,2,3]); # about 1min
    @time serialize("x-data-corner-transformed.jls", x);
end

@time x = deserialize("x-data-corner-transformed.jls");

sample_size = round(Int, size(x, 4)*0.99)

train_x = x[:, :, :, 1:sample_size] |> cu;
test_x = x[:, :, :, sample_size+1:end] |> cu;

train_y = y[:, 1:sample_size] |> cu;
test_y = y[:, sample_size+1:end] |> cu;

x=nothing; y=nothing;GC.gc() # free up memory

loss(x, y) = logitbinarycrossentropy(model_group_cnn(x), y)

@time softmax(model_group_cnn(test_x))
@time softmax(model_group_cnn(test_x))
@time loss(test_x, test_y)
@time loss(test_x, test_y)

opt=ADAM()

dl = DataLoader((train_x, train_y); batchsize=2^5, shuffle=true);

accuracy(x, y) = mean(onecold(softmax(model_group_cnn(x))) .== onecold(y))

# @time accuracy(train_x, train_y)

last_loss = loss(test_x, test_y)
last_acc = accuracy(test_x, test_y)
worse_cnt = 0

function evalcb()
    global last_loss, worse_cnt, last_acc
    new_loss = loss(test_x, test_y)
    new_acc = accuracy(test_x, test_y)
    # println("train loss: $(loss(train_x, train_y)) train accuracy: $(accuracy(train_x, train_y))")
    println("test loss: $(new_loss) test accuracy: $(new_acc)")

    if (new_loss < last_loss) | (new_acc > last_acc)
        last_loss = new_loss
        last_acc = new_acc
        worse_cnt = 0
    else
        worse_cnt += 1
        if worse_cnt == 8
            Flux.stop()
        end
    end
end

@time evalcb()

cb = throttle(evalcb, 88)

dl2 = DataLoader((train_x[:,:,:,1:100], train_y[:, 1:100]); batchsize=32, shuffle=true);

@time Flux.train!(loss, Flux.params(model_group_cnn), dl2, opt; cb); # compile
@time Flux.train!(loss, Flux.params(model_group_cnn), dl, opt; cb); # compile
@time Flux.@epochs 8 Flux.train!(loss, Flux.params(model_group_cnn), dl, opt; cb);
@time Flux.@epochs 88 Flux.train!(loss, Flux.params(model_group_cnn), dl, opt; cb);
@time Flux.@epochs 888 Flux.train!(loss, Flux.params(model_group_cnn), dl, opt; cb);
# @time Flux.@epochs 8888 Flux.train!(loss, Flux.params(model_group_cnn), dl, opt; cb);

### let's play the game vs random policy without search and see if it's better than random already

# function compareonce()
#     bboard = initbboard()
#     m1 = maximum(Game2048.play_game_with_policy(randompolicy, bboard))
#     m2 =maximum(Game2048.play_game_with_policy(group_cnn_policy, bboard))
#     m1, m2
# end

# yy = [compareonce() for _ in 1:1000]

# using StatsBase
# mean([m1-m2 for (m1, m2) in yy])

# countmap([m1 for (m1, m2) in yy])
# countmap([m2 for (m1, m2) in yy])

function go_back(game1)
    i = 1
    m1 = maximum(Game2048.play_via_monte_carlo_wo_recorder(randompolicy; n=100, bitboard = game1.states[end-i]))
    while m1 < 11
        i += 1
        m1 = maximum(Game2048.play_via_monte_carlo_wo_recorder(randompolicy; n=100, bitboard = game1.states[end-i]))
    end
    return i
end

function make_new_training_data(bitboard)
    @chain bitboard begin
        bitboard_to_array
        transform_4corners_to_4cols
        reshape(_, 4, 4, 1, :) ./ Float32(11)
    end
end

function train_n_times(games)
    model_group_cnn_cpu = model_group_cnn |> cpu

    function group_cnn_policy(bitboard)
        @chain bitboard begin
            bitboard_to_array
            transform_4corners_to_4cols
            reshape(_, 4, 4, 1, :) ./ Float32(11)
            model_group_cnn_cpu
            softmax
            reshape(4)
        end
    end
    ## which one is it predicting better
    game1 = Game2048.play_game_with_policy_w_record(group_cnn_policy, initbboard());

    # For each state compute their policy
    policies = reduce((x,y)->cat(x,y, dims=2), group_cnn_policy.(game1.states))

    # progress move backwards until we can't reach 11
    goback_i = go_back(game1)

    # assess each state by finding the best move
    best_moves_according_to_random = Game2048.find_best_move.(@view(game1.states[1:end-goback_i]), 100, randompolicy)
    # countmap(best_moves_according_to_random)

    push!(games, (game1, best_moves_according_to_random))

    lossval = logitbinarycrossentropy(policies[:,1:end-goback_i], onehotbatch(Int.(best_moves_according_to_random) .+ 1, 1:4))
    acc = mean(onecold(policies[:,1:end-goback_i]) .== Int.(best_moves_according_to_random) .+ 1)
    println("accuracy: $(acc); loss $(lossval); best: $(maximum(game1.final_state))")

    # find the move that it's the surest about
    most_sure_move = argmax(policies)[2]
    println("")
    display(game1.states[most_sure_move])
    println("move made: $(game1.moves[most_sure_move])")
    println("policy $(policies[:, most_sure_move])")
    println("maxpolicy $(maximum(policies[:, most_sure_move]))")
    println("best policy move: $(DIRS[argmax(policies[:, most_sure_move])])")
    println("best move: $(best_moves_according_to_random[most_sure_move])")

    # trim out that many moves and make it learn bitches
    new_train_x = reduce((x,y)->cat(x, y, dims=4), make_new_training_data.(@view(game1.states[1:end-goback_i]))) |> cu
    new_train_y = onehotbatch(Int.(best_moves_according_to_random) .+ 1, 1:4) |> cu

    Flux.train!(loss, Flux.params(model_group_cnn), [(new_train_x, new_train_y)], opt);
end

games = []
@time train_n_times(games)

while true
    train_n_times(games)
end

model_group_cnn_cpu = model_group_cnn |> cpu
pcpu = Flux.params(model_group_cnn_cpu)
serialize("group_pcpu", pcpu)
pcpu = deserialize("group_pcpu")
Flux.loadparams!(model_group_cnn_cpu, pcpu)

@time x = [maximum(Game2048.play_via_monte_carlo_wo_recorder(randompolicy; n=10)) for _ in 1:10]
@time y = [maximum(Game2048.play_via_monte_carlo_wo_recorder(group_cnn_policy; n=30)) for _ in 1:10]

@time xx1 = [Game2048.play_via_monte_carlo_wo_recorder(randompolicy; n=1) for _ in 1:10]
@time xx2 = [Game2048.play_via_monte_carlo_wo_recorder(group_cnn_policy; n=1) for _ in 1:10]

maximum.(xx1)
maximum.(xx2)

b = initbboard()

group_cnn_policy(b)
group_cnn_policy(rotl90(b))
group_cnn_policy(rotl90(b) |> rotl90)
group_cnn_policy(rotl90(b) |> rotl90 |> rotl90)

function rotate4(board)
    reduce((x,y)->cat(x, y, dims=3), (board, rotl90(board), rot180(board), rotr90(board)))
end

board = bitboard_to_array(initbboard())
rotate4(board)

function normalize4(board)
    reshape(board[1], 4, 4, 1, :) ./ Float32(11),
    reshape(board[2], 4, 4, 1, :) ./ Float32(11),
    reshape(board[3], 4, 4, 1, :) ./ Float32(11),
    reshape(board[4], 4, 4, 1, :) ./ Float32(11)
end

cat(board, rotr90(board), dims=3)

function rotated_cnn(bitboard)
    @chain bitboard begin
        bitboard_to_array
        rotate4
        mapslices(transform_4corners_to_4cols, _, dims=(1,2))
        reshape(_, 4, 4, 1, 4) ./ Float32(11)
        model_group_cnn_cpu
    end
end

b = initbboard()

function firstpart(bitboard)
    @chain bitboard begin
        bitboard_to_array
        rotate4
    end
end

function circular!(x)
    x[1], x[2], x[3], x[4] = x[2], x[3], x[4], x[1]
    x
end

rb = rotated_cnn(b)

rb[:, 2] .= circular!(rb[:, 2])
rb[:, 3] .= circular!(circular!(rb[:, 3]))
rb[:, 4] .= circular!(circular!(circular(rb[:, 4])))

using StatsBase
mapslices(rb, dims=2) do slice
    sum((slice .- mean(slice)).^2)
end |> sum

# @time BSON.@save "model_cpu" model_cpu

pcpu = Flux.params(model_group_cpu)
serialize("group_pcpu", pcpu)
pcpu = deserialize("group_pcpu")
Flux.loadparams!(model_group_cpu, pcpu)

# I want to see which move it's the most sure about

full_results = mapreduce(deserialize, vcat, readdir("c:/data/game2048/raw/", join=true)[2:2]);

one_game = full_results[1]

move_probs1 = move_probs(one_game)
# find the move that regret most
moves_idx = Int.(one_game.moves) .+ 1

regret_scores = [maximum(r)/r[idx] for (r, idx) in zip(eachslice(move_probs1, dims=2), moves_idx)]
pos = argmax(regret_scores)
state, actual_move pos =  one_game.states[pos], one_game.moves[pos] # the move that was made

pos, state, actual_move = regrets(one_game)

state

# best move according to random play
countmap([Game2048.find_best_move(one_game.states[pos], 1000, Game2048.randompolicy) for i in 1:100])

# this approach makes the
prob = @chain one_game.states[pos] begin
    bitboard_to_array
    transform_4corners_to_4cols
    reshape(_, 4, 4, 1, :) ./ 11
    cu
    model_group_cnn
    softmax
    cpu
    reshape(4)
end

DIRS[argmax(prob)] #right is best move

model_group_cnn_cpu = model_group_cnn |> cpu

@time group_cnn_policy(one_game.states[pos])

@time [play_game_with_policy(group_cnn_policy, one_game.states[pos]) |> maximum for i in 1:1000] |> countmap
@time [play_game_with_policy(randompolicy, one_game.states[pos]) |> maximum for i in 1:1000] |> countmap

@time meh = Game2048.play_via_monte_carlo_w_recorder(randompolicy; n=1000, bitboard=one_game.states[pos])
meh.final_state
@time meh2 = play_via_monte_carlo_w_recorder(group_cnn_policy; n=1000, bitboard=one_game.states[pos])

# the rest is just checking to see how the recommendations change in relations to recommendations
best_move = prob |> onecold |> cpu |>
    (x->DIRS[x[1]])

# rotate left; up is now the best move should be right
one_game.states[pos] |>
    bitboard_to_array |>
    rotl90 |>
    transform_4corners_to_4cols |>
    (x->reshape(x, 4, 4, 1, :) ./ 11) |>
    cu |> model_group_cnn |> softmax

# up should be the best move
one_game.states[pos] |>
    bitboard_to_array |>
    rot180 |>
    transform_4corners_to_4cols |>
    (x->reshape(x, 4, 4, 1, :) ./ 11) |>
    cu |> model_group_cnn |> softmax

one_game.states[pos] |>
    bitboard_to_array |>
    rotr90 |>
    transform_4corners_to_4cols |>
    (x->reshape(x, 4, 4, 1, :) ./ 11) |>
    cu |> model_group_cnn |> softmax

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

function makey(state1)
    dir_cnt = zeros(Int, 4)
    for i in [Int(Game2048.find_best_move(state1, 1000, Game2048.randompolicy)) + 1 for i in 1:10]
        dir_cnt[i] += 1
    end
    dir_cnt./100
end

function normal_to_model(state1)
    @chain state1 begin
        bitboard_to_array
        transform_4corners_to_4cols
        reshape(_, 4, 4, 1, :) ./ 11
    end
end

state = gamerecord_nn.states[1]

model_group_cnn_cpu = model_group_cnn |> cpu

function model(x)
    # x = normal_to_model(state)

    xs = Ref(reshape(x,4,4)) .|> (identity, rotr90, rot180, rotl90) .|> (x->reshape(x, 4, 4, 1, :))

    result = model_group_cnn_cpu.(xs)

    for (i, r) in enumerate(result)
        r .= @view r[mod1.((1:4) .+ i, 4)]
    end

    # compute distance to mean
    rr = hcat(result...)
    map(eachrow(rr)) do row
        mean(row)
    end
end

loss(x, y) = Flux.logitbinarycrossentropy(model(x), y)

loss(x, [0.0, 1.0, 0.0, 0.0])

opt=ADAM()
@time Flux.train!(loss, Flux.params(model_group_cnn), [(x, [0.0, 1.0, 0.0, 0.0])], opt);

while true
    # play a game from scratch
    @time gamerecord_nn = play_via_monte_carlo_w_recorder(group_cnn_policy; n=1, bitboard=initbboard())

    move_probs = get_move_probs(gamerecord_nn)

    # look for the top 32 moves that need correction since computing every move is taking a long time
    # seems to be a better strategy to play 100 and then stop if maximum exceeds 1000
    # @time best_play = play_via_monte_carlo_w_recorder(randompolicy; n=1000, bitboard = gamerecord_nn.states[end-1])
    best_play100 = play_via_monte_carlo_w_recorder(randompolicy; n=100, bitboard = gamerecord_nn.states[end-1])
    while maximum(best_play100.final_state) < 11
        best_play100 = play_via_monte_carlo_w_recorder(randompolicy; n=100, bitboard = gamerecord_nn.states[end-1])
    end

    best_play100

    Flux.train!(loss, Flux.params(model_group_cnn), [(new_train_x, new_train_y)], opt);

    # use random play to figure out the best moves
    x = cat(normal_to_model.(gamerecord_nn.states)..., dims=4)
    wy = reduce(hcat, makey.(gamerecord_nn.states))

    obj = (gamerecorder = gamerecord_nn, x=x, y=y)
    hashofit = hash(obj)

    using Serialization: serialize
    # serialize(raw"C:\data\game2048\self-play-data\$hashofit", obj)

    @time Flux.train!(loss, Flux.params(model_group_cnn), [(cu(x), cu(y))], opt; cb); # compile

    model_group_cpu=model_group_cnn |> cpu
    pcpu = Flux.params(model_group_cpu)
    serialize("group_pcpu", pcpu)
    # pcpu = deserialize("group_pcpu")
    # Flux.loadparams!(model_group_cpu, pcpu)
    # model_group_cnn = model_group_cpu |> gpu
end


text:
@testset "compiler.jl" begin
    @testset "assume" begin
        @model function test_assume()
            x ~ Bernoulli(1)
            y ~ Bernoulli(x / 2)
            return x, y
        end

        smc = SMC()
        pg = PG(10)

        res1 = sample(test_assume(), smc, 1000)
        res2 = sample(test_assume(), pg, 1000)

        check_numerical(res1, [:y], [0.5]; atol=0.1)
        check_numerical(res2, [:y], [0.5]; atol=0.1)

        # Check that all xs are 1.
        @test all(isone, res1[:x])
        @test all(isone, res2[:x])
    end
    @testset "beta binomial" begin
        prior = Beta(2, 2)
        obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]
        exact = Beta(prior.α + sum(obs), prior.β + length(obs) - sum(obs))
        meanp = exact.α / (exact.α + exact.β)

        @model function testbb(obs)
            p ~ Beta(2, 2)
            x ~ Bernoulli(p)
            for i in 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            return p, x
        end

        smc = SMC()
        pg = PG(10)
        gibbs = Gibbs(HMC(0.2, 3, :p), PG(10, :x))

        chn_s = sample(testbb(obs), smc, 1000)
        chn_p = sample(testbb(obs), pg, 2000)
        chn_g = sample(testbb(obs), gibbs, 1500)

        check_numerical(chn_s, [:p], [meanp]; atol=0.05)
        check_numerical(chn_p, [:x], [meanp]; atol=0.1)
        check_numerical(chn_g, [:x], [meanp]; atol=0.1)
    end
    @testset "forbid global" begin
        xs = [1.5 2.0]
        # xx = 1

        @model function fggibbstest(xs)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            # xx ~ Normal(m, sqrt(s)) # this is illegal

            for i in 1:length(xs)
                xs[i] ~ Normal(m, sqrt(s))
                # for xx in xs
                # xx ~ Normal(m, sqrt(s))
            end
            return s, m
        end

        gibbs = Gibbs(PG(10, :s), HMC(0.4, 8, :m))
        chain = sample(fggibbstest(xs), gibbs, 2)
    end
    @testset "new grammar" begin
        x = Float64[1 2]

        @model function gauss(x)
            priors = TArray{Float64}(2)
            priors[1] ~ InverseGamma(2, 3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            return priors
        end

        chain = sample(gauss(x), PG(10), 10)
        chain = sample(gauss(x), SMC(), 10)

        @model function gauss2(::Type{TV}=Vector{Float64}; x) where {TV}
            priors = TV(undef, 2)
            priors[1] ~ InverseGamma(2, 3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            return priors
        end

        chain = sample(gauss2(; x=x), PG(10), 10)
        chain = sample(gauss2(; x=x), SMC(), 10)

        chain = sample(gauss2(Vector{Float64}; x=x), PG(10), 10)
        chain = sample(gauss2(Vector{Float64}; x=x), SMC(), 10)
    end
    @testset "new interface" begin
        obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

        @model function newinterface(obs)
            p ~ Beta(2, 2)
            for i in 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            return p
        end

        chain = sample(
            newinterface(obs), HMC{Turing.ForwardDiffAD{2}}(0.75, 3, :p, :x), 100
        )
    end
    @testset "no return" begin
        @model function noreturn(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            for i in 1:length(x)
                x[i] ~ Normal(m, sqrt(s))
            end
        end

        chain = sample(noreturn([1.5 2.0]), HMC(0.15, 6), 1000)
        check_numerical(chain, [:s, :m], [49 / 24, 7 / 6])
    end
    @testset "observe" begin
        @model function test()
            z ~ Normal(0, 1)
            x ~ Bernoulli(1)
            1 ~ Bernoulli(x / 2)
            0 ~ Bernoulli(x / 2)
            return x
        end

        is = IS()
        smc = SMC()
        pg = PG(10)

        res_is = sample(test(), is, 10000)
        res_smc = sample(test(), smc, 1000)
        res_pg = sample(test(), pg, 100)

        @test all(isone, res_is[:x])
        @test res_is.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_smc[:x])
        @test res_smc.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_pg[:x])
    end
    @testset "sample" begin
        alg = Gibbs(HMC(0.2, 3, :m), PG(10, :s))
        chn = sample(gdemo_default, alg, 1000)
    end
    @testset "vectorization @." begin
        @model function vdemo1(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            @. x ~ Normal(m, sqrt(s))
            return s, m
        end

        alg = HMC(0.01, 5)
        x = randn(100)
        res = sample(vdemo1(x), alg, 250)

        @model function vdemo1b(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            @. x ~ Normal(m, $(sqrt(s)))
            return s, m
        end

        res = sample(vdemo1b(x), alg, 250)

        D = 2
        @model function vdemo2(x)
            μ ~ MvNormal(zeros(D), ones(D))
            @. x ~ $(MvNormal(μ, ones(D)))
        end

        alg = HMC(0.01, 5)
        res = sample(vdemo2(randn(D, 100)), alg, 250)

        # Vector assumptions
        N = 10
        setchunksize(N)
        alg = HMC(0.2, 4)

        @model function vdemo3()
            x = Vector{Real}(undef, N)
            for i in 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo3(), alg, 1000)

        # Test for vectorize UnivariateDistribution
        @model function vdemo4()
            x = Vector{Real}(undef, N)
            @. x ~ Normal(0, 2)
        end

        t_vec = @elapsed res = sample(vdemo4(), alg, 1000)

        @model vdemo5() = x ~ MvNormal(zeros(N), 2 * ones(N))

        t_mv = @elapsed res = sample(vdemo5(), alg, 1000)

        println("Time for")
        println("  Loop : ", t_loop)
        println("  Vec  : ", t_vec)
        println("  Mv   : ", t_mv)

        # Transformed test
        @model function vdemo6()
            x = Vector{Real}(undef, N)
            @. x ~ InverseGamma(2, 3)
        end

        sample(vdemo6(), alg, 1000)

        N = 3
        @model function vdemo7()
            x = Array{Real}(undef, N, N)
            @. x ~ [InverseGamma(2, 3) for i in 1:N]
        end

        sample(vdemo7(), alg, 1000)
    end
    @testset "vectorization .~" begin
        @model function vdemo1(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            x .~ Normal(m, sqrt(s))
            return s, m
        end

        alg = HMC(0.01, 5)
        x = randn(100)
        res = sample(vdemo1(x), alg, 250)

        D = 2
        @model function vdemo2(x)
            μ ~ MvNormal(zeros(D), ones(D))
            return x .~ MvNormal(μ, ones(D))
        end

        alg = HMC(0.01, 5)
        res = sample(vdemo2(randn(D, 100)), alg, 250)

        # Vector assumptions
        N = 10
        setchunksize(N)
        alg = HMC(0.2, 4)

        @model function vdemo3()
            x = Vector{Real}(undef, N)
            for i in 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo3(), alg, 1000)

        # Test for vectorize UnivariateDistribution
        @model function vdemo4()
            x = Vector{Real}(undef, N)
            return x .~ Normal(0, 2)
        end

        t_vec = @elapsed res = sample(vdemo4(), alg, 1000)

        @model vdemo5() = x ~ MvNormal(zeros(N), 2 * ones(N))

        t_mv = @elapsed res = sample(vdemo5(), alg, 1000)

        println("Time for")
        println("  Loop : ", t_loop)
        println("  Vec  : ", t_vec)
        println("  Mv   : ", t_mv)

        # Transformed test
        @model function vdemo6()
            x = Vector{Real}(undef, N)
            return x .~ InverseGamma(2, 3)
        end

        sample(vdemo6(), alg, 1000)

        @model function vdemo7()
            x = Array{Real}(undef, N, N)
            return x .~ [InverseGamma(2, 3) for i in 1:N]
        end

        sample(vdemo7(), alg, 1000)
    end
    @testset "Type parameters" begin
        N = 10
        setchunksize(N)
        alg = HMC(0.01, 5)
        x = randn(1000)
        @model function vdemo1(::Type{T}=Float64) where {T}
            x = Vector{T}(undef, N)
            for i in 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo1(), alg, 250)
        t_loop = @elapsed res = sample(vdemo1(Float64), alg, 250)

        vdemo1kw(; T) = vdemo1(T)
        t_loop = @elapsed res = sample(vdemo1kw(; T=Float64), alg, 250)

        @model function vdemo2(::Type{T}=Float64) where {T<:Real}
            x = Vector{T}(undef, N)
            @. x ~ Normal(0, 2)
        end

        t_vec = @elapsed res = sample(vdemo2(), alg, 250)
        t_vec = @elapsed res = sample(vdemo2(Float64), alg, 250)

        vdemo2kw(; T) = vdemo2(T)
        t_vec = @elapsed res = sample(vdemo2kw(; T=Float64), alg, 250)

        @model function vdemo3(::Type{TV}=Vector{Float64}) where {TV<:AbstractVector}
            x = TV(undef, N)
            @. x ~ InverseGamma(2, 3)
        end

        sample(vdemo3(), alg, 250)
        sample(vdemo3(Vector{Float64}), alg, 250)

        vdemo3kw(; T) = vdemo3(T)
        sample(vdemo3kw(; T=Vector{Float64}), alg, 250)
    end
end


text:
using Flux, DiffEqFlux, DiffEqFlux.GalacticOptim, Plots, LinearAlgebra, Optim
using KitBase, KitBase.FastGaussQuadrature, KitBase.ProgressMeter
import KitML

function flux_wall!(
    ff::T1,
    f::T2,
    u::T3,
    dt,
    rot = 1,
) where {
    T1<:AbstractVector{<:AbstractFloat},
    T2<:AbstractVector{<:AbstractFloat},
    T3<:AbstractVector{<:AbstractFloat},
}
    δ = heaviside.(u .* rot)
    fWall = 0.5 .* δ .+ f .* (1.0 .- δ)
    @. ff = u * fWall * dt

    return nothing
end

begin
    # setup
    set = Setup("radiation", "linesource", "1d1f1v", "kfvs", "bgk", 1, 2, "vanleer", "extra", 0.3, 0.8)

    # physical space
    x0 = 0
    x1 = 1
    nx = 100
    nxg = 0
    ps = PSpace1D(x0, x1, nx, nxg)

    # velocity space
    nu = 28
    points, weights = gausslegendre(nu)
    vs = VSpace1D(points[1], points[end], nu, points, ones(nu) .* (points[end] - points[1]) / (nu - 1), weights)

    # material
    σs = ones(Float32, nx)
    σa = zeros(Float32, nx)
    σt = σs + σa
    σq = zeros(Float32, nx)

    # moments
    L = 1
    ne = 2
    m = eval_sphermonomial(points, L)

    # time
    cfl = 0.5
    dt = cfl * ps.dx[1]
    nt = set.maxTime / dt |> floor |> Int
    global t = 0.0

    # solution
    f0 = 0.0001 * ones(nu)
    phi = zeros(ne, nx)
    for i = 1:nx
        phi[:, i] .= m * f0
    end
    α = zeros(Float32, ne, nx)

    # NN
    αT = zeros(Float32, nx, ne)
    phiT = zeros(Float32, nx, ne)
    phi_old = zeros(Float32, ne, nx)
    phi_temp = deepcopy(phi_old)

    opt = optimize_closure(zeros(Float32, 2), m, weights, phi[:, 1], KitBase.maxwell_boltzmann_dual)

    global X = zeros(Float32, ne, 1)
    X[:, 1] .= phi[:, 1]
    global Y = zeros(Float32, 1, 1) # h
    Y[1, 1] = kinetic_entropy(opt.minimizer, m, weights)
    #global Y = zeros(Float32, ne, 1) # α
    #Y[:, 1] = opt.minimizer
end

begin
    # initial condition
    f0 = 0.0001 * ones(nu)
    phi = zeros(ne, nx)
    for i = 1:nx
        phi[:, i] .= m * f0
    end
    α = zeros(Float32, ne, nx)
    flux = zeros(Float32, ne, nx + 1)
    fη = zeros(nu)
end

#anim = @animate
for iter = 1:nt
    println("iteration $iter of $nt")

    # mathematical optimizer
    @inbounds for i = 1:nx
        opt = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
        α[:, i] .= opt.minimizer
        phi[:, i] .= KitBase.realizable_reconstruct(opt.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

        #X = hcat(X, phi[:, i])
        #Y = hcat(Y, kinetic_entropy(α[:, i], m, weights))
    end

    flux_wall!(fη, maxwell_boltzmann_dual.(α[:, 1]' * m)[:], points, dt, 1.0)
    for k in axes(flux, 1)
        flux[k, 1] = sum(m[k, :] .* weights .* fη)
    end

    @inbounds for i = 2:nx
        KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, i-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i]' * m)[:], points, dt)

        for k in axes(flux, 1)
            flux[k, i] = sum(m[k, :] .* weights .* fη)
        end
    end

    @inbounds for i = 1:nx-1
        for q = 1:1
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end

        for q = 2:ne
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (-σt[i] * phi[q, i]) * dt
        end
    end
    phi[:, nx] .=  phi[:, nx-1]

    global t += dt

    #plot(ps.x[1:nx], phi[1, :])
end
#gif(anim, "newton1d.gif")

for i = 1:nx
    X = hcat(X, phi[:, i])
    Y = hcat(Y, kinetic_entropy(α[:, i], m, weights))
    #Y = hcat(Y, α[:, i])
end

X_old = deepcopy(X)
Y_old = deepcopy(Y)

fnn = FastChain(FastDense(2, 10, tanh), FastDense(10, 10, tanh), FastDense(10, 10, tanh), FastDense(10, 10, tanh), FastDense(10, 1))
res = KitML.sci_train(fnn, (X, Y); maxiters = 2000)
res = KitML.sci_train(fnn, (X, Y), res.u, LBFGS(); maxiters=1000)
#=
p = initial_params(fnn)
function loss(p)
    #norm(KitML.neural_closure(fnn, p, X) .- Y, 2)
    l = 0.0
    for i in axes(X, 2)
        l += sum(fnn(X[:, i], p) .- Y[1, i])
    end
    return l /= size(X, 2)
end
cb = function (p, l)
    println("loss: $(loss(p))")
    return false
end
res = KitML.sci_train(loss, p, ADAM(); maxiters=200, cb=cb)
res = KitML.sci_train(loss, res.u, LBFGS(); maxiters=400, cb=cb)
res = KitML.sci_train(loss, res.u, ADAM(); maxiters=4000, cb=cb)
=#
#cd(@__DIR__)
#tmodel = KitML.load_model("best_model.h5"; mode = :tf)
#tmodel.fit(permutedims(X), permutedims(Y), epochs=10)
#tmodel.predict(permutedims(phi[:, 1]))
#KitML.neural_closure(tmodel, permutedims(phi[:, 1]))

X = deepcopy(X_old)
Y = deepcopy(Y_old)
begin
    # initial condition
    f0 = 0.0001 * ones(nu)
    phi = zeros(ne, nx)
    for i = 1:nx
        phi[:, i] .= m * f0
    end
    α = zeros(Float32, ne, nx)
    flux = zeros(Float32, ne, nx + 1)
    fη = zeros(nu)
end

anim = @animate for iter = 1:nt
#for iter = 1:1#nt
    println("iteration $iter of $nt")
    phi_old .= phi

    # regularization
    #α .= KitML.neural_closure(fnn, res.u, phi_old)
    for i = 1:nx
        α[:, i] .= KitML.neural_closure(fnn, res.u, phi_old[:, i])
    end

    #α .= fnn(phi_old, res.u)
    @inbounds Threads.@threads for i = 1:nx
        phi_temp[:, i] .= KitBase.realizable_reconstruct(α[:, i], m, weights, KitBase.maxwell_boltzmann_dual_prime)
    end

    counter = 0
    @inbounds for i = 1:nx
        if norm(phi_temp[:, i] .- phi_old[:, i], 1) / phi_old[1, i] > 1e-2#2e-3
            counter +=1

            opt = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
            α[:, i] .= opt.minimizer
            phi[:, i] .= KitBase.realizable_reconstruct(opt.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

            X = hcat(X, phi[:, i])
            Y = hcat(Y, kinetic_entropy(α[:, i], m, weights))
            #Y = hcat(Y, α[:, i])
        else
            phi[:, i] .= phi_temp[:, i]
        end
    end
    println("newton: $counter of $nx")

    flux_wall!(fη, maxwell_boltzmann_dual.(α[:, 1]' * m)[:], points, dt, 1.0)
    for k in axes(flux, 1)
        flux[k, 1] = sum(m[k, :] .* weights .* fη)
    end

    @inbounds for i = 2:nx
        KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, i-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i]' * m)[:], points, dt)

        for k in axes(flux, 1)
            flux[k, i] = sum(m[k, :] .* weights .* fη)
        end
    end

    @inbounds for i = 1:nx-1
        for q = 1:1
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end

        for q = 2:ne
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (-σt[i] * phi[q, i]) * dt
        end
    end
    phi[:, nx] .=  phi[:, nx-1]

    global t += dt

    if iter%9 == 0
        #res = KitML.sci_train(fnn, (X, Y), res.u, ADAM(); maxiters=200)
        res = KitML.sci_train(fnn, (X, Y), res.u, LBFGS(); maxiters=400)
        #=
        X = zeros(Float32, ne, 1)
        #Y = zeros(Float32, 1, 1)
        Y = zeros(Float32, ne, 1)
        opt = KitBase.optimize_closure(α[:, nx], m, weights, phi[:, nx], KitBase.maxwell_boltzmann_dual)
        X[:, 1] .= phi[:, nx]
        #Y[1, 1] = kinetic_entropy(α[:, nx], m, weights)
        Y[:, 1] = opt.minimizer=#
    end

    plot(ps.x[1:nx], phi[1, :])
end

gif(anim, "unified_1d.gif")
#gif(anim, "neural_1d.gif")
plot(ps.x[1:nx], phi[1, :])

#=using BenchmarkTools

#@btime KitML.neural_closure(tmodel, phiT[:, :])
@btime fnn(phi, res.u)

@btime for i = 1:nx
    KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual)
end
=#

max_stars_repo_path: example/1D propagation/map_h.jl

max_issues_repo_path: example/1D propagation/map_h.jl

max_forks_repo_path: example/1D propagation/map_h.jl


text:
sumhappy(n) = sum(x->x^2, digits(n))

function ishappy(x, mem = [])
  x == 1?   true :
  x in mem? false :
  ishappy(sumhappy(x),[mem ; x])
end

nexthappy (x) = ishappy(x+1) ? x+1 : nexthappy(x+1)

happy(n) = [z = 1 ; [z = nexthappy(z) for i = 1:n-1]]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using MyExample1
using Test

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


text:
#!/usr/bin/env julia

function main()
    tokens = Channel{String}(32)
    Task() do
        for line in eachline(@static VERSION < v"0.6" ? STDIN : stdin)
            for token in split(chomp(line))
                put!(tokens, token)
            end
        end
        close(tokens)
    end |> schedule
    A = parse(Int, take!(tokens))
    B = parse(Int, take!(tokens))
    C = parse(Int, take!(tokens))
    if A + B >= C
        println("Yes")
    else
        println("No")
    end
end

isempty(ARGS) && main()


text:
using NaiveNASlib
using Test
import JuMP

include("testutil.jl")

@testset "NaiveNASlib.jl" begin

        @info "Testing computation"

        include("vertex.jl")
        include("compgraph.jl")

        include("prettyprint.jl")

        @info "Testing mutation"

        include("mutation/op.jl")
        include("mutation/vertex.jl")
        include("mutation/graph.jl")

        @info "Testing size mutation"

        include("mutation/size.jl")

        @info "Testing index mutation"

        include("mutation/apply.jl")
        include("mutation/select.jl")

        @info "Testing structural mutation"

        include("mutation/structure.jl")

        @info "Testing sugar"

        include("mutation/sugar.jl")

        @info "Testing README examples"

        include("examples.jl")
end


text:
function conv_bn(
    kernelsize,
    inplanes,
    outplanes,
    activation = relu;
    rev = false,
    initβ = Flux.zeros32,
    initγ = Flux.ones32,
    ϵ = 1f-5,
    momentum = 1f-1,
    kwargs...,
)
    layers = []

    if rev
        activations = (conv = activation, bn = identity)
        bnplanes = inplanes
    else
        activations = (conv = identity, bn = activation)
        bnplanes = outplanes
    end

    push!(
        layers,
        Conv(kernelsize, Int(inplanes) => Int(outplanes), activations.conv; kwargs...),
    )
    push!(
        layers,
        BatchNorm(
            Int(bnplanes),
            activations.bn;
            initβ = initβ,
            initγ = initγ,
            ϵ = ϵ,
            momentum = momentum,
        ),
    )

    return rev ? reverse(layers) : layers
end

basicblock(inplanes, outplanes, downsample = false) =
    downsample ?
    Chain(
        conv_bn((3, 3), inplanes, outplanes[1]; stride = 2, pad = 1, bias = false)...,
        conv_bn(
            (3, 3),
            outplanes[1],
            outplanes[2],
            identity;
            stride = 1,
            pad = 1,
            bias = false,
        )...,
    ) :
    Chain(
        conv_bn((3, 3), inplanes, outplanes[1]; stride = 1, pad = 1, bias = false)...,
        conv_bn(
            (3, 3),
            outplanes[1],
            outplanes[2],
            identity;
            stride = 1,
            pad = 1,
            bias = false,
        )...,
    )

bottleneck(inplanes, outplanes, downsample = false) =
    downsample ?
    Chain(
        conv_bn((1, 1), inplanes, outplanes[1]; stride = 2, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, bias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = 1, bias = false)...,
    ) :
    Chain(
        conv_bn((1, 1), inplanes, outplanes[1]; stride = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, bias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = 1, bias = false)...,
    )

skip_projection(inplanes, outplanes, downsample = false) =
    downsample ?
    Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 2, bias = false)...) :
    Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 1, bias = false)...)

function skip_identity(inplanes, outplanes)
    if outplanes > inplanes
        return Chain(
            MaxPool((1, 1), stride = 2),
            y -> cat(
                y,
                zeros(eltype(y), size(y, 1), size(y, 2), outplanes - inplanes, size(y, 4));
                dims = 3,
            ),
        )
    else
        return identity
    end
end
skip_identity(inplanes, outplanes, downsample) = skip_identity(inplanes, outplanes)

function resnet(
    block,
    residuals::NTuple{2,Any},
    connection = (x, y) -> @. relu(x) + relu(y);
    channel_config,
    block_config,
    nclasses = 1000,
)
    inplanes = 64
    baseplanes = 64
    layers = []
    append!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = (3, 3)))
    push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
    for (i, nrepeats) in enumerate(block_config)
        # output planes within a block
        outplanes = baseplanes .* channel_config
        # push first skip connection on using first residual
        # downsample the residual path if this is the first repetition of a block
        push!(
            layers,
            Parallel(
                connection,
                block(inplanes, outplanes, i != 1),
                residuals[1](inplanes, outplanes[end], i != 1),
            ),
        )
        # push remaining skip connections on using second residual
        inplanes = outplanes[end]
        for _ = 2:nrepeats
            push!(
                layers,
                Parallel(
                    connection,
                    block(inplanes, outplanes, false),
                    residuals[1](inplanes, outplanes[end], false),
                ),
            )
            inplanes = outplanes[end]
        end
        # next set of output plane base is doubled
        baseplanes *= 2
    end

    return Chain(
        Chain(layers...),
        Chain(AdaptiveMeanPool((1, 1)), flatten, Dense(inplanes, nclasses)),
    )
end

resnet(block, shortcut_config::Symbol, args...; kwargs...) =
    (shortcut_config == :A) ?
    resnet(block, (skip_identity, skip_identity), args...; kwargs...) :
    (shortcut_config == :B) ?
    resnet(block, (skip_projection, skip_identity), args...; kwargs...) :
    (shortcut_config == :C) ?
    resnet(block, (skip_projection, skip_projection), args...; kwargs...) :
    error(
        "Unrecognized shortcut config == $shortcut_config passed to resnet (use :A, :B, or :C).",
    )

const resnet_config = Dict(
    :resnet18 => ([1, 1], [2, 2, 2, 2], :A),
    :resnet34 => ([1, 1], [3, 4, 6, 3], :A),
    :resnet50 => ([1, 1, 4], [3, 4, 6, 3], :B),
    :resnet101 => ([1, 1, 4], [3, 4, 23, 3], :B),
    :resnet152 => ([1, 1, 4], [3, 8, 36, 3], :B),
)

struct ResNet
    layers
end

function ResNet(channel_config, block_config, shortcut_config; block, nclasses = 1000)
    layers = resnet(
        block,
        shortcut_config;
        channel_config = channel_config,
        block_config = block_config,
        nclasses = nclasses,
    )

    ResNet(layers)
end

@functor ResNet

(m::ResNet)(x) = m.layers(x)

backbone(m::ResNet) = m.layers[1]
classifier(m::ResNet) = m.layers[2]

function ResNet18(; pretrain = false, nclasses = 1000)
    model = ResNet(resnet_config[:resnet18]...; block = basicblock, nclasses = nclasses)
    return model
end

function ResNet34(; pretrain = false, nclasses = 1000)
    model = ResNet(resnet_config[:resnet34]...; block = basicblock, nclasses = nclasses)
    return model
end


text:
import ProximalOperators: gradient!, gradient # this can be removed when moved to Prox

export PrecomposeNonlinear

struct PrecomposeNonlinear{P <: ProximableFunction,
                           T <: AbstractOperator,
                           D <: AbstractArray,
                           C <: AbstractArray
                          } <: ProximableFunction
	g::P
  G::T
  bufD::D
  bufC::C
  bufC2::C
end

function PrecomposeNonlinear(g::P, G::T) where {P, T}
  t, s = domainType(G), size(G,2)
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s))
  t, s = codomainType(G), size(G,1)
  bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s))
  bufC2 = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s))
  PrecomposeNonlinear{P, T, typeof(bufD), typeof(bufC)}(g, G, bufD, bufC, bufC2)
end

is_smooth(f::PrecomposeNonlinear) = is_smooth(f.g)

function (f::PrecomposeNonlinear)(x)
    return f.g(f.G*x)
end

function gradient(f::PrecomposeNonlinear, x::ArrayPartition)
  y = zero(x)
  fy = gradient!(y,f,x)
  return y, fy
end

#TODO simplify this
function gradient!(y::D, f::PrecomposeNonlinear{P,T,D,C}, x::D) where {P,T,D <: ArrayPartition,C}
    mul!(f.bufC, f.G, x)
    v = gradient!(f.bufC2, f.g, f.bufC)
    J = Jacobian(f.G, x)
    y = mul!(y, J', f.bufC2)
    return v
end

function gradient!(y::D, f::PrecomposeNonlinear{P,T,D,C}, x::D) where {P,T,D <: AbstractArray,C}
    mul!(f.bufC, f.G, x)
    v = gradient!(f.bufC2, f.g, f.bufC)
    J = Jacobian(f.G, x)
    y = mul!(y, J', f.bufC2)
    return v
end


text:
##################################################
#MMMMMWxkOOOOOOOOOOOkxxxxxxxxxxxxxxkkkkkkkkkkkkNM#
#MMMMMN.'O0000000000d..........................0M#
#MMMMMX.'O0000000000x..........................0M#
#MMMMMX..xkxxxxxxxxxl.........................'0M#
#MMMMMX........................................0M#
#MMMMMX........................................0M#
#MMMMMX........................................0M#
#MMMMMX...................''...................0M#
#MMMMMX................,cx00kl,'...............0M#
#MMMMMX...............,xMMMMMM0,'..............0M#
#MMMMMX..............',0MMMMMMN,'..............0M#
#MMMMMX...............';kXWMNO:'...............0M#
#MMMMMX.................',,,,'.................0M#
#MMMMMX........................................0M#
#MMMMMX........................................0M#
#MMMMMN...................''...................0M#
#MMMMMN...................''...................OM#
#MMMMMN........................................OM#
#MMMMMN........................................OM#
#MMMMMW,,,,,,,,,,,,,;;;;:;;:;::;;;;;::::::::::cKM#
##################################################

# Caching.jl - a memory/disk memoizer written at 0x0α Research
#              by Corneliu Cofaru, 2018

module Caching

using Serialization
using DataStructures
using TranscodingStreams
using CodecZlib, CodecBzip2
using MacroTools
using Random
using Base: summarysize

import Base: show, empty!, length
import Core.Compiler: return_type

const MAX_CACHE_SIZE = typemax(Int)

export AbstractCache,
       Cache,
       AbstractSize,
       CountSize,
       MemorySize,
       object_size,
       max_cache_size,
       arghash,
       syncache!,
       persist!,
       empty!,
       @cache,
       @syncache!,
       @persist!,
       @empty!

include("cache.jl")
include("utils.jl")
include("file.jl")

end  # module


text:
α = 0.025

@testset "helpers" begin
  Random.seed!(43)
  @test pvalue(ExactOneSampleKSTest(rand2cop(rand(500000), 0.5, "clayton"), Uniform(0,1))) > α
  Random.seed!(43)
  @test pvalue(ExactOneSampleKSTest(rand2cop(rand(500000), 0.5, "frank"), Uniform(0,1))) > α
  Random.seed!(42)
  @test pvalue(ExactOneSampleKSTest(rand2cop(rand(500000), 0.5, "amh"), Uniform(0,1))) > α
  Random.seed!(43)
  @test rand2cop([0.815308, 0.894269], 0.5, "clayton") ≈ [0.292041, 0.836167] atol=1.0e-5
end

@testset "exceptions" begin
  @test_throws DomainError testbivθ(-2., "clayton")
  @test_throws DomainError usebivρ(-.9, "amh", "Spearman")
  @test_throws DomainError usebivρ(-.25, "amh", "Kendall")
  @test_throws AssertionError chaincopulagen(100000, [1.1, 1.6], "gumbel")
end

@testset "chain of Archimedean copulas" begin
  Random.seed!(43)
  cops = ["clayton", "clayton", "clayton", "frank", "amh", "amh"]
  x = chaincopulagen(500000, [-0.9, 3., 2, 4., -0.3, 1.], cops)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,7], Uniform(0,1))) > α
  @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "l") ≈ 0.5 atol=1.0e-1
  @test tail(x[:,5], x[:,6], "r", 0.0001) ≈ 0
  @test tail(x[:,3], x[:,4], "l") ≈ 1/(2^(1/2)) atol=1.0e-1
  @test tail(x[:,3], x[:,4], "r", 0.0001) ≈ 0
  @test corkendall(x)[1,2] ≈ -0.9/(2-0.9) atol=1.0e-3
  @test corkendall(x)[2,3] ≈ 3/(2+3) atol=1.0e-3
end
@testset "correlations" begin
  Random.seed!(43)
  x = chaincopulagen(500000, [0.6, -0.2], "clayton"; cor = "Spearman")
  @test corspearman(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-2
  @test corspearman(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-2
  Random.seed!(43)
  x = chaincopulagen(500000, [0.6, -0.2], "clayton"; cor = "Kendall")
  @test corkendall(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-3
  @test corkendall(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-3
end
@testset "rev copula" begin
  Random.seed!(43)
  x = chaincopulagen(500000, [-0.9, 3., 2.], "clayton"; rev = true)
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test tail(x[:,3], x[:,4], "r") ≈ 1/(2^(1/2)) atol=1.0e-1
  @test tail(x[:,3], x[:,4], "l", 0.0001) ≈ 0
end

@testset "chain of Frechet copulas" begin
  @test fncopulagen([0.2, 0.4], [0.1, 0.1], [0.2 0.4 0.6; 0.3 0.5 0.7]) == [0.6 0.4 0.2; 0.7 0.5 0.3]
  Random.seed!(43)
  x = chainfrechetcopulagen(500000, [0.9, 0.6, 0.2])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test corspearman(x) ≈ [1. 0.9 0.6 0.2; 0.9 1. 0.6 0.2; 0.6 0.6 1. 0.2; 0.2 0.2 0.2 1.] atol=1.0e-2
  @test tail(x[:,1], x[:,2], "r") ≈ 0.9 atol=1.0e-1
  @test tail(x[:,1], x[:,2], "l") ≈ 0.9 atol=1.0e-1
  @test tail(x[:,1], x[:,4], "r") ≈ 0.2 atol=1.0e-1
  Random.seed!(43)
  x = chainfrechetcopulagen(500000, [0.8, 0.5], [0.2, 0.3]);
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test corspearman(x) ≈ [1. 0.6 0.2; 0.6 1. 0.2; 0.2 0.2 1.] atol=1.0e-3
  @test tail(x[:,1], x[:,2], "r") ≈ 0.8 atol=1.0e-1
  @test tail(x[:,2], x[:,3], "r") ≈ 0.5 atol=1.0e-1
end


text:
using ArgParse,Match
using NIPALS_PCA

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin

        "--xfile"
            help = "x data in delimited file with the first row containing variable names and the first column containing observation names"
            required = true

        "--yfile"
            help = "y data in delimited file with the first row containing variable names and the first column containing observation names"

        "--ycontinous"
            help = "y variable names to include for model calibration. Multiple entries are separated by ;, example:\"DV200;RIN\""

        "--ycategorical"
            help = "y variable names to include for model calibration"

        "--mode"
            help = "calibrate or correct"
            required = true

        "--uvscale"
            help = "scale variables for unit variance"
            arg_type = Bool
            default = true

        "--modelfile"
            help = "name of model file to read or write"
            required = true

        "--components", "-A"
            help = "specify the number of components to use"
            arg_type = Int
            default = 3

        "--mvcutoff"
            help = "cutoff for excluding variables from modelling. Values with missingness <= mvcutoff is gets included"
            arg_type = Float64
            default = 0.25

        "--minvalues"
            help = "the number of values that a variable needs to have to not be excluded, numvalues >= minvalues gets included"
            arg_type = Int
            default = 5

         "--outfile", "-o"
            help = "output corrected data filename"
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    parsed_args |> println

    xfile = parsed_args["xfile"]

    @match parsed_args["mode"] begin
        "calibrate" => calibrate_model(parsed_args)
        "correct" => correct(parsed_args)
    end
end

main()


text:
status,len,point,buffer = Grisu.fastfixedtoa(1.0, 0,1, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.0, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.0, 0,0, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0xFFFFFFFF, 0,5, buffer)
@test "4294967295" == bytestring(pointer(buffer))
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(4294967296.0, 0,5, buffer)
@test "4294967296" == bytestring(pointer(buffer)) #todo
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e21, 0,5, buffer)
@test "1" == bytestring(pointer(buffer)) #todo extra '0's
@test 22 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(999999999999999868928.00, 0,2, buffer)
@test "999999999999999868928" == bytestring(pointer(buffer)) #todo extra '0'
@test 21 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(6.9999999999999989514240000e+21, 0,5, buffer)
@test "6999999999999998951424" == bytestring(pointer(buffer)) #todo short several '9's
@test 22 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.5, 0,5, buffer)
@test "15" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.55, 0,5, buffer)
@test "155" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.55, 0,1, buffer)
@test "16" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.00000001, 0,15, buffer)
@test "100000001" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.1, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test 0 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.01, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.001, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0001, 0,10, buffer) #todo
@test "1" == bytestring(pointer(buffer))
@test -3 == point
fill!(buffer,0);

v = 0.00001; fractional_count = 10
status,len,point,buffer = Grisu.fastfixedtoa(0.00001, 0,10, buffer) #todo
@test "1" == bytestring(pointer(buffer))
@test -4 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000001, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -5 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000001, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -6 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000001, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -7 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000001, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000001, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -9 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000001, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000001, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -11 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000001, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -12 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000000001, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -13 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -14 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -15 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -16 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -17 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -18 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000000000000001, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -19 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.10000000004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test 0 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.01000000004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00100000004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00010000004, 0,10, buffer) #todo
@test "1" == bytestring(pointer(buffer))
@test -3 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00001000004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -4 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000100004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -5 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000010004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -6 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000001004, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -7 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000104, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test -8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000001000004, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -9 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000100004, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000010004, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -11 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000001004, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -12 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000000104, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test -13 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000001000004, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -14 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000100004, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -15 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000010004, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -16 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000001004, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -17 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000000104, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -18 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000000014, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -19 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.10000000006, 0,10, buffer)
@test "1000000001" == bytestring(pointer(buffer))
@test 0 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.01000000006, 0,10, buffer)
@test "100000001" == bytestring(pointer(buffer))
@test -1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00100000006, 0,10, buffer)
@test "10000001" == bytestring(pointer(buffer))
@test -2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00010000006, 0,10, buffer)
@test "1000001" == bytestring(pointer(buffer))
@test -3 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00001000006, 0,10, buffer)
@test "100001" == bytestring(pointer(buffer))
@test -4 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000100006, 0,10, buffer)
@test "10001" == bytestring(pointer(buffer))
@test -5 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000010006, 0,10, buffer)
@test "1001" == bytestring(pointer(buffer))
@test -6 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000001006, 0,10, buffer)
@test "101" == bytestring(pointer(buffer))
@test -7 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000000106, 0,10, buffer)
@test "11" == bytestring(pointer(buffer))
@test -8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000001000006, 0,15, buffer)
@test "100001" == bytestring(pointer(buffer))
@test -9 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000100006, 0,15, buffer)
@test "10001" == bytestring(pointer(buffer))
@test -10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000010006, 0,15, buffer)
@test "1001" == bytestring(pointer(buffer))
@test -11 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000001006, 0,15, buffer)
@test "101" == bytestring(pointer(buffer))
@test -12 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000000000000106, 0,15, buffer)
@test "11" == bytestring(pointer(buffer))
@test -13 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000001000006, 0,20, buffer)
@test "100001" == bytestring(pointer(buffer))
@test -14 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000100006, 0,20, buffer)
@test "10001" == bytestring(pointer(buffer))
@test -15 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000010006, 0,20, buffer)
@test "1001" == bytestring(pointer(buffer))
@test -16 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000001006, 0,20, buffer)
@test "101" == bytestring(pointer(buffer))
@test -17 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000000106, 0,20, buffer)
@test "11" == bytestring(pointer(buffer))
@test -18 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000000000000000016, 0,20, buffer)
@test "2" == bytestring(pointer(buffer))
@test -19 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.6, 0,0, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.96, 0,1, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.996, 0,2, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.9996, 0,3, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.99996, 0,4, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.999996, 0,5, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.9999996, 0,6, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.99999996, 0,7, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.999999996, 0,8, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.9999999996, 0,9, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.99999999996, 0,10, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.999999999996, 0,11, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.9999999999996, 0,12, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.99999999999996, 0,13, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.999999999999996, 0,14, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.9999999999999996, 0,15, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00999999999999996, 0,16, buffer)
@test "1" == bytestring(pointer(buffer))
@test -1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000999999999999996, 0,17, buffer)
@test "1" == bytestring(pointer(buffer))
@test -2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.0000999999999999996, 0,18, buffer)
@test "1" == bytestring(pointer(buffer))
@test -3 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.00000999999999999996, 0,19, buffer)
@test "1" == bytestring(pointer(buffer))
@test -4 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.000000999999999999996, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -5 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(323423.234234, 0,10, buffer)
@test "323423234234" == bytestring(pointer(buffer))
@test 6 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(12345678.901234, 0,4, buffer)
@test "123456789012" == bytestring(pointer(buffer))
@test 8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(98765.432109, 0,5, buffer)
@test "9876543211" == bytestring(pointer(buffer))
@test 5 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(42, 0,20, buffer)
@test "42" == bytestring(pointer(buffer))
@test 2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(0.5, 0,0, buffer)
@test "1" == bytestring(pointer(buffer))
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-23, 0,10, buffer)
@test "" == bytestring(pointer(buffer))
@test -10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-123, 0,2, buffer)
@test "" == bytestring(pointer(buffer))
@test -2 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-123, 0,0, buffer)
@test "" == bytestring(pointer(buffer))
@test 0 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-23, 0,20, buffer)
@test "" == bytestring(pointer(buffer))
@test -20 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-21, 0,20, buffer)
@test "" == bytestring(pointer(buffer))
@test -20 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1e-22, 0,20, buffer)
@test "" == bytestring(pointer(buffer))
@test -20 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(6e-21, 0,20, buffer)
@test "1" == bytestring(pointer(buffer))
@test -19 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(9.1193616301674545152000000e+19, 0,0,buffer)
@test "91193616301674545152" == bytestring(pointer(buffer))
@test 20 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(4.8184662102767651659096515e-04, 0,19,buffer)
@test "4818466210276765" == bytestring(pointer(buffer))
@test -3 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1.9023164229540652612705182e-23, 0,8,buffer)
@test "" == bytestring(pointer(buffer))
@test -8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastfixedtoa(1000000000000000128.0, 0,0,buffer)
@test "1000000000000000128" == bytestring(pointer(buffer))
@test 19 == point
fill!(buffer,0);