text:
#!/home/fng514/bin/julia
# Christian Theil Have, 2017.
# Merging of imputed data.
#
# Merges a selected set of individuals from a pair of imputations on the same panel.
# If the same individual is present in more both imputations, then the script will
# select the best-possible imputed genotypes for the individuals from either imputation
# on a per genotype basis. INFO scores as well as AN (allele number) and (AC) allele
# count is recalculated and the output VCF file is written with the same fields.
#
# It is not as generic as I would like it to be:
# Currently, designed to work with the format produced by Sangers imputation server.
# Adaptation will be need for other formats.

using ArgParse
function main(args)
	s = ArgParseSettings("Example for merge.jl: " *  # description
		"flags, options help, " *
		"required arguments.")

	@add_arg_table s begin
		"--vcfs"
		nargs = '*'
		help = "A list of VCF files to merge"
		"--merge-mode"
		help = "One of best,mean,weightedmean. Default is best."
		default = "best"
		"--out"
	end

	parsed_args = parse_args(s)
	println("Parsed args:")
	for (key,val) in parsed_args
		println("  $key  =>  $(repr(val))")
	end

	@assert length(parsed_args["vcfs"]) > 0

	# Merge method:
	if parsed_args["merge-mode"] == "best"
		mergefun = best_gp
	elseif parsed_args["merge-mode"] == "mean"
		mergefun = mean_gp
	elseif parsed_args["merge-mode"] == "weightedmean"
		mergefun = confidence_weighted_mean_gp
	elseif parsed_args["merge-mode"] == "bestinfo"
		mergefun = best_info
	else
		throw("Invalid merge-mode : $(parsed_args["merge-mode"])")
	end

	merged_imputed(parsed_args["vcfs"], parsed_args["out"], mergefun)
end

main(ARGS)

max_stars_repo_name: cth/merge-imputations

max_issues_repo_name: cth/merge-imputations

max_forks_repo_name: cth/merge-imputations


text:
# show.jl

function showfull(io::IO, f::Union{TaylorModel1, RTaylorModel1, TaylorModelN})
    print(io,
    """Taylor model of degree $(get_order(f)):
    - x0:  $(f.x0)
    -  I:  $(f.I)
    -  p:  $(f.pol)
    -  Δ:  $(f.rem)
    """
    )
end
showfull(x) = showfull(stdout::IO, x)

function show(io::IO, a::Union{TaylorModel1, RTaylorModel1, TaylorModelN})
    if TaylorSeries._show_default[end]
        return Base.show_default(io, a)
    else
        return print(io, pretty_print(a))
    end
end

for T in (:TaylorModel1, :TaylorModelN, :RTaylorModel1)
    @eval function pretty_print(a::$T)
        _bigOnotation = TaylorSeries.bigOnotation[end]
        _bigOnotation && TaylorSeries.displayBigO(false)
        strout = $T == TaylorModelN ?
            string(a.pol, " + ", a.rem) : string(a.pol, "+ ", a.rem)
        if $T == RTaylorModel1
            _order = get_order(a)
            strout = strout * " t" * TaylorSeries.superscriptify(_order+1)
        end
        _bigOnotation && TaylorSeries.displayBigO(true)
        strout
    end
end

ext: jl


text:
function Topo(x,y,zeta,Param)
if Param.TopoS == "AgnesiCartX"
    h=Param.hC/(((x-Param.x0C)/Param.aC)^2+1);
elseif Param.TopoS == "AgnesiCartY"
    h=Param.hC/(((y-Param.y0C)/Param.aC)^2+1);
else
    h=0;
end
Z=zeta+(Param.H-zeta)*h/Param.H;
dZ=1-h/Param.H;
return (Z,dZ)
end


text:
"""
A linear operator corresponding to convolution with either s or g
"""
function operator(pa, attrib)

	fw=(y,x)->F!(y, x, pa, attrib)
	bk=(y,x)->Fadj!(y, x, pa, attrib)

	return LinearMap{collect(typeof(pa).parameters)[1]}(fw, bk,
		  length(pa.d),  # length of output
		  length(pa.g),  # length of output
		  ismutating=true)
end

function F!(y, x, pa, ::S)
	copyto!(pa.g,x)
	mod!(pa, D())
	copyto!(y,pa.d)
	return nothing
end

function Fadj!(y, x, pa, ::S)
	copyto!(pa.d,x)
	mod!(pa, G())
	copyto!(y,pa.g)
	return nothing
end

function F!(y, x, pa, ::G)
	copyto!(pa.s,x)
	mod!(pa, D())
	copyto!(y,pa.d)
	return nothing
end

function Fadj!(y, x, pa, ::G)
	copyto!(pa.d,x)
	mod!(pa, S())
	copyto!(y,pa.s)
	return nothing
end

function F!(y, x, pa, ::D)
	copyto!(pa.g,x)
	mod!(pa, S())
	copyto!(y,pa.s)
	return nothing
end

function Fadj!(y, x, pa, ::D)
	copyto!(pa.s,x)
	mod!(pa, G())
	copyto!(y,pa.g)
	return nothing
end


text:
@testset "broadcast" begin
    @testset "Misc. Tests" begin
        @testset "sin.(x)" begin
            x = rand(3, 3)
            y, (dsin, dx) = rrule(broadcast, sin, x)

            @test y == sin.(x)
            @test extern(dx(One())) == cos.(x)

            x̄, ȳ = rand(), rand()
            @test extern(accumulate(x̄, dx, ȳ)) == x̄ .+ ȳ .* cos.(x)

            x̄, ȳ = Zero(), rand(3, 3)
            @test extern(accumulate(x̄, dx, ȳ)) == ȳ .* cos.(x)

            x̄, ȳ = Zero(), cast(rand(3, 3))
            @test extern(accumulate(x̄, dx, ȳ)) == extern(ȳ) .* cos.(x)
        end
    end
end


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

mutable struct IoK8sApiFlowcontrolV1alpha1LimitResponse <: SwaggerModel
    queuing::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1QueuingConfiguration } # spec name: queuing
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiFlowcontrolV1alpha1LimitResponse(;queuing=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitResponse, Symbol("queuing"), queuing)
        setfield!(o, Symbol("queuing"), queuing)
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitResponse, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1LimitResponse

const _property_map_IoK8sApiFlowcontrolV1alpha1LimitResponse = Dict{Symbol,Symbol}(Symbol("queuing")=>Symbol("queuing"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiFlowcontrolV1alpha1LimitResponse = Dict{Symbol,String}(Symbol("queuing")=>"IoK8sApiFlowcontrolV1alpha1QueuingConfiguration", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1LimitResponse }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1LimitResponse))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1LimitResponse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1LimitResponse[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1LimitResponse }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1LimitResponse[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1LimitResponse)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

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


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using StatsBase: sample, Weights

function episode!(actor, episode, env, maxt)
    reset!(env)
    for _ in 1:maxt
        s = copy(env.current)
        a = sample(Weights(actor(s)))
        _, r, done = step!(env, a)
        append!(episode, s, a, r)
        done && break
    end
end

function run!(actor, opt, env; epochs=500, maxt=100, γ=99f-2)
    rewards = zeros(Float32, div(epochs, 10))
    episode = PGEpisode{Float32, Int64, Float32}(length(env.observationspace), 1, maxt, γ=γ)

    for i in 1:epochs
        episode!(actor, episode, env, maxt)
        discount!(episode)

        Flux.train!(params(actor), [batch(episode)], opt) do s, a, r, _
            -sum(log.(gather(actor(s), a)) .* r)
        end
        i % 10 == 0 && (rewards[div(i, 10)] = test(s -> Flux.onecold(actor(s)), env))
    end
    return rewards
end

env = CartPole{Float32}()

actor = Chain(
    Dense(length(env.observationspace), 32, relu),
    Dense(32, length(env.actionspace)),
    softmax
)

r = run!(actor, ADAM(0.01, (0.99, 0.999)), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "reinforce")


text:
type SCF
    Etol::Float64
    rhotol::Float64
    maxsteps::Int  # max num of iterations
    n::Int  # how many previous iterations to store
    grid  # integration grid
    i::Int  # iteration count
    converged::Bool
    E::Array{Float64}  #
    rho::Array{AbstractArray}  #
    dE::Float64
    drho::Float64
end

SCF(;Etol=1e-6, rhotol=1e-8, maxsteps=100, n::Int=1, grid=nothing) = SCF(Etol,rhotol,maxsteps,n,grid,1,false,[],[],NaN,NaN)

function Base.call(scf::SCF, E::Float64, rho::AbstractArray)
    if length(scf.E) > 0
        scf.dE = abs(E - scf.E[end])
        scf.drho = integrate(scf.grid, (rho - scf.rho[end]).^2)
    end

    if scf.dE < scf.Etol && scf.drho < scf.rhotol
        scf.converged = true
    end

    scf.i += 1

    if scf.i > scf.maxsteps
        warn("max # of iterations exceeded")
        scf.converged = true
    end

    push!(scf.rho, rho)
    push!(scf.E, E)

    if scf.n > 0 && length(scf.rho) > scf.n
        shift!(scf.rho)
        shift!(scf.E)
    end
end

simple_logger = scf -> !isnan(scf.dE) ? @printf("|δE| = %8.2e     |δρ| = %8.2e\n", scf.dE, scf.drho) : return

max_stars_repo_name: johncsnyder/Element


text:
#It would be nice to automate downloading of these files
#also may want some routines to process the metadata and flat files
#does pre-processing of Open990 files

#=
uses constants
const OPEN990_PATH = DATA_PATH * "\\open990"
const OPEN990_FIELD_PATH = OPEN990_PATH * "\\fieldfiles"
const OPEN990_DATE_FORMAT = "yyyymm"
const OPEN990_FIELD_INDEX_NAME = "fieldfileindex"
=#

#loads the Open990 index
function loadfieldfileindex(;path::String = OPEN990_PATH,
  indexname::String = OPEN990_FIELD_INDEX_NAME)

  return CSV.read("$path\\$indexname.csv") |> DataFrame
end

#reads a single file
#takes in the file name, the targeted field name (for grouping)
#and the precedence (higher numbers have greater precedence)
function readfieldfile(filename::String;
  istreamer::Function = I_CSV_STREAMER,
  path::String = OPEN990_FIELD_PATH,
  dateformat::String = OPEN990_DATE_FORMAT,
  #submittedondateformat::String = OPEN990_SUBMITTED__ON_DATE_FORMAT,
  checkunique=true)::DataFrame

  fullpath::String = "$path\\$filename.csv.gz"
  b = istreamer(fullpath) do istream
    b = IOBuffer(read(istream))
    (b)
  end

  try
    df = CSV.read(b, strict=false, header=2,
    types = Dict(:last_updated=>String)) |> DataFrame
  catch err
    println("ERROR: $filename.csv.gz")
    error(err)
  end
  #print(describe(df))

  #now do some light pre-processing

  #strips the month off the fiscal year
  df[:fisyr] = (d::Int->d ÷ 100).(df[:tax_period])
  df[:datefiscal] = (d::Int->Date("$d", dateformat)).(df[:tax_period])

  #properly format the ein
  rename!(df, :ein=>:ein_old)
  df[:ein] = (i::Int->lpad(string(i),9,'0')).(df[:ein_old])

  rename!(df, :submitted_on=>:submitted_on_old)
  df[:submitted_on] = Vector{Date}(df[:submitted_on_old])

  #get unique ropws by taking the newest
  dfmaxsubmitted = by(df, [:ein, :fisyr], submitted_on = :submitted_on => maximum)
  df = join(df, dfmaxsubmitted, on = [:ein, :fisyr, :submitted_on])

  println("Read $(size(df,1)) records from $filename")
  #delete columns which will not be used

  deletecols!(df, setdiff(names(df), [:ein, :value, :fisyr, :submitted_on, :datefiscal]))

  #a spot test of the remaining dups shows ver similar, mostly identical
  #values, so just take the mean
  df = by(df, [:ein, :fisyr, :submitted_on, :datefiscal], value = :value => mean)

  #record the file-level metadata
  df[:filename] = filename

  return df
end

#reads in multiple files
function readfieldfiles(indexdf::DataFrame)::DataFrame

  local futures::Vector{Future} = Vector{Future}()
  local fielddfs::Vector{DataFrame} = Vector{DataFrame}()
  local mappingnames = unique(indexdf[:mappingname])

  for r ∈ eachrow(indexdf)
    push!(futures,
      @spawn readfieldfile(r[:filename]))
  end

  fielddflong::DataFrame = [(fetch).(futures)...;]

  println("Total records read: $(size(fielddflong))")
  fielddf::DataFrame = unstack(fielddflong, :filename, :value)

  #get rid of some dups by taking the newest entries
  dfmax = by(fielddf, [:ein, :fisyr], submitted_on = :submitted_on => maximum)
  fielddf = join(fielddf, dfmax, on = [:ein, :fisyr, :submitted_on])

  #get rid of the remainder by taking the latest fiscal year
  dfmax = by(fielddf, [:ein, :fisyr], datefiscal = :datefiscal=>maximum)
  fielddf = join(fielddf, dfmax, on = [:ein, :fisyr, :datefiscal])

  deletecols!(fielddf, [:datefiscal, :submitted_on])

  #NOTE: I checked and there are a few entries with multiple fisyrs. I elect
  #to select the entries with the newest fiscal year

  @assert size(unique(fielddf[[:ein, :fisyr]]),1) == size(fielddf,1)

  println("Total org-years: $(size(fielddf))")

  #println("uniquesize: $(size(unique(fielddf[[:ein, :fisyr]])))")

  return fielddf
end

function mapfieldfiles!(indexdf::DataFrame, fielddf::DataFrame)::Nothing

  #first build up indices of the mapping

  #println(describe(fielddf))
  mappingfields::Vector{Symbol} = (Symbol).(unique(indexdf[:mappingname]))
  Nmap::Int = length(mappingfields)

  #index that provides teh precedence
  #=precedenceindex::Dict =
    Dict(indexdf[i,:filename]=>indexdf[i, :precedence]
      for i ∈ 1:size(indexdf,1))=#

  #index that gives all of the filenames in order of precedence
  sort!(indexdf, (:mappingname, order(:precedence, rev=true)))
  filenameindex::Dict = Dict{Symbol, Vector{Symbol}}()
  lowprecedencevalues::Vector{Float64} = Vector{Float64}(undef, Nmap)
  for subdf ∈ groupby(indexdf, :mappingname)
    mapfield::Symbol = Symbol(subdf[1,:mappingname])
    filenameindex[mapfield] = (Symbol).(subdf[:filename])

    imapping::Int = findfirst(isequal(mapfield), mappingfields)
    lowprecedencevalues[imapping] = subdf[1,:lowprecedencevalue]
  end

  #allocate the mappingfields. Initially store the highest precedence
  for f ∈ mappingfields
    fielddf[f] = deepcopy(fielddf[filenameindex[f][1]])
  end

  #now do the mapping. Deafult to the first value.
  for r ∈ eachrow(fielddf) #iterate down each row in the new data
    for imapping ∈ 1:Nmap
      mapfield::Symbol = mappingfields[imapping]

      for filename ∈ filenameindex[mapfield][2:end] #iterate in order of precedence

        #if the candidate field is not missing, replace the old field if the old field is missing
        #or its the low precedence value
        if (!ismissing(r[filename])) && (ismissing(r[mapfield]) ||
          (lowprecedencevalues[imapping] == r[mapfield]))

          r[mapfield] = r[filename]
        end
      end
    end
  end

  #CSV.write("$WORKING_PATH\\open990test.csv", fielddf[1:100_000, :])

  #delete the old columns
  for f ∈ mappingfields
    deletecols!(fielddf, filenameindex[f])
  end

  #println(describe(fielddf))
  return nothing
end

function avg(v)::MFloat64
  local f::MFloat64
  if !isempty(skipmissing(v))
    f = mean(skipmissing(v))
    f = isfinite(f) ? f : missing
  else
    f = missing
  end

  return f
end

#does any aggregation functions that are needed and rejoins them to the main dataframe if desired
function aggregatefielddf(fielddf::DataFrame; rejoin::Bool=true,
  functions::Vector = [avg])

  local avgdf::DataFrame = aggregate(
    fielddf[setdiff(names(fielddf), [:fisyr])], :ein, functions)

  if rejoin
    fielddf = join(fielddf, avgdf, on=[:ein], kind=:left)
    return fielddf
  end

  return avgdf
end

#creates the full dataframe if needed
function createfielddf()
  indexdf::DataFrame = loadfieldfileindex()
  fielddf = readfieldfiles(indexdf)
  mapfieldfiles!(indexdf, fielddf)
  fielddf = aggregatefielddf(fielddf)

  return fielddf
end

#spits out the computationally heavy componeny sof the work into a compressed jls
function getfielddf(; path::String = OPEN990_PATH,
  fielddataname::String = OPEN990_FIELD_DATA_NAME,
  refreshopen990data::Bool = true,
  ostreamer::Function = O_JLS_STREAMER,
  istreamer::Function = I_JLS_STREAMER)

  local fielddf::DataFrame

  local fulldatapath::String = "$path\\$fielddataname.jls.gz"

  #recreate the file if none exists or a refresh is called
  if (!isfile(fulldatapath)) || refreshopen990data
    fielddf = createfielddf()

    ostream = ostreamer(fulldatapath)
    serialize(ostream, fielddf)
    close(ostream)
  else
    istream = istreamer(fulldatapath)
    fielddf = deserialize(istream)
    close(istream)
  end

  return fielddf
end

function mergeopen990data!(data::NCCSData; refreshopen990data= true)::Nothing

  fielddf::DataFrame = getfielddf(refreshopen990data=refreshopen990data)
  println(describe(fielddf))
  data.df = join(data.df, fielddf, on=[:ein, :fisyr], kind=:left)

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# interface for skipping null entries
import Compat: @functorize

function skipnull_init(f, op, X::NullableArray,
                       ifirst::Int, ilast::Int)
    # Get first non-null element
    ifirst = Base.findnext(x -> x == false, X.isnull, ifirst)
    @inbounds v1 = X.values[ifirst]

    # Get next non-null element
    ifirst = Base.findnext(x -> x == false, X.isnull, ifirst + 1)
    @inbounds v2 = X.values[ifirst]

    # Reduce first two elements
    return op(f(v1), f(v2)), ifirst
end

# sequential map-reduce
function mapreduce_seq_impl_skipnull(f, op, X::NullableArray,
                                     ifirst::Int, ilast::Int)
    # initialize first reduction
    v, i = skipnull_init(f, op, X, ifirst, ilast)

    while i < ilast
        i += 1
        @inbounds isnull = X.isnull[i]
        isnull && continue
        @inbounds entry = X.values[i]
        v = op(v, f(entry))
    end
    return Nullable(v)
end

# pairwise map-reduce
function mapreduce_pairwise_impl_skipnull{T}(f, op, X::NullableArray{T},
                                             ifirst::Int, ilast::Int,
                                            #  n_notnull::Int, blksize::Int)
                                            blksize::Int)
    if ifirst + blksize > ilast
        # fall back to Base implementation if no nulls in block
        # if anynull(slice(X, ifirst:ilast))
            return mapreduce_seq_impl_skipnull(f, op, X, ifirst, ilast)
        # else
            # Nullable(Base.mapreduce_seq_impl(f, op, X.values, ifirst, ilast))
        # end
    else
        imid = (ifirst + ilast) >>> 1
        # n_notnull1 = imid - ifirst + 1 - countnz(X.isnull[ifirst:imid])
        # n_notnull2 = ilast - imid - countnz(X.isnull[imid+1:ilast])
        v1 = mapreduce_pairwise_impl_skipnull(f, op, X, ifirst, imid,
                                              blksize)
        v2 = mapreduce_pairwise_impl_skipnull(f, op, X, imid+1, ilast,
                                              blksize)
        return op(v1, v2)
    end
end

# from comment: https://github.com/JuliaLang/julia/pull/16217#issuecomment-223768129
if VERSION < v"0.5.0-dev+4441"
    sum_pairwise_blocksize = Base.sum_pairwise_blocksize
else
    sum_pairwise_blocksize(T) = Base.pairwise_blocksize(T, +)
end

mapreduce_impl_skipnull{T}(f, op, X::NullableArray{T}) =
    mapreduce_seq_impl_skipnull(f, op, X, 1, length(X.values))
mapreduce_impl_skipnull(f, op::typeof(@functorize(+)), X::NullableArray) =
    mapreduce_pairwise_impl_skipnull(f, op, X, 1, length(X.values),
                                   max(128, sum_pairwise_blocksize(f)))

# general mapreduce interface

function _mapreduce_skipnull{T}(f, op, X::NullableArray{T}, missingdata::Bool)
    n = length(X)
    !missingdata && return Nullable(Base.mapreduce_impl(f, op, X.values, 1, n))

    nnull = countnz(X.isnull)
    nnull == n && return Base.mr_empty(f, op, T)
    nnull == n - 1 && return Base.r_promote(op, f(X.values[findnext(x -> x == false), X, 1]))
    # nnull == 0 && return Base.mapreduce_impl(f, op, X, 1, n)

    return mapreduce_impl_skipnull(f, op, X)
end

function Base._mapreduce(f, op, X::NullableArray, missingdata)
    missingdata && return Base._mapreduce(f, op, X)
    Nullable(Base._mapreduce(f, op, X.values))
end

# to fix ambiguity warnings
function Base.mapreduce(f, op::Union{typeof(@functorize(&)), typeof(@functorize(|))},
                        X::NullableArray, skipnull::Bool = false)
    missingdata = anynull(X)
    if skipnull
        return _mapreduce_skipnull(f, op, X, missingdata)
    else
        return Base._mapreduce(f, op, X, missingdata)
    end
end

if VERSION >= v"0.5.0-dev+3701"
    const specialized_binary = identity
else
    const specialized_binary = Base.specialized_binary
end

@doc """
`mapreduce(f, op::Function, X::NullableArray; [skipnull::Bool=false])`

Map a function `f` over the elements of `X` and reduce the result under the
operation `op`. One can set the behavior of this method to skip the null entries
of `X` by setting the keyword argument `skipnull` equal to true. If `skipnull`
behavior is enabled, `f` will be automatically lifted over the elements of `X`.
Note that, in general, mapreducing over a `NullableArray` will return a
`Nullable` object regardless of whether `skipnull` is set to `true` or `false`.
""" ->
function Base.mapreduce(f, op::Function, X::NullableArray;
                        skipnull::Bool = false)
    missingdata = anynull(X)
    if skipnull
        return _mapreduce_skipnull(f, specialized_binary(op),
                                   X, missingdata)
    else
        return Base._mapreduce(f, specialized_binary(op), X, missingdata)
    end
end

function Base.mapreduce(f, op, X::NullableArray; skipnull::Bool = false)
    missingdata = anynull(X)
    if skipnull
        return _mapreduce_skipnull(f, op, X, missingdata)
    else
        return Base._mapreduce(f, op, X, missingdata)
    end
end

@doc """
`mapreduce(f, op::Function, X::NullableArray; [skipnull::Bool=false])`

Reduce `X`under the operation `op`. One can set the behavior of this method to
skip the null entries of `X` by setting the keyword argument `skipnull` equal
to true. If `skipnull` behavior is enabled, `f` will be automatically lifted
over the elements of `X`. Note that, in general, mapreducing over a
`NullableArray` will return a `Nullable` object regardless of whether `skipnull`
is set to `true` or `false`.
""" ->
Base.reduce(op, X::NullableArray; skipnull::Bool = false) =
    mapreduce(@functorize(identity), op, X; skipnull = skipnull)

# standard reductions

for (fn, op) in ((:(Base.sum), @functorize(+)),
                 (:(Base.prod), @functorize(*)),
                 (:(Base.minimum), @functorize(scalarmin)),
                 (:(Base.maximum), @functorize(scalarmax)))
    @eval begin
        # supertype(typeof(@functorize(abs))) returns Func{1} on Julia 0.4,
        # and Function on 0.5
        $fn(f::Union{Function,supertype(typeof(@functorize(abs)))},
            X::NullableArray;
            skipnull::Bool = false) =
                mapreduce(f, $op, X; skipnull = skipnull)
        $fn(X::NullableArray; skipnull::Bool = false) =
            mapreduce(@functorize(identity), $op, X; skipnull = skipnull)
    end
end

for (fn, f, op) in ((:(Base.sumabs), @functorize(abs), @functorize(+)),
                    (:(Base.sumabs2), @functorize(abs2), @functorize(+)))
    @eval $fn(X::NullableArray; skipnull::Bool = false) =
        mapreduce($f, $op, X; skipnull=skipnull)
end

# internal methods for Base.minimum and Base.maximum
for op in (@functorize(scalarmin), @functorize(scalarmax))
    @eval begin
        function Base._mapreduce{T}(::typeof(@functorize(identity)), ::$(typeof(op)),
                                    X::NullableArray{T}, missingdata)
            missingdata && return Nullable{T}()
            Nullable(Base._mapreduce(@functorize(identity), $op, X.values))
        end
    end
end

function Base.mapreduce_impl{T}(f, op::typeof(@functorize(scalarmin)), X::NullableArray{T},
                                first::Int, last::Int)
    i = first
    v = f(X[i])
    i += 1
    while i <= last
        @inbounds x = f(X[i])
        if x.isnull || v.isnull
            return Nullable{eltype(x)}()
        elseif x.value < v.value
            v = x
        end
        i += 1
    end
    return v
end

function Base.mapreduce_impl{T}(f, op::typeof(@functorize(scalarmax)), X::NullableArray{T},
                                first::Int, last::Int)
    i = first
    v = f(X[i])
    i += 1
    while i <= last
        @inbounds x = f(X[i])
        if x.isnull || v.isnull
            return Nullable{eltype(x)}()
        elseif x.value > v.value
            v = x
        end
        i += 1
    end
    return v
end


text:
using PGFPlots
using CSV
using DataFrames

res = CSV.read("results.csv")

pcoll = Axis()
for (i,pol) in enumerate(res[:policy])
    push!(pcoll, Plots.BarChart([Float64.(res[:c_rate][i])], style="blue", legendentry=pol))
end

pavg = Axis(Plots.BarChart(res[:policy], Float64.(res[:avg_steps])),
            ylabel = "Average steps")

ptout = Axis(Plots.BarChart(res[:policy], Float64.(res[:t_out])),
             ylabel = "Percentage of time out failures")

g= GroupPlot(1,3, groupStyle="horizontal sep=1cm, vertical sep=1cm");
push!(g, pcoll)
push!(g, pavg)
push!(g, ptout)

save("barplots.tex", g)


text:
module GoGameGraphs

using LinearAlgebra

export Board, GameGraph, export_graph, unique_graphs, smallest_isomorphic_id,
       is_board_graph_connected, remap

using Combinatorics

# Load LightGraphs integration whenever both packages are active.
using Requires
function __init__()
    @require LightGraphs="093fc24a-ae57-5d10-9952-331d41423f4d" include("lightgraphs.jl")
end

# --------- Board code --------------------------------
"""
    Board(id::Integer)

Generate a representation of a go board graph with the given `id`.
The number of nodes is the smallest `n` such that `2^(n*(n-1)/2)>id`.
The edges are found in the following way:
 * Enumerate all potential edges in the order `(1,2), (1,3), ..., (1,n), (2,3), (2,4), ..., (2,n), ..., (n,n)`.
 * Convert `id` to binary with `n*(n-1)/2` digits.
 * The edges in the list matching ones in the binary sequence are included in the graph.
Note: Many of these graphs will be isomorphic. Boards with 2 to 12
nodes are supported (corresponding to `1 <= id < 2^66`).

    Board(edges::Vector{Vector{Int}})

Construct a board graph from a specified list of outgoing edges from
each node. Note, if this does not describe an undirected graph, things
will go bad if it is used to construct a game graph.
"""
mutable struct Board
    edges::Vector{Vector{Int}}
    board::Vector{Int}
end

function num_nodes(id::Integer)
    for n = 1:16
        if Int128(1) << (n * (n - 1) ÷ 2) > id
            return n
        end
    end
    error("id $(id) out of range")
end

function Board(id::Integer)
    id < 1 && error("id must be at least 1.")
    n = num_nodes(id)
    # The code does work with more nodes but the memory requirements
    # start to become cumbersome. Relax this constraint at your own
    # risk.
    n > 12 && error("Board graphs with more than 12 nodes are not supported.")
    @assert(1 <= id < (Int128(1) << (n * (n - 1) ÷ 2)))
    edges = Vector{Int}[Int[] for k = 1:n]
    pairs = collect([i, j] for j = 1:n, i = 1:n if i < j)
    for p in reverse(pairs)
        if (id & 1) == 1
            i, j = p
            push!(edges[i], j)
            push!(edges[j], i)
        end
        id >>= 1
    end
    foreach(reverse!, edges)

    return Board(edges, zeros(Int, length(edges)))
end

Board(edges::Vector{Vector{Int}}) = Board(edges, zeros(Int, length(edges)))

function Base.show(io::IO, board::Board)
    print(io, "Board with $(length(board.edges)) nodes and $(sum(length.(board.edges)) ÷ 2) edges")
end

function clear!(board::Board)
    fill!(board.board, 0)
end

function mark_string_and_find_liberty!(board::Board, node, color)
    board.board[node] |= 4
    for neighbor in board.edges[node]
        if board.board[neighbor] == 0
            return true
        elseif board.board[neighbor] == color
            if mark_string_and_find_liberty!(board, neighbor, color)
                return true
            end
        end
    end

    return false
end

function remove_mark!(board::Board)
    for i = 1:length(board.board)
        board.board[i] &= 3
    end
end

function remove_marked_string!(board::Board)
    for i = 1:length(board.board)
        if board.board[i] & 4 != 0
            board.board[i] = 0
        end
    end
end

# Play a stone on the board. Remove neighbor stones without liberties.
# Remove itself and connected stones if they end up without liberties.
# Return whether the move was a suicide.
function play!(board::Board, node, color)
    board.board[node] = color
    for neighbor in board.edges[node]
        if board.board[neighbor] == 3 - color
            if mark_string_and_find_liberty!(board, neighbor, 3 - color)
                remove_mark!(board)
            else
                remove_marked_string!(board)
            end
        end
    end
    if mark_string_and_find_liberty!(board, node, color)
        remove_mark!(board)
    else
        remove_marked_string!(board)
        return true
    end

    return false
end

# Encode the board in base 3.
ternary_encode(board::Board) = 1 + dot(board.board, 3 .^ (length(board.board)-1:-1:0))

# Set up stones on a board to match a base 3 encoding. This is done by
# playing the stones one by one, starting from the empty board. If and
# only if the position is illegal some stones will be captured or
# suicided in the process, and the obtained position will not match
# the requested base 3 encoding. In that case return false.
function setup_position!(board::Board, pos::Int)
    clear!(board)
    x = string(pos - 1, base = 3, pad = length(board.board))
    for i = 1:length(board.board)
        c = parse(Int, x[i])
        if c > 0
            play!(board, i, c)
        end
    end
    return ternary_encode(board) == pos
end

# Try to make a move on a board and return the base 3 encoding of the
# obtained position. Undo the move before returning.
function trymove(board::Board, j::Int, allow_suicide = true)
    c = 1 + (j - 1) ÷ length(board.board)
    i = mod1(j, length(board.board))
    if board.board[i] != 0
        return 0
    end
    b = copy(board.board)
    is_suicide = play!(board, i, c)
    if !allow_suicide && is_suicide
        pos = 0
    else
        pos = ternary_encode(board)
    end
    board.board = b
    return pos
end

"""
    GameGraph(board::Board)

Compute a representation of the graph of positions for `board`.

    GameGraph(board::Board, allow_suicide = false)

Compute a representation of the graph of positions for `board` with no suicide moves.
"""
mutable struct GameGraph
    edges::Vector{Vector{Int}}
    board::Vector{Vector{Int}}
end

# B is the transition matrix indexed by move number. Illegal moves are
# represented by 0 entries.
# edges splits the B columns into a vector of vectors, and removes the
# 0 entries from each vector.
function GameGraph(board::Board; allow_suicide = true)
    N = length(board.edges)
    B = zeros(Int, 2 * N, 3^N)
    boards = fill(Int[], 3^N)
    for j = 1:3^N
        if setup_position!(board, j)
            boards[j] = copy(board.board)
            for i = 1:2*N
                B[i,j] = trymove(board, i, allow_suicide)
                # Filter out single stone suicides to get rid of
                # self-loops in the graph representation.
                if B[i,j] == j
                    B[i,j] = 0
                end
            end
        end
    end

    B, boards = compress_board_graph(B, boards)
    edges = [sort(vec(filter(x -> x > 0, B[:,k]))) for k = 1:size(B, 2)]

    return GameGraph(edges, boards)
end

function Base.show(io::IO, graph::GameGraph)
    num_nodes = length(graph.edges)
    num_edges = sum(length.(graph.edges))
    print("Game graph with $(num_nodes) nodes and $(num_edges) edges.")
end

# Remove illegal positions. Illegal positions are detected by not
# having any successors in the B matrix since they have never been
# reached. However, on the 1x1 board (and other boards whose graphs
# have no edges) the empty position only have self-suicides and thus
# no successors, breaking the assumption that only illegal positions
# lack successors in B. We handle this by assuming that first column
# in B corresponds to the empty board and never remove it.
function compress_board_graph(B::Matrix{Int}, boards::Vector{Vector{Int}})
    Bsum = vec(sum(B, dims = 1))
    Bsum[1] = 1
    I = findall(Bsum .!= 0)
    J = zeros(Int, size(B, 2))
    J[I] = 1:length(I)
    C = B[:,I]
    C[C.>0] = J[C[C.>0]]
    return C, boards[I]
end

"""
    export_graph(filename, graph::GameGraph)
    export_graph(filename, edges::Vector{Vector{Int}})

Export a graph to file. Format is simplest possible. On line `n` is
listed the outgoing edges from node `n`, separated by spaces.
"""
function export_graph(filename::String, edges::Vector{Vector{Int}})
    f = open(filename, "w")
    for e in edges
        println(f, join(string.(e), " "))
    end
    close(f)
end

function export_graph(filename::String, graph::GameGraph)
    export_graph(filename, graph.edges)
end

"""
    unique_graphs(n)

Identify the set of unique undirected graphs on `N` unlabeled vertices.
"""
function unique_graphs(n)
    n > 8 && error("You probably don't have enough memory to use this function for n > 8.")
    m = n * (n - 1) ÷ 2
    @assert(m <= 62)
    M = 1 << m
    isomorphisms = fill(-1, M)
    a = 1 .<< ((m-1):-1:0)
    for k = 1:M
        if isomorphisms[k] >= 0
            continue
        end
        x = map(x->parse(Int, x),
                split(string(k - 1, base = 2, pad = m), ""))
        mask = .!triu(trues(n, n))
        C = zeros(Int, n, n)
        C[mask] = x
        C += C'
        id = Int[]
        for p in permutations(1:n)
            i = dot(C[p,p][mask], a)
            push!(id, i)
        end
        isomorphisms[id .+ 1] .= minimum(id)
    end
    return unique(isomorphisms)
end

"""
    smallest_isomorphic_id(board::Board)

Find smallest id generating a board graph isomorphic to `board`.
"""
function smallest_isomorphic_id(board::Board)
    N = length(board.edges)
    m = N * (N - 1) ÷ 2
    a = [Int128(1) << x for x in ((m-1):-1:0)]
    C = zeros(Int, N, N)
    for i = 1:N
        for j in board.edges[i]
            C[i,j] = C[j,i] = 1
        end
    end
    mask = .!triu(trues(N, N))
    id = typemax(Int128)
    for p in permutations(1:N)
        i = dot(C[p,p][mask], a)
        id = min(id, i)
    end
    return id
end

"""
    is_board_graph_connected(id, N)

Determine whether the board graph with the given `id` is connected on
`N` nodes.

    is_board_graph_connected(id)

Like above when `N` is the smallest number of nodes that can represent
the board graph `id`, i.e. with no isolated nodes.
"""
function is_board_graph_connected(id, N = num_nodes(id))
    m = N * (N - 1) ÷ 2
    x = map(x->parse(Int, x), split(string(id, base = 2, pad = m), ""))
    mask = .!triu(trues(N, N))
    C = zeros(Int, N, N)
    C[mask] = x
    C += C'
    C += I
    return all(C^(N-1) .> 0)
end

"""
    remap(source_board::Board, target_board::Board)

Compute a remapping of the vertices of the Game Graph generated from
`source_board` to the Game Graph generated from `target_board`.
`source_board` and `target_board` must be isomorphic.

    remap(board::Board)

Like above where `target_board` is the `Board` with smallest id that
is isomorphic to `source_board`.
"""
function remap(source_board::Board,
               target_board::Board = Board(smallest_isomorphic_id(source_board)))
    if length(source_board.edges) != length(target_board.edges)
        error("source_board and target_board are not isomorphic")
    end
    N = length(source_board.edges)
    m = N * (N - 1) ÷ 2
    a = [Int128(1) << x for x in ((m-1):-1:0)]
    C_source = zeros(Int, N, N)
    C_target = zeros(Int, N, N)
    for i = 1:N
        for j in source_board.edges[i]
            C_source[i,j] = C_source[j,i] = 1
        end
        for j in target_board.edges[i]
            C_target[i,j] = C_target[j,i] = 1
        end
    end
    permutation = Int[]
    for p in permutations(1:N)
        if C_source[p, p] == C_target
            permutation = p
            break
        end
    end
    if isempty(permutation)
        error("source_board and target_board are not isomorphic")
    end

    B = Int[]
    for i = 1:3^N
        if setup_position!(source_board, i)
            j = 1 + parse(Int,
                          string(i - 1, base = 3, pad = N)[permutation],
                          base = 3)
            push!(B, j)
        end
    end

    BB = sort(B)
    M = Int[]
    for i = 1:length(B)
        push!(M, findall(BB .== B[i])[1])
    end

    return M
end

end

max_stars_repo_name: GunnarFarneback/GoGameGraphs

max_issues_repo_name: GunnarFarneback/GoGameGraphs

max_forks_repo_name: GunnarFarneback/GoGameGraphs


text:
include("population.jl")

function coverage(bestFits)
    len = length(bestFits)
    if bestFits[len] == 0
        return true
    elseif len < 3
        return false
    else
        return bestFits[len-1] == bestFits[len]
    end
end

function geneticAlgorithm(data, crossoverPoint1,crossoverPoint2,values,weights,maxCapacity, mutationPercentage,elitePercentage)
    calculatePopulationFitness!(data,values,weights,maxCapacity)
    populationSize = length(data)
    eliteSize = Int(trunc(elitePercentage*populationSize))
    eliteSize = eliteSize + (populationSize - eliteSize)%2
    bestFits = [data[1].fitness]
    while !coverage(bestFits)
        elite = deepcopy(selectPopulation(data,eliteSize))
        data = selectPopulation(data, populationSize-eliteSize)
        data = crossoverPopulation(data,crossoverPoint1,crossoverPoint2)
        mutatePopulation!(data,mutationPercentage)
        data = [data;elite]
        calculatePopulationFitness!(data,values,weights,maxCapacity)
        bestFits = [bestFits;data[1].fitness]
    end
    return length(bestFits), data
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
module StrayCopses

using StatsBase
using Random

include("impurity.jl")
include("decision_tree.jl")
include("copse.jl")

export predict, fit, StrayCopse, create_tree

end # module


text:
const SIDE_COEF_SQUARE = 1 / 2
const SIDE_COEF_CIRCLE = 1 / sqrt(pi)
const SIDE_COEF_DIAMOND = sqrt(2) / 2
const SIDE_COEF_TRIANGLE = sqrt(2 / (3 * sqrt(3 / 2)))
const SIDE_COEF_PENTAGON = sqrt(2 / (5 * sin(1.2566)))
const SIDE_COEF_HEXAGON = sqrt(2 / (6 * sin(1.0472)))
const SIDE_COEF_OCTAGON = sqrt(2 / (8 * sin(0.7854)))

function _polygon_radius(a::Real, n::Integer, r::Real, start::Real=-pi / 2)
    a -= start
    while true
        if a > 2pi / n
            a -= 2pi / n
        elseif a < 0
            a += 2pi / n
        else
            break
        end
    end
    return r * cos(pi / n) / cos(pi / n - a)
end

"""
    node_radius(::Type{Val{:shape}}, s::Real, a::Real)

Get the distance between the center and the edge of a node at angle `a`.
The node is of shape `:shape` and size `s`.
"""
node_radius(_, s::Real, a::Real) = s / 2

node_radius(::Type{Val{:square}}, s::Real, a::Real) =
    s * SIDE_COEF_SQUARE / max(abs(sin(a)), abs(cos(a)))

node_radius(::Type{Val{:circle}}, s::Real, a::Real) = s * SIDE_COEF_CIRCLE

node_radius(::Type{Val{:diamond}}, s::Real, a::Real) =
    s * SIDE_COEF_DIAMOND / (1 + abs(sin(a) / cos(a))) / abs(cos(a))

node_radius(::Type{Val{:triangle}}, s::Real, a::Real) =
    _polygon_radius(a, 3, s * SIDE_COEF_TRIANGLE)

node_radius(::Type{Val{:pentagon}}, s::Real, a::Real) =
    _polygon_radius(a, 5, s * SIDE_COEF_PENTAGON)

node_radius(::Type{Val{:hexagon}}, s::Real, a::Real) =
    _polygon_radius(a, 6, s * SIDE_COEF_HEXAGON)

node_radius(::Type{Val{:octagon}}, s::Real, a::Real) =
    _polygon_radius(a, 8, s * SIDE_COEF_OCTAGON, pi / 8)

function _polygon_path(context::CairoContext, x, y, n, r, start=-pi / 2)
    a = start
    move_to(context, x + r * cos(a), y + r * sin(a))
    for i = 1:n - 1
        a += 2pi / n
        line_to(context, x + r * cos(a), y + r * sin(a))
    end
    close_path(context)
end

"""
    outline_node!(::Type{Val{:shape}}, context, x, y, s)

Create a path on a Cairo `context` that corresponds to `:shape` of size `s`.
"""
outline_node!(_, context, x, y, r) = nothing

outline_node!(::Type{Val{:circle}}, context, x, y, r) = circle(context, x, y, r * SIDE_COEF_CIRCLE)

outline_node!(::Type{Val{:square}}, context, x, y, r) = rectangle(context,
    x - r * SIDE_COEF_SQUARE, y - r * SIDE_COEF_SQUARE, r, r)

outline_node!(::Type{Val{:diamond}}, context, x, y, side) = begin
    move_to(context, x - side * SIDE_COEF_DIAMOND, y)
    line_to(context, x, y - side * SIDE_COEF_DIAMOND)
    line_to(context, x + side * SIDE_COEF_DIAMOND, y)
    line_to(context, x, y + side * SIDE_COEF_DIAMOND)
    close_path(context)
end

outline_node!(::Type{Val{:triangle}}, context, x, y, side) =
    _polygon_path(context, x, y, 3, side * SIDE_COEF_TRIANGLE)

outline_node!(::Type{Val{:pentagon}}, context, x, y, side) =
    _polygon_path(context, x, y, 5, side * SIDE_COEF_PENTAGON)

outline_node!(::Type{Val{:hexagon}}, context, x, y, side) =
    _polygon_path(context, x, y, 6, side * SIDE_COEF_HEXAGON)

outline_node!(::Type{Val{:octagon}}, context, x, y, side) =
    _polygon_path(context, x, y, 8, side * SIDE_COEF_OCTAGON, pi / 8)


text:
using SequenceInference
using Test

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

max_stars_repo_name: guchchcug/SequenceInference

max_issues_repo_name: guchchcug/SequenceInference

max_forks_repo_name: guchchcug/SequenceInference


text:
export test_jacobian, test_gradients, linedata, lineview, meshdata,
    meshview, gradview, jacview, PCLview, pcolormeshview,
    animate, saveanim, test_hessian, require_optim

function require_pyplot()
    if !isdefined(Main, :PyPlot)
        error("You must load PyPlot to use this function, e.g., `using PyPlot` or `import PyPlot`")
    end
    Main.PyPlot
end

function require_optim()
    if !isdefined(Main, :Optim)
        error("You must load Optim to use this function, e.g., `import Optim` or `using Optim`")
    end
    Main.Optim
end

"""
    test_gradients(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0, showfig::Bool = true)

Testing the gradients of a vector function `f`:
`y, J = f(x)` where `y` is a scalar output and `J` is the vector gradient.
"""
function test_gradients(f::Function, x0::Array{Float64, 1};
    scale::Float64 = 1.0, showfig::Bool = true, mpi::Bool = false)

    v0 = rand(Float64,length(x0))
    γs = scale ./10 .^(1:5)
    err2 = Float64[]
    err1 = Float64[]
    f0, J = f(x0)
    for i = 1:5
        f1, _ = f(x0+γs[i]*v0)
        push!(err1, abs(f1-f0))
        push!(err2, abs(f1-f0-γs[i]*sum(J.*v0)))
    end
    if showfig
        if mpi && mpi_rank()>0
            return
        end
        mp = require_pyplot()
        mp.close("all")
        mp.loglog(γs, err1, label="Finite Difference")
        mp.loglog(γs, err2, label="Automatic Differentiation")
        mp.loglog(γs, γs * 0.5*abs(err1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
        mp.loglog(γs, γs.^2 * 0.5*abs(err2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
        mp.plt.gca().invert_xaxis()
        mp.legend()
        mp.savefig("test.png")
        @info "Results saved to test.png"
        println("Finite difference: $err1")
        println("Automatic differentiation: $err2")
    end
    return err1, err2
end

"""
    test_jacobian(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0, showfig::Bool = true)

Testing the gradients of a vector function `f`:
`y, J = f(x)` where `y` is a vector output and `J` is the Jacobian.
"""
function test_jacobian(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0, showfig::Bool = true)

    v0 = rand(Float64,size(x0))
    γs = scale ./10 .^(1:5)
    err2 = Float64[]
    err1 = Float64[]
    f0, J = f(x0)
    for i = 1:5
        f1, _ = f(x0+γs[i]*v0)
        push!(err1, norm(f1-f0))
        push!(err2, norm(f1-f0-γs[i]*J*v0))
    end
    if showfig
        mp = require_pyplot()
        mp.close("all")
        mp.loglog(γs, err1, label="Finite Difference")
        mp.loglog(γs, err2, label="Automatic Differentiation")
        mp.loglog(γs, γs * 0.5*abs(err1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
        mp.loglog(γs, γs.^2 * 0.5*abs(err2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
        mp.plt.gca().invert_xaxis()
        mp.legend()
        mp.savefig("test.png")
        @info "Results saved to test.png"
        println("Finite difference: $err1")
        println("Automatic differentiation: $err2")
    end
    return err1, err2
end

"""
    test_hessian(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0)

Testing the Hessian of a scalar function `f`:
`g, H = f(x)` where `y` is a scalar output, `g` is a vector gradient output, and `H` is the Hessian.
"""
test_hessian(f::Function, x0::Array{Float64, 1}; kwargs...) = test_jacobian(f, x0;  kwargs...)

function linedata(θ1, θ2=nothing; n::Integer = 10)
    if θ2 === nothing
        θ2 = θ1 .* (1 .+ randn(size(θ1)...))
    end
    T = []
    for x in LinRange(0.0,1.0,n)
        push!(T, θ1*(1-x)+θ2*x)
    end
    T
end

function lineview(losses::Array{Float64})
    mp = require_pyplot()
    n = length(losses)
    v = collect(LinRange(0.0,1.0,n))
    mp.close("all")
    mp.plot(v, losses)
    mp.xlabel("\$\\alpha\$")
    mp.ylabel("loss")
    mp.grid("on")
end

@doc raw"""
    lineview(sess::PyObject, pl::PyObject, loss::PyObject, θ1, θ2=nothing; n::Integer = 10)

Plots the function

$$h(α) = f((1-α)θ_1 + αθ_2)$$

# Example
```julia
pl = placeholder(Float64, shape=[2])
l = sum(pl^2-pl*0.1)
sess = Session(); init(sess)
lineview(sess, pl, l, rand(2))
```
"""
function lineview(sess::PyObject, pl::PyObject, loss::PyObject, θ1, θ2=nothing; n::Integer = 10)
    mp = require_pyplot()
    dat = linedata(θ1, θ2, n=n)
    V = zeros(length(dat))
    for i = 1:length(dat)
        @info i, n
        V[i] = run(sess, loss, pl=>dat[i])
    end
    lineview(V)
end

function meshdata(θ;
     a::Real=1, b::Real=1, m::Integer=10, n::Integer=10 ,direction::Union{Array{<:Real}, Missing}=missing)
    local δ, γ
    as = LinRange(-a, a, m)
    bs = LinRange(-b, b, n)
    α = zeros(m, n)
    β = zeros(m, n)
    θs = Array{Any}(undef, m, n)
    if !ismissing(direction)
        δ = direction - θ
        γ = randn(size(θ)...)
        γ = γ/norm(γ)*norm(δ)
    else
        δ = randn(size(θ)...)
        γ = randn(size(θ)...)
    end
    for i = 1:m
        for j = 1:n
            α[i,j] = as[i]
            β[i,j] = bs[j]
            θs[i,j] = θ + δ*as[i] + γ*bs[j]
        end
    end
    return θs
end

function meshview(losses::Array{Float64}, a::Real=1, b::Real=1)
    mp = require_pyplot()
    m, n = size(losses)
    α = zeros(m, n)
    β = zeros(m, n)
    as = LinRange(-a, a, m)
    bs = LinRange(-b, b, n)
    for i = 1:m
        for j = 1:n
            α[i,j] = as[i]
            β[i,j] = bs[j]
        end
    end
    mp.close("all")
    mp.mesh(α, β, losses)
    mp.xlabel("\$\\alpha\$")
    mp.ylabel("\$\\beta\$")
    mp.zlabel("loss")
    mp.scatter3D(α[(m+1)÷2, (n+1)÷2], β[(m+1)÷2, (n+1)÷2], losses[(m+1)÷2, (n+1)÷2], color="red", s=100)
    return α, β, losses
end

function pcolormeshview(losses::Array{Float64}, a::Real=1, b::Real=1)
    mp = require_pyplot()
    m, n = size(losses)
    α = zeros(m, n)
    β = zeros(m, n)
    as = LinRange(-a, a, m)
    bs = LinRange(-b, b, n)
    for i = 1:m
        for j = 1:n
            α[i,j] = as[i]
            β[i,j] = bs[j]
        end
    end
    mp.close("all")
    mp.pcolormesh(α, β, losses)
    mp.xlabel("\$\\alpha\$")
    mp.ylabel("\$\\beta\$")
    mp.scatter(0.0,0.0, marker="*", s=100)
    mp.colorbar()
    return α, β, losses
end

function meshview(sess::PyObject, pl::PyObject, loss::PyObject, θ;
        a::Real=1, b::Real=1, m::Integer=9, n::Integer=9,
        direction::Union{Array{<:Real}, Missing}=missing)
    dat = meshdata(θ; a=a, b=b, m=m, n=n, direction=direction)
    m, n = size(dat)
    V = zeros(m, n)
    for i = 1:m
        @info i, m
        for j = 1:n
            V[i,j] = run(sess, loss, pl=>dat[i,j])
        end
    end
    meshview(V, a, b)
end

function pcolormeshview(sess::PyObject, pl::PyObject, loss::PyObject, θ;
    a::Real=1, b::Real=1, m::Integer=9, n::Integer=9,
    direction::Union{Array{<:Real}, Missing}=missing)
    dat = meshdata(θ; a=a, b=b, m=m, n=n, direction=direction)
    m, n = size(dat)
    V = zeros(m, n)
    for i = 1:m
        @info i, m
        for j = 1:n
            V[i,j] = run(sess, loss, pl=>dat[i,j])
        end
    end
    pcolormeshview(V, a, b)
end

function gradview(sess::PyObject, pl::PyObject, loss::PyObject, u0, grad::PyObject;
    scale::Float64 = 1.0, mpi::Bool = false)
    mp = require_pyplot()
    local v
    if length(size(u0))==0
        v = rand()
    else
        v = rand(length(u0))
    end
    γs = scale ./ 10 .^ (1:5)
    v1 = Float64[]
    v2 = Float64[]
    L_ = run(sess, loss, pl=>u0)
    J_ = run(sess, grad, pl=>u0)
    for i = 1:5
        @info i
        L__ = run(sess, loss, pl=>u0+v*γs[i])
        push!(v1, norm(L__-L_))
        if size(J_)==size(v)
            if length(size(J_))==0
                push!(v2, norm(L__-L_-γs[i]*sum(J_*v)))
            else
                push!(v2, norm(L__-L_-γs[i]*sum(J_.*v)))
            end
        else
            push!(v2, norm(L__-L_-γs[i]*J_*v))
        end

    end
    if !(mpi) || (mpi && mpi_rank()==0)
        mp.close("all")
        mp.loglog(γs, abs.(v1), "*-", label="finite difference")
        mp.loglog(γs, abs.(v2), "+-", label="automatic linearization")
        mp.loglog(γs, γs.^2 * 0.5*abs(v2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
        mp.loglog(γs, γs * 0.5*abs(v1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
        mp.plt.gca().invert_xaxis()
        mp.legend()
        mp.xlabel("\$\\gamma\$")
        mp.ylabel("Error")
        if mpi
            mp.savefig("mpi_gradview.png")
        end
    end

    return v1, v2
end

"""
    gradview(sess::PyObject, pl::PyObject, loss::PyObject, u0; scale::Float64 = 1.0)

Visualizes the automatic differentiation and finite difference convergence converge. For correctly implemented
differentiable codes, the convergence rate for AD should be 2 and for FD should be 1 (if not evaluated at stationary point).

- `scale`: you can control the step size for perturbation.
"""
function gradview(sess::PyObject, pl::PyObject, loss::PyObject, u0; scale::Float64 = 1.0, mpi::Bool = false)
    grad = tf.convert_to_tensor(gradients(loss, pl))
    gradview(sess, pl, loss, u0, grad, scale = scale, mpi = mpi)
end

@doc raw"""
    jacview(sess::PyObject, f::Function, θ::Union{Array{Float64}, PyObject, Missing},
    u0::Array{Float64}, args...)

Performs gradient test for a vector function. `f` has the signature
```
f(θ, u) -> r, J
```
Here `θ` is a nuisance  parameter, `u` is the state variables (w.r.t. which the Jacobian is computed),
`r` is the residual vector, and `J` is the Jacobian matrix (a dense matrix or a [`SparseTensor`](@ref)).

# Example 1
```julia
u0 = rand(10)
function verify_jacobian_f(θ, u)
    r = u^3+u - u0
    r, spdiag(3u^2+1.0)
end
jacview(sess, verify_jacobian_f, missing, u0)
```

# Example 2
```
u0 = rand(10)
rs = rand(10)
function verify_jacobian_f(θ, u)
    r = [u^2;u] - [rs;rs]
    r, [spdiag(2*u); spdiag(10)]
end
jacview(sess, verify_jacobian_f, missing, u0); close("all")
```
"""
function jacview(sess::PyObject, f::Function, θ::Union{Array{Float64}, PyObject, Missing},
                u0::Array{Float64}, args...)
    mp = require_pyplot()
    u = placeholder(Float64, shape=[length(u0)])
    L, J = f(θ, u)
    init(sess)
    L_ = run(sess, L, u=>u0, args...)
    J_ = run(sess, J, u=>u0, args...)
    v = rand(length(u0))
    γs = 1.0 ./ 10 .^ (1:5)
    v1 = Float64[]
    v2 = Float64[]
    for i = 1:5
        L__ = run(sess, L, u=>u0+v*γs[i], args...)
        push!(v1, norm(L__-L_))
        push!(v2, norm(L__-L_-γs[i]*J_*v))
    end
    mp.close("all")
    mp.loglog(γs, abs.(v1), "*-", label="finite difference")
    mp.loglog(γs, abs.(v2), "+-", label="automatic linearization")
    mp.loglog(γs, γs.^2 * 0.5*abs(v2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    mp.loglog(γs, γs * 0.5*abs(v1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
    mp.plt.gca().invert_xaxis()
    mp.legend()
    mp.xlabel("\$\\gamma\$")
    mp.ylabel("Error")
end

function PCLview(sess::PyObject, f::Function, L::Function, θ::Union{PyObject,Array{Float64,1}, Float64},
    u0::Union{PyObject, Array{Float64}}, args...; options::Union{Dict{String, T}, Missing}=missing) where T<:Real
    mp = require_pyplot()
    if isa(θ, PyObject)
        θ = run(sess, θ, args...)
    end
    x = placeholder(Float64, shape=[length(θ)])
    l, u, g = NonlinearConstrainedProblem(f, L, x, u0; options=options)
    init(sess)
    L_ = run(sess, l, x=>θ, args...)
    J_ = run(sess, g, x=>θ, args...)
    v = rand(length(x))
    γs = 1.0 ./ 10 .^ (1:5)
    v1 = Float64[]
    v2 = Float64[]
    for i = 1:5
        @info i
        L__ = run(sess, l, x=>θ+v*γs[i], args...)
        # @show L__,L_,J_, v
        push!(v1, L__-L_)
        push!(v2, L__-L_-γs[i]*sum(J_.*v))
    end
    mp.close("all")
    mp.loglog(γs, abs.(v1), "*-", label="finite difference")
    mp.loglog(γs, abs.(v2), "+-", label="automatic linearization")
    mp.loglog(γs, γs.^2 * 0.5*abs(v2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    mp.loglog(γs, γs * 0.5*abs(v1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
    mp.plt.gca().invert_xaxis()
    mp.legend()
    mp.xlabel("\$\\gamma\$")
    mp.ylabel("Error")
end

@doc raw"""
    animate(update::Function, frames; kwargs...)

Creates an animation using update function `update`.
# Example
```julia
θ = LinRange(0, 2π, 100)
x = cos.(θ)
y = sin.(θ)
pl, = plot([], [], "o-")
t = title("0")
xlim(-1.2,1.2)
ylim(-1.2,1.2)
function update(i)
    t.set_text("$i")
    pl.set_data([x[1:i] y[1:i]]'|>Array)
end
animate(update, 1:100)
```

!!! info
    `clf()` can be used to clear figures within `update` function.
"""
function animate(update::Function, frames; kwargs...)
    mp = require_pyplot()
    animation_ = pyimport("matplotlib.animation")
    if !isa(frames, Integer)
        frames = Array(frames)
    end
    animation_.FuncAnimation(mp.gcf(), update, frames=frames; kwargs...)
end

"""
    saveanim(anim::PyObject, filename::String; kwargs...)

Saves the animation produced by [`animate`](@ref)
"""
function saveanim(anim::PyObject, filename::String; kwargs...)
    if Sys.iswindows()
        anim.save(filename, "ffmpeg"; kwargs...)
    else
        anim.save(filename, "imagemagick"; kwargs...)
    end
end


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/abc_run.jl
## (c) 2018 Danley C. Hsu & Eric B. Ford
# Script for producing Q1-Q16 FGK planet candidate occurrence rate estimates

using ExoplanetsSysSim
include(joinpath(pwd(), "abc_setup.jl"))

using .SysSimABC
using ExoplanetsSysSim
using StatsBase
using Profile

out2txt = false # Write occurrence rates & densities to text files
expandpart = true # Expand final generation for robust posteriors
profile_code = false

global output
if profile_code
   println("Setting up simulation...")
   @time abc_plan = setup_abc(0,max_generations=1)
   println("Running simulation...")
   @time output = run_abc(abc_plan)
   println("Setting up simulation...")
   @time abc_plan = setup_abc(0,max_generations=5)
   Profile.init(n = 10^7, delay = 0.01)
   Profile.clear_malloc_data()
   println("Running simulation...")
   @profile output = run_abc(abc_plan)
   open("profile.txt.corbits.reallocated", "w") do s
       Profile.print(IOContext(s, :displaysize => (24, 500)))
   end
else
    println("Setting up simulation...")
    @time abc_plan = setup_abc(0,max_generations=50)
    println("")
    println("Running simulation...")
    @time output = run_abc(abc_plan)
    println("")
end

if expandpart
    println("Expanding to large generation...")
    @time theta_largegen, weights_largegen = run_abc_largegen(output, EvalSysSimModel.get_ss_obs(), output.accept_log.epsilon[end])
    println("")
end

if out2txt
    file_rate = open("rate_output.txt", "w")
    file_dens = open("dens_output.txt", "w")
end

limitP = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})::Array{Float64,1}
limitR = get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1})::Array{Float64,1}
const r_dim = length(limitR)-1

if expandpart
    weight_vec = pweights(weights_largegen)
else
    weight_vec = pweights(output.weights)
end

for p_ind = 1:(length(limitP)-1)
    col_ind = (p_ind-1)*(r_dim+1)+1
    for r_ind = 1:(length(limitR)-1)
        dens_denom = 1.0/(log2(limitP[p_ind+1])-log2(limitP[p_ind]))/(log2(limitR[r_ind+1])-log2(limitR[r_ind]))

        bin_ind = (p_ind-1)*(r_dim+1)+1+r_ind
        if expandpart
            quant_arr = quantile(theta_largegen[bin_ind,:].*theta_largegen[col_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
        else
            quant_arr = quantile(output.theta[bin_ind,:].*output.theta[col_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
        end

        println("-----------------------------")
        println("Orbital Period (day) = ", string(limitP[p_ind:p_ind+1]), " / Planet Radius (R_earth) = ", string(limitR[r_ind:r_ind+1]/ExoplanetsSysSim.earth_radius))
        println("")
        println("Rate = ", string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1]))
        println("Density = ", string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom))

        if out2txt
            write(file_rate, string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1], "\n"))
            write(file_dens, string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom, "\n"))
        end
    end
end

if out2txt
    close(file_rate)
    close(file_dens)
end

println("-----------------------------")
println("")
println(EvalSysSimModel.get_ss_obs())
if expandpart
    save(string("test-pop-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs(), "theta_largegen", theta_largegen, "weights_largegen", weights_largegen)
else
    save(string("test-pop-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs())
end


text:
using Test

using MLJBase
using LinearAlgebra
using Statistics
using MLJGLMInterface
import GLM

import Distributions
import StableRNGs
using Tables

expit(X) = 1 ./ (1 .+ exp.(-X))

###
### OLSREGRESSOR
###

X, y = @load_boston

train, test = partition(eachindex(y), 0.7)

atom_ols = LinearRegressor()

Xtrain = selectrows(X, train)
ytrain = selectrows(y, train)
Xtest  = selectrows(X, test)

fitresult, _, report = fit(atom_ols, 1, Xtrain, ytrain)
θ = MLJBase.fitted_params(atom_ols, fitresult)

p = predict_mean(atom_ols, fitresult, Xtest)

# hand made regression to compare

Xa = MLJBase.matrix(X) # convert(Matrix{Float64}, X)
Xa1 = hcat(Xa, ones(size(Xa, 1)))
coefs = Xa1[train, :] \ y[train]
p2 = Xa1[test, :] * coefs

@test p ≈ p2

model = atom_ols
@test name(model) == "LinearRegressor"
@test package_name(model) == "GLM"
@test is_pure_julia(model)
@test is_supervised(model)
@test package_license(model) == "MIT"
@test prediction_type(model) == :probabilistic
@test hyperparameters(model) == (:fit_intercept, :allowrankdeficient, :offsetcol)
@test hyperparameter_types(model) == ("Bool", "Bool", "Union{Nothing, Symbol}")

p_distr = predict(atom_ols, fitresult, selectrows(X, test))

@test p_distr[1] == Distributions.Normal(p[1], GLM.dispersion(fitresult))

###
### Logistic regression
###

rng = StableRNGs.StableRNG(0)

N = 100
X = MLJBase.table(rand(rng, N, 4));
ycont = 2*X.x1 - X.x3 + 0.6*rand(rng, N)
y = (ycont .> mean(ycont)) |> categorical;

lr = LinearBinaryClassifier()
fitresult, _, report = fit(lr, 1, X, y)

yhat = predict(lr, fitresult, X)
@test mean(cross_entropy(yhat, y)) < 0.25

pr = LinearBinaryClassifier(link=GLM.ProbitLink())
fitresult, _, report = fit(pr, 1, X, y)
yhat = predict(lr, fitresult, X)
@test mean(cross_entropy(yhat, y)) < 0.25

fitted_params(pr, fitresult)

###
### Count regression
###

rng = StableRNGs.StableRNG(123)

X = randn(rng, 500, 5)
θ = randn(rng, 5)
y = map(exp.(X*θ)) do mu
    rand(rng, Distributions.Poisson(mu))
end

XTable = MLJBase.table(X)

lcr = LinearCountRegressor(fit_intercept=false)
fitresult, _, _ = fit(lcr, 1, XTable, y)

θ̂ = fitted_params(lcr, fitresult).coef

@test norm(θ̂ .- θ)/norm(θ) ≤ 0.03

modeltypes = [LinearRegressor, LinearBinaryClassifier, LinearCountRegressor]
@testset "Test prepare_inputs" begin
    @testset "intercept/offsetcol" for mt in modeltypes
            X = (x1=[1,2,3], x2=[4,5,6])
            m = mt(fit_intercept=true, offsetcol=:x2)
            Xmatrix, offset = MLJGLMInterface.prepare_inputs(m, X)

            @test offset == [4, 5, 6]
            @test Xmatrix== [1 1;
                             2 1;
                             3 1]
    end

    @testset "no intercept/no offsetcol" for mt in modeltypes
        X = (x1=[1,2,3], x2=[4,5,6])
        m = mt(fit_intercept=false)
        Xmatrix, offset = MLJGLMInterface.prepare_inputs(m, X)

        @test offset == []
        @test Xmatrix == [1 4;
                          2 5;
                          3 6]
    end

end

@testset "Test offsetting models" begin
    @testset "Test split_X_offset" begin
        X = (x1=[1,2,3], x2=[4,5,6])
        @test MLJGLMInterface.split_X_offset(X, nothing) == (X, Float64[])
        @test MLJGLMInterface.split_X_offset(X, :x1) == ((x2=[4,5,6],), [1,2,3])

        X = MLJBase.table(rand(rng, N, 3))
        Xnew, offset = MLJGLMInterface.split_X_offset(X, :x2)
        @test offset isa Vector
        @test length(Xnew) == 2
    end

    # In the following:
    # The second column is taken as an offset by the model
    # This is equivalent to assuming the coef is 1 and known

    @testset "Test Logistic regression with offset" begin
        N = 1000
        rng = StableRNGs.StableRNG(0)
        X = MLJBase.table(rand(rng, N, 3))
        y = rand(rng, Distributions.Uniform(0,1), N) .< expit(2*X.x1 + X.x2 - X.x3)
        y = categorical(y)

        lr = LinearBinaryClassifier(fit_intercept=false, offsetcol=:x2)
        fitresult, _, report = fit(lr, 1, X, y)
        fp = fitted_params(lr, fitresult)

        @test fp.coef ≈ [2, -1] atol=0.03
    end
    @testset "Test Linear regression with offset" begin
        N = 1000
        rng = StableRNGs.StableRNG(0)
        X = MLJBase.table(rand(rng, N, 3))
        y = 2*X.x1 + X.x2 - X.x3 + rand(rng, Distributions.Normal(0,1), N)

        lr = LinearRegressor(fit_intercept=false, offsetcol=:x2)
        fitresult, _, report = fit(lr, 1, X, y)
        fp = fitted_params(lr, fitresult)

        @test fp.coef ≈ [2, -1] atol=0.07
    end
    @testset "Test Count regression with offset" begin
        N = 1000
        rng = StableRNGs.StableRNG(0)
        X = MLJBase.table(rand(rng, N, 3))
        y = map(exp.(2*X.x1 + X.x2 - X.x3)) do mu
            rand(rng, Distributions.Poisson(mu))
        end

        lcr = LinearCountRegressor(fit_intercept=false, offsetcol=:x2)
        fitresult, _, _ = fit(lcr, 1, X, y)
        fp = fitted_params(lcr, fitresult)

        @test fp.coef ≈ [2, -1] atol=0.04
    end
end


text:
module Outputs

using QuantumESPRESSOParser.Outputs: SubroutineError

export SubroutineError

module PWscf

using QuantumESPRESSOParser.Outputs.PWscf:
    Diagonalization,
    Preamble,
    Diagonalization,
    Davidson,
    ConjugateGradient,
    ProjectedPreconditionedConjugateGradient,
    parse_fft_base_info,
    parse_symmetries,
    parse_ibz,
    parse_stress,
    parse_iteration_time,
    parse_bands,
    parse_all_electron_energy,
    parse_energy_decomposition,
    parse_paw_contribution,
    parse_smearing_energy,
    parse_version,
    parse_parallel_info,
    parse_fft_dimensions,
    parse_iteration_head,
    parse_electrons_energies,
    parse_clock,
    # whatinput,
    isoptimized,
    isjobdone,
    tryparsefirst,
    parsefirst,
    tryparseall,
    parseall,
    tryparselast,
    parselast,
    tryparsenext,
    parsenext,
    tryparsefinal,
    parsefinal

export Diagonalization,
    Preamble,
    Davidson,
    ConjugateGradient,
    ProjectedPreconditionedConjugateGradient,
    parse_fft_base_info,
    parse_symmetries,
    parse_ibz,
    parse_stress,
    parse_iteration_time,
    parse_bands,
    parse_all_electron_energy,
    parse_energy_decomposition,
    parse_paw_contribution,
    parse_smearing_energy,
    parse_version,
    parse_parallel_info,
    parse_fft_dimensions,
    parse_iteration_head,
    parse_electrons_energies,
    parse_clock,
    # whatinput,
    isoptimized,
    isjobdone,
    tryparsefirst,
    parsefirst,
    tryparseall,
    parseall,
    tryparselast,
    parselast,
    tryparsenext,
    parsenext,
    tryparsefinal,
    parsefinal

end # module PWscf

module PHonon

using QuantumESPRESSOParser.Outputs.PHonon: parse_frequency, parse_dos

export parse_frequency, parse_dos

end # module PHonon

end # module Outputs


text:
# lstm weights initialization
# w[2k-1], w[2k] : weight and bias for kth layer respectively
function initweights(atype, hiddens, embedding, vocab, init=xavier)
    weights = Array(Any, 2length(hiddens))
    input = embedding
    for k = 1:length(hiddens)
        weights[2k-1] = init(input+hiddens[k], 4hiddens[k])
        weights[2k] = zeros(1, 4hiddens[k])
        weights[2k][1:hiddens[k]] = 1 # forget gate bias
        input = hiddens[k]
    end
    return map(w->convert(atype, w), weights)
end

# state initialization
# s[2k-1], s[2k] : hidden and cell respectively
function initstate(atype, hiddens, batchsize)
    state = Array(Any, 2length(hiddens))
    for k=1:length(hiddens)
        state[2k-1] = atype(zeros(batchsize, hiddens[k]))
        state[2k] = atype(zeros(batchsize, hiddens[k]))
    end
    return state
end

function initmodel(atype, hiddens, embedding, vocabsize, init=xavier)
    model = Dict{Symbol, Any}()
    model[:forw] = initweights(atype, hiddens, embedding, vocabsize, init)
    model[:back] = initweights(atype, hiddens, embedding, vocabsize, init)
    model[:fembed] = atype(init(vocabsize, embedding))
    model[:bembed] = atype(init(vocabsize, embedding))
    model[:soft] = [ atype(init(2hiddens[end], vocabsize)), atype(init(1, vocabsize)) ]
    return model
end

function lstm(weight,bias,hidden,cell,input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

# multilayer lstm forward, returns the final hidden
function forward(weight, states, input)
    x = input
    for i=1:2:length(states)
        (states[i], states[i+1]) = lstm(weight[i], weight[i+1], states[i], states[i+1], x)
        x = states[i]
    end
    return x
end

# dropout layer
function dropout(x, p)
    if p > 0
        return x .* (rand!(similar(AutoGrad.getval(x))) .> p) * (1/(1-p))
    else
        return x
    end
end

function logprob(output, ypred)
    nrows,ncols = size(ypred)
    index = similar(output)
    @inbounds for i=1:length(output)
        index[i] = i + (output[i]-1)*nrows
    end
    o1 = logp(ypred,2)
    o2 = o1[index]
    o3 = sum(o2)
    return o3
end

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdwordcloud

"""
    ''_antdwordcloud(;kwargs...)

An AntdWordCloud component.

Keyword arguments:
- `id` (String; optional)
- `annotations` (optional)
- `appendPadding` (Real | Array of Reals; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `color` (optional): . color has the following type: String | Array of Strings | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `colorField` (String; optional)
- `data` (Array of Dicts; optional)
- `height` (Real; optional)
- `imageMask` (String; optional)
- `label` (optional)
- `legend` (optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-CN', 'en-US'; optional)
- `meta` (optional)
- `padding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `placementStrategy` (optional): . placementStrategy has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `recentlyWordClickRecord` (optional): . recentlyWordClickRecord has the following type: lists containing elements 'timestamp', 'data'.
Those elements have the following types:
  - `timestamp` (Real; optional)
  - `data` (Dict; optional)
- `renderer` (a value equal to: 'canvas', 'svg'; optional)
- `spiral` (a value equal to: 'archimedean', 'rectangular'; optional)
- `style` (Dict; optional)
- `tooltip` (optional)
- `weightField` (String; required)
- `width` (Real; optional)
- `wordField` (String; required)
- `wordStyle` (optional): . wordStyle has the following type: lists containing elements 'fontFamily', 'fontWeight', 'padding', 'fontSize'.
Those elements have the following types:
  - `fontFamily` (String; optional)
  - `fontWeight` (String; optional)
  - `padding` (Real; optional)
  - `fontSize` (Array of Reals; optional)
"""
function ''_antdwordcloud(; kwargs...)
        available_props = Symbol[:id, :annotations, :appendPadding, :autoFit, :className, :color, :colorField, :data, :height, :imageMask, :label, :legend, :loading_state, :locale, :meta, :padding, :placementStrategy, :recentlyWordClickRecord, :renderer, :spiral, :style, :tooltip, :weightField, :width, :wordField, :wordStyle]
        wild_props = Symbol[]
        return Component("''_antdwordcloud", "AntdWordCloud", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
# Generate a toy problem without any obvious structure in the mean, precision, or noise std.
# Important to ensure that the unit tests don't just pass for a special case by accident.
# Everything should be reasonably well conditioned.
function generate_toy_problem(rng, N, D, ::Type{<:Matrix})
    X = randn(rng, D, N)
    B = randn(rng, D, D)
    C = 0.1 * randn(rng, N, N)
    mw, Λw, Σy = randn(rng, D), B * B' + I, C * C' + I
    return X, BayesianLinearRegressor(mw, Λw), Σy
end

function generate_toy_problem(rng, N, D, ::Type{<:ColVecs})
    X, f, Σy = generate_toy_problem(rng, N, D, Matrix)
    return ColVecs(X), f, Σy
end

function generate_toy_problem(rng, N, D, ::Type{<:RowVecs})
    X, f, Σy = generate_toy_problem(rng, N, D, Matrix)
    return RowVecs(collect(X')), f, Σy
end

# Some type-piracy.
Base.isapprox(dx::NamedTuple{(:X,)}, dy::ColVecs) = isapprox(dx.X, dy.X)

Base.isapprox(dx::NamedTuple{(:X,)}, dy::RowVecs) = isapprox(dx.X, dy.X)

@testset "bayesian_linear_regression" begin
    @testset "blr $Tx" for Tx in [Matrix, ColVecs, RowVecs]
        @testset "consistency" begin
            rng, N, D, samples = MersenneTwister(123456), 11, 3, 1_000_000
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)

            AbstractGPs.TestUtils.test_finitegp_primary_and_secondary_public_interface(
                rng, f(X, Σy)
            )
        end
        @testset "rand" begin
            rng, N, D, samples = MersenneTwister(123456), 11, 3, 10_000_000
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)

            # Roughly test the statistical properties of rand.
            Y = rand(rng, f(X, Σy), samples)
            m_empirical = mean(Y; dims=2)
            Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
            @test mean(f(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
            @test cov(f(X, Σy)) ≈ Σ_empirical atol = 1e-3 rtol = 1e-3

            @testset "Zygote (everything dense)" begin
                function rand_blr(X, A_Σy, mw, A_Λw)
                    Σy, Λw = Symmetric(A_Σy * A_Σy' + I), Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    return rand(MersenneTwister(123456), f(X, Σy), 3)
                end
                mw, A_Σy, A_Λw = f.mw, 0.1 .* randn(rng, N, N), 0.1 .* randn(rng, D, D)

                # Run the model forwards and check that output agrees with non-Zygote output
                z, back = Zygote.pullback(rand_blr, X, A_Σy, mw, A_Λw)
                @test z == rand_blr(X, A_Σy, mw, A_Λw)

                # Compute adjoints using Zygote.
                z̄ = randn(rng, size(z))
                dX, dA_Σy, dmw, dA_Λw = back(z̄)

                # Verify adjoints via finite differencing.
                fdm = central_fdm(5, 1)
                @test dX ≈ first(j′vp(fdm, X -> rand_blr(X, A_Σy, mw, A_Λw), z̄, X))
                @test dA_Σy ≈
                    first(j′vp(fdm, A_Σy -> rand_blr(X, A_Σy, mw, A_Λw), z̄, A_Σy))
                @test dmw ≈ first(j′vp(fdm, mw -> rand_blr(X, A_Σy, mw, A_Λw), z̄, mw))
                @test dA_Λw ≈
                    first(j′vp(fdm, A_Λw -> rand_blr(X, A_Σy, mw, A_Λw), z̄, A_Λw))
            end
        end
        @testset "logpdf" begin
            rng, N, D = MersenneTwister(123456), 13, 7
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)
            y = rand(rng, f(X, Σy))

            # Construct MvNormal using a naive but simple computation for the mean / cov.
            function naive_normal_stats(X::Matrix)
                return (X' * f.mw, Symmetric(X' * (cholesky(f.Λw) \ X) + Σy))
            end
            naive_normal_stats(X::ColVecs) = naive_normal_stats(X.X)
            naive_normal_stats(X::RowVecs) = naive_normal_stats(collect(X.X'))
            m, Σ = naive_normal_stats(X)

            # Check that logpdf agrees between distributions and BLR.
            @test logpdf(f(X, Σy), y) ≈ logpdf(MvNormal(m, Σ), y)

            @testset "Zygote (everything dense)" begin
                function logpdf_blr(X, A_Σy, y, mw, A_Λw)
                    Σy, Λw = Symmetric(A_Σy * A_Σy' + I), Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    return logpdf(f(X, Σy), y)
                end
                mw, A_Σy, A_Λw = f.mw, 0.1 .* randn(rng, N, N), 0.1 .* randn(rng, D, D)

                z, back = Zygote.pullback(logpdf_blr, X, A_Σy, y, mw, A_Λw)
                @test z == logpdf_blr(X, A_Σy, y, mw, A_Λw)

                # Compute gradients using Zygote.
                z̄ = randn(rng)
                dX, dA_Σy, dy, dmw, dA_Λw = back(z̄)

                # Check correctness via finite differencing.
                fdm = central_fdm(5, 1)
                @test dX ≈ first(j′vp(fdm, X -> logpdf_blr(X, A_Σy, y, mw, A_Λw), z̄, X))
                @test dA_Σy ≈
                    first(j′vp(fdm, A_Σy -> logpdf_blr(X, A_Σy, y, mw, A_Λw), z̄, A_Σy))
                @test dy ≈ first(j′vp(fdm, y -> logpdf_blr(X, A_Σy, y, mw, A_Λw), z̄, y))
                @test dmw ≈ first(j′vp(fdm, mw -> logpdf_blr(X, A_Σy, y, mw, A_Λw), z̄, mw))
                @test dA_Λw ≈
                    first(j′vp(fdm, A_Λw -> logpdf_blr(X, A_Σy, y, mw, A_Λw), z̄, A_Λw))
            end
        end
        @testset "posterior" begin
            @testset "low noise" begin
                rng, N, D = MersenneTwister(123456), 13, 7
                X, f, Σy = generate_toy_problem(rng, N, D, Tx)
                y = rand(rng, f(X, eps()))

                f′_low_noise = posterior(f(X, eps()), y)
                @test mean(f′_low_noise(X, eps())) ≈ y
                @test all(cov(f′_low_noise(X, eps())) .< 1_000 * eps())
            end
            @testset "repeated conditioning" begin
                rng, N, D = MersenneTwister(123456), 13, 7
                X, f, Σy = generate_toy_problem(rng, N, D, Tx)
                X′ = randn(rng, D, N)
                y = rand(rng, f(X, Σy))

                # Chop up the noise because we can't condition on noise that's correlated
                # between things.
                N1 = N - 3
                Σ1, Σ2 = Σy[1:N1, 1:N1], Σy[(N1 + 1):end, (N1 + 1):end]
                Σy′ = vcat(hcat(Σ1, zeros(N1, N - N1)), hcat(zeros(N - N1, N1), Σ2))

                X1 = X isa Matrix ? X[:, 1:N1] : X[1:N1]
                X2 = X isa Matrix ? X[:, (N1 + 1):end] : X[(N1 + 1):end]
                y1, y2 = y[1:N1], y[(N1 + 1):end]

                f′1 = posterior(f(X1, Σ1), y1)
                f′2 = posterior(f′1(X2, Σ2), y2)
                f′ = posterior(f(X, Σy′), y)
                @test mean(f′(X′, Σy)) ≈ mean(f′2(X′, Σy))
                @test cov(f′(X′, Σy)) ≈ cov(f′2(X′, Σy))
            end
        end

        @testset "sampling functions" begin
            rng, N, D = MersenneTwister(123456), 11, 5
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)

            g = rand(rng, f)
            @test g(X) == g(X)  # check the sample doesn't change between evaluations

            if X isa Matrix
                Xc = ColVecs(X)
                Xr = RowVecs(X')
                @test g(X) == g(Xc)
                @test g(X) == g(Xr)
            end

            # test the Random interface
            @test rand(rng, Random.Sampler(rng, f, Val(Inf))) isa
                BayesianLinearRegressors.BLRFunctionSample

            samples1, samples2 = 1000, 1000
            samples = samples1 * samples2
            gs = rand(rng, f, samples1, samples2)
            @test size(gs) == (samples1, samples2)

            # test statistical properties of the sampled functions
            let
                Y = reduce(hcat, map(h -> h(X), reshape(gs, :)))
                m_empirical = mean(Y; dims=2)
                Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
                @test mean(f(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
                @test cov(f(X, Σy)) ≈ Σ_empirical + Σy atol = 1e-2 rtol = 1e-2
            end

            # test statistical properties of in-place rand
            let
                A = Array{BayesianLinearRegressors.BLRFunctionSample,2}(
                    undef, samples1, samples2
                )
                A = rand!(rng, A, f)
                Y = reduce(hcat, map(h -> h(X), reshape(gs, :)))
                m_empirical = mean(Y; dims=2)
                Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
                @test mean(f(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
                @test cov(f(X, Σy)) ≈ Σ_empirical + Σy atol = 1e-2 rtol = 1e-2
            end

            @testset "Zygote (everything dense)" begin
                function test_rand_funcs_adjoints(sample_function)
                    rng, N, D = MersenneTwister(123456), 11, 5
                    X, f, _ = generate_toy_problem(rng, N, D, Tx)
                    mw, A_Λw = f.mw, 0.1 .* randn(rng, D, D)

                    # Run the model forwards and check that output agrees with non-Zygote.
                    z, back = Zygote.pullback(sample_function, X, mw, A_Λw)
                    @test z == sample_function(X, mw, A_Λw)

                    # Compute adjoints using Zygote.
                    z̄ = randn(rng, size(z))
                    dX, dmw, dA_Λw = back(z̄)

                    # Verify adjoints via finite differencing.
                    fdm = central_fdm(5, 1)
                    @test dX ≈ first(j′vp(fdm, X -> sample_function(X, mw, A_Λw), z̄, X))
                    @test dmw ≈ first(j′vp(fdm, mw -> sample_function(X, mw, A_Λw), z̄, mw))
                    @test dA_Λw ≈
                        first(j′vp(fdm, A_Λw -> sample_function(X, mw, A_Λw), z̄, A_Λw))
                end

                function rand_funcs_single(X, mw, A_Λw)
                    Λw = Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    g = rand(MersenneTwister(123456), f)
                    return g(X)
                end

                function rand_funcs_multi(X, mw, A_Λw)
                    Λw = Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    gs = rand(MersenneTwister(123456), f, 1, 1)
                    return reduce(hcat, map(h -> h(X), reshape(gs, :)))
                end

                test_rand_funcs_adjoints(rand_funcs_single)
                test_rand_funcs_adjoints(rand_funcs_multi)
            end
        end
    end
    @testset "unrecognised AbstractVector" begin
        rng = MersenneTwister(123456)
        N, D = 11, 5
        x = collect(eachrow(randn(rng, N, D)))
        _, f, Σy = generate_toy_problem(rng, N, D, ColVecs)
        @test_throws ErrorException rand(f(x, Σy))
    end
end


text:
"""
Dict-like data structure which maps function signature to a value.
Unlike real dict, getindex(rsv, sig) returns either exact match, or
closest matching function signature. Example:

    rsv = FunctionResolver{Symbol}()
    rsv[Tuple{typeof(sin), Float64}] = :Float64
    rsv[Tuple{typeof(sin), Real}] = :Real
    rsv[Tuple{typeof(sin), Number}] = :Number

    rsv[Tuple{typeof(sin), Float64}]   # ==> :Float64
    rsv[Tuple{typeof(sin), Float32}]   # ==> :Real
"""
mutable struct FunctionResolver{T}
    signatures::Dict{Symbol, Vector{Pair{Any, T}}}
    ordered::Bool
    FunctionResolver{T}() where T = new{T}(Dict(), false)
end

# function FunctionResolver{T}(pairs::Vector{Pair{S, T}}) where {S, T}
# function FunctionResolver{T}(pairs::Vector{Pair{S, T} where S}) where T
function FunctionResolver{T}(pairs::Vector) where T
    rsv = FunctionResolver{T}()
    for (sig, val) in pairs
        rsv[sig] = val
    end
    order!(rsv)
    return rsv
end

Base.show(io::IO, rsv::FunctionResolver) = print(io, "FunctionResolver($(length(rsv.signatures)))")

function_type(@nospecialize sig) = sig isa UnionAll ? function_type(sig.body) : sig.parameters[1]
function_type_key(fn_typ) = Symbol("$(Base.parentmodule(fn_typ)).$(Base.nameof(fn_typ))")

function Base.setindex!(rsv::FunctionResolver{T}, val::T, @nospecialize sig::Type{<:Tuple}) where T
    fn_typ = function_type(sig)
    key = function_type_key(fn_typ)
    if !haskey(rsv.signatures, key)
        rsv.signatures[key] = Pair{Type, T}[]
    end
    pairs = rsv.signatures[key]
    # if such signature already exists, just replace the value
    updated = false
    for (i, (old_sig, _)) in enumerate(pairs)
        if sig == old_sig
            pairs[i] = sig => val
            updated = true
            break
        end
    end
    # otherwise push new pair and mark as unordered
    if !updated
        push!(pairs, sig => val)
        rsv.ordered = false
    end
    return val
end

function Base.getindex(rsv::FunctionResolver{T}, @nospecialize sig::Type{<:Tuple}) where T
    rsv.ordered || order!(rsv)
    fn_typ = function_type(sig)
    key = function_type_key(fn_typ)
    if haskey(rsv.signatures, key)
        for (TT, val) in rsv.signatures[key]
            if sig <: TT
                return val
            end
        end
    end
    return nothing
end

if Vararg isa Type
    # Vararg changes between julia 1.6 and julia 1.7
    is_Vararg(T) = T isa Type && T <: Vararg
else
    is_Vararg(T) = T isa typeof(Vararg)
end

function isless_signature(sig1, sig2)
    # signatures with Varargs should go last
    if any(is_Vararg, get_type_parameters(sig2))
        return true
    else
        return sig1 <: sig2
    end
end

function order!(rsv::FunctionResolver)
    for (fn_typ, sigs) in rsv.signatures
        sort!(sigs, lt=(p1, p2) -> isless_signature(p1[1], p2[1]))
    end
    rsv.ordered = true
end

Base.haskey(rsv::FunctionResolver, sig::Type{<:Tuple}) = (rsv[sig] !== nothing)
Base.in(sig::Type{<:Tuple}, rsv::FunctionResolver) = haskey(rsv, sig)
Base.empty!(rsv::FunctionResolver) = empty!(rsv.signatures)

function find_signatures_for(rsv::FunctionResolver, f::Union{Function, DataType})
    return rsv.signatures[function_type_key(typeof(f))]
end


text:
# Julia wrapper for header: /usr/include/libavutil/hwcontext_vdpau.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export


text:
using Base.Test
using Gumbo

include("basics.jl")
include("comparison.jl")
include("parsing.jl")
include("traversal.jl")


text:
@testcase "subtyping" begin
    @lift abstract type A{a} end
    @lift abstract type B{a} <: A{a} end

    @testset "adt subtying" begin
        @lift @data C{A} <: B{A} begin
            C1(A, Int)
        end

        @test C1(1, 2) isa C{Int}
    end
end

@testcase "adt list" begin
    # 1. define the List data type
    @lift @data List{T} begin
        Nil()
        Cons(head :: T, tail :: List{T})
    end
    # 2. define interfaces
    len(xs::List{T}) where T = @match xs begin
        Nil{T}() => 0
        Cons{T}(_, tail) => 1 + len(tail)
    end
    @testset "adt List" begin
        @test len(Nil{Any}()) == 0
        xs = Cons(3,Cons(2, Cons(1, Nil{Int}())))
        @test len(xs) == 3
    end
end

@testcase "adt arith" begin
    # 1. define Arith
    @lift @data Arith begin
        Num(v :: Int)
        Minus(fst :: Arith, snd :: Arith)
        Add(fst :: Arith, snd :: Arith)
        Mult(fst :: Arith, snd :: Arith)
        Divide(fst :: Arith, snd :: Arith)
    end
    # 2. define interfaces
    function eval_arith(arith :: Arith)
        @match arith begin
            Num(v)       => v
            Add(fst, snd) => eval_arith(fst) + eval_arith(snd)
            Minus(fst, snd) => eval_arith(fst) - eval_arith(snd)
            Mult(fst, snd)   => eval_arith(fst) * eval_arith(snd)
            Divide(fst, snd) => eval_arith(fst) / eval_arith(snd)
        end
    end

    @testset "adt Arith" begin
        Number = Num
        @test eval_arith(
            Add(Number(1),
                Minus(Number(2),
                    Divide(Number(20),
                            Mult(Number(2),
                                Number(5)))))) == 1
    end
end

@testcase "share data with several modules" begin
    @lift @data CD begin
        D(a, b)
        C{T} :: (a :: Int, b::T) => CD
    end
    @lift @data A begin
        E()
    end

    @testset "case" begin
        @test E <: A
        @test fieldnames(D) == (:a, :b)
        @test_throws MethodError C(3.0, :abc)
    end

    @lift module ADummy
        using MLStyle
    end

    @lift module BDummy
        using MLStyle
    end

    @lift @data visible in [ADummy] SSS begin
        SSS_1(Int)
    end

    ADummy.eval(:(SSS_1 = $SSS_1; SSS = $SSS))

    @test ADummy.eval(quote
        @match SSS_1(2) begin
            SSS_1(_) => :ok
        end
    end) == :ok

    BDummy.eval(:(SSS_1 = $SSS_1; SSS = $SSS))

    @test_skip BDummy.eval(quote
        @match SSS_1(2) begin
            SSS_1(_) => :ok
        end
    end)

    @testset "enum: #85" begin
        @lift @data Es{T} begin
            E1 :: Es{Int}
            E2 :: Es{Char}
        end

        @lift @data Es′ begin
            E3
            E4
        end
        e1 = E1
        e2 = E2
        e3 = E3
        e4 = E4

        @test @match 1 begin
            E1 => false
            E2 => false
            E3 => false
            E4 => false
            _ => true
        end

        @test @match E1 begin
            E1 => true
            _ => false
        end

        @test @match E2 begin
            E2 => true
            _ => false
        end

        @test @match E3 begin
            E3 => true
            _ => false
        end

        @test @match E4 begin
            E4 => true
            _ => false
        end

        @test @match E1 begin
            E2 => false
            _ => true
        end

        @test @match E3 begin
            E4 => false
            _ => true
        end
    end
end


text:
using LightGraphs, DifferentialEquations, Plots, LinearAlgebra

lil(a) = vcat(a, circshift(reverse(a),1))

function circulant_kura(a::Vector) # a is half vec
	n = 2length(a)
    g = SimpleGraph(reshape(vcat(map(x->circshift(lil(a), x), 0:(n-1))...), (n, n)))
    dgs = degree(g)
    @assert all(x->x==dgs[1], dgs)
    println("μ: $(dgs[1]/(n-1))")
    # println("fake connectivity: ≈$((round(ne / binomial(nv, 2), digits=2)))")
    g
end

# shitty wrapper
function kura_prob(u0, tspan, p)::ODEProblem
    ODEProblem(kuramoto!, u0, tspan, p)
end

function kuramoto!(du, u, p, t)
    n = length(u)
    A, ω = p # A is adjacency_matrix, and ω is
    u .+= ω * t
    for i in 1:n
        du[i] = sum(A[i, j] * sin(u[j] - u[i]) for j in 1:n)
    end
end

function anim_kura(u0, g)
    θn = θo = u0
    # p = scatter(θo)
    # p = plot(20)
    p = scatter(u0)
    anim = @animate for _ in 1:100
       θn = kuramoto(g, θo)
       println("$(θn - θo)")
       push!(p, 20, θn)
       θo = θn
    end
end

function demo(g)

    u0 = rand(nv(g))
    tspan = (0., 10.)
    p = (adjacency_matrix(g), 10.)
    ODEProblem(kuramoto!, u0, tspan, p)

end


text:
"""
    get_scaling_val(Nt, swap_strategy)

Calculates the correct scaling factor for polynomial step size between temperatures
"""
function get_scaling_val(Nt, swap_strategy)
    # Why these?
    if swap_strategy == :standard
        scaling_val = Nt - 1
    elseif swap_strategy == :nonrev
        scaling_val = 2
    else
        scaling_val = 1
    end
    return scaling_val
end

"""
    generate_Δ(Nt, swap_strategy)

Returns a temperature ladder `Δ` containing `Nt` temperatures,
generated in accordance with the chosen `swap_strategy`
"""
function generate_Δ(Nt, swap_strategy)
    scaling_val = get_scaling_val(Nt, swap_strategy)
    Δ = zeros(Real, Nt)
    Δ[1] = 1.0
    β′ = Δ[1]
    for i ∈ 1:(Nt - 1)
        β′ += exp(scaling_val)
        Δ[i + 1] = Δ[1] / β′
    end
    return Δ
end

"""
    check_Δ(Δ)

Checks and returns a sorted `Δ` containing `{β₀, ..., βₙ}` conforming such that `1 = β₀ > β₁ > ... > βₙ ≥ 0`
"""
function check_Δ(Δ)
    if !all(zero.(Δ) .≤ Δ .≤ one.(Δ))
        error("Temperature schedule provided has values outside of the acceptable range, ensure all values are in [0, 1].")
    end
    Δ = sort(Δ; rev=true)
    if Δ[1] != one(Δ[1])
        error("Δ must contain 1, as β₀.")
    end
    return Δ
end


text:
using Vann
using Base.Test

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

# Test of Gr4j

# Read data

filename = joinpath(Pkg.dir("Vann"), "data", "airgr", "test_data.txt")

data = readdlm(filename, ',', header = true)

prec  = data[1][:,1]
epot  = data[1][:,2]
q_obs = data[1][:,3]

prec = transpose(prec)
epot = transpose(epot)

frac = zeros(Float64, 1)

param  = [257.238, 1.012, 88.235, 2.208]

tstep = 24.0

time = DateTime(2000,1,1)

# Select model

st_gr4j = Gr4j(tstep, time, param)

# Run model

q_sim = run_model(st_gr4j, prec, epot)

# Compute largest error

err_max = maximum(map(abs, q_sim - q_obs))

@test err_max < 1.

println("Gr4j maximum error = " * string(err_max))

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

# Test of TinBasic and Gr4j

# Load data

filename = joinpath(Pkg.dir("Vann"), "data", "atnasjo")

date, tair, prec, q_obs, frac = load_data(filename, "Q_ref.txt")

tstep = 24.0

time = date[1]

# Compute potential evapotranspiration

epot = epot_zero(date)

# Parameters

param_snow  = [0.0, 3.69, 1.02]
param_hydro = [74.59, 0.81, 214.98, 1.24]

# Select model

st_snow  = TinBasic(tstep, time, param_snow, frac)
st_hydro = Gr4j(tstep, time, param_hydro)

q_sim = run_model(st_snow, st_hydro, tair, prec, epot)

# Compute largest error

err_max = maximum(map(abs, q_sim - q_obs))

println("Gr4j + TinBasic maximum error = " * string(err_max))

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

# Test of TinStandard and Gr4j

# Load data

filename = joinpath(Pkg.dir("Vann", "data", "atnasjo"))

date, tair, prec, q_obs, frac = load_data(filename, "Q_ref.txt")

tstep = 24.0

time = date[1]

# Compute potential evapotranspiration

epot = epot_zero(date)

# Parameters

param_snow  = [0.0, 3.69, 3.69, 0.0, 1.02]
param_hydro = [74.59, 0.81, 214.98, 1.24]

# Select model

st_snow  = TinStandard(tstep, time, param_snow, frac)
st_hydro = Gr4j(tstep, time, param_hydro)

q_sim = run_model(st_snow, st_hydro, tair, prec, epot)

# Compute largest error

err_max = maximum(map(abs, q_sim - q_obs))

println("Gr4j + TinStandard maximum error = " * string(err_max))

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann