text:
##########################################
## define dictionary inverting function ##
##########################################

function invertDict(dict::Dict)
    invDict = Dict()
    for val in unique(values(dict))
        ## for each possible value, find all keys
        ## @show val
        keysWithGivenValue = Array(Symbol, 0)
        ## counter = 1
        for key in dict
            ## @show key
            if key[2] == val
                push!(keysWithGivenValue, key[1])
                ## counter = counter + 1
            end
        end
        invDict[val] = keysWithGivenValue
    end
    return invDict
end

#########################################
## split-apply-combine functionalities ##
#########################################

## grouping columns:
## - grouping is determined through dictionary
## - grouping occurs with respect to columns as opposed to rows
## - grouping will return dimensions nObs x nGroups
## - for each group, dictionary lists all members of the group
## - function should return scalar value and require first input

function groupAssets(invs::Investments, groups::Dict{Any, Any})
    ## aggregate weights with respect to group membership
    ##
    ## Inputs:
    ##	invs 				nObs x nAss Investments
    ## 	groups 			nGroups x 1 Dict{Symbol, Array{Symbol}}
    ## 						specifying group membership
    ##
    ## Outputs:
    ## 	groupedInvs 	nObs x nGroups group Investments

    wgts = invs.vals
    groupWgts = DataFrame()
    for group in groups
        ## group[2] is Array of symbols
        coreWgts = array(wgts[group[2]])
        groupWgts[group[1]] = sum(coreWgts, 2)[:]
    end
    AssetMgmt.Investments(groupWgts, idx(invs))
end

function groupAssets(invs::Investments, groups::Dict{Any, Any},
                     f::Function)
    ## apply function to groups of assets
    ##
    ## Inputs:
    ##	invs 				nObs x nAss Investments
    ## 	groups 			nGroups x 1 Dict{Symbol, Array{Symbol}}
    ## 						specifying group membership
    ## 	f					function to be applied on groups, return
    ## 						column vector of scalar values per group
    ##
    ## Outputs:
    ## 	groupedInvs 	nObs x nGroups Timematr / Timedata

    wgts = invs.vals
    groupWgts = DataFrame()
    for group in groups
        ## group[2] is Array of symbols
        coreWgts = array(wgts[group[2]])
        groupWgts[group[1]] = f(coreWgts)[:]
    end
    AssetMgmt.Timematr(groupWgts, idx(invs))
end

function groupVars(tm::Timematr, groups::Dict{Any, Any},
                     f::Function)
    ## apply function to groups of assets
    ##
    ## Inputs:
    ##	tm 				nObs x nAss Timematr
    ## 	groups 			nGroups x 1 Dict{Symbol, Array{Symbol}}
    ## 						specifying group membership
    ## 	f					function to be applied on groups, return
    ## 						column vector of scalar values per group
    ##
    ## Outputs:
    ## 	groupedInvs 	nObs x nGroups Timematr / Timedata

    wgts = tm.vals
    groupWgts = DataFrame()
    for group in groups
        ## group[2] is Array of symbols
        coreWgts = array(wgts[group[2]])
        groupWgts[group[1]] = f(coreWgts)[:]
    end
    AssetMgmt.Timematr(groupWgts, idx(tm))
end

function sectorPerformances(invs::Investments, discRet::Timematr,
                           assetsInSector::Dict{Any, Any})
    ## calculate (normalized) sector return
    ## for each sector, pretend that one is only invested in this
    ## sector. hence, relative weights within sector are rescaled, and
    ## the result is the performance of the sector if one is invested
    ## according to the proportions chosen by the larger investment
    ## strategy. this could be compared to sector performances with
    ## equal weights: does chosen strategy chose assets within each
    ## sector sophisticated enough to improve sector performance?
    ##
    ## PROBLEM: how to cope with situations where all weights of a
    ## given sector are negative? what would be the associated
    ## normalization?
    ##
    ## Inputs:
    ## 	invs 				nObs x nAss Investments of portfolio weights
    ## 	discRet			nObs x nAss Timematr of asset returns
    ## 	assetsInSector	nGroups x 1 Dict{Symbol, Array{Symbol}}
    ## 						specifying group membership
    ##
    ## Outputs:
    ## 	nObs x nGroups Timematr of sector returns

    ## get overall sector weights for normalization
    sectorWgts = AssetMgmt.groupAssets(invs, assetsInSector)

    ## transform investments to Timematr to allow mathematical
    ## functions
    invsTm = convert(Timematr, invs)

    ## preallocation
    sectorPortfolios = DataFrame(size(sectorWgts)...)
    rename!(sectorPortfolios, names(sectorPortfolios), names(sectorWgts))

    for sect in assetsInSector
        currentSectorWgts = invsTm[sect[2]]
        associatedReturns = discRet[:, sect[2]]
        unscaledSectorReturn = rowsums(currentSectorWgts .*
                                       associatedReturns)
        sectorPortfolios[sect[1]] = (core(unscaledSectorReturn) ./
                                     core(sectorWgts[sect[1]]))[:]
    end
    return Timematr(sectorPortfolios, idx(invs))
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 TunnelConnectionHealth <: SwaggerModel
    tunnel::Any # spec type: Union{ Nothing, String } # spec name: tunnel
    connectionStatus::Any # spec type: Union{ Nothing, String } # spec name: connectionStatus
    ingressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: ingressBytesTransferred
    egressBytesTransferred::Any # spec type: Union{ Nothing, Int64 } # spec name: egressBytesTransferred
    lastConnectionEstablishedUtcTime::Any # spec type: Union{ Nothing, String } # spec name: lastConnectionEstablishedUtcTime

    function TunnelConnectionHealth(;tunnel=nothing, connectionStatus=nothing, ingressBytesTransferred=nothing, egressBytesTransferred=nothing, lastConnectionEstablishedUtcTime=nothing)
        o = new()
        validate_property(TunnelConnectionHealth, Symbol("tunnel"), tunnel)
        setfield!(o, Symbol("tunnel"), tunnel)
        validate_property(TunnelConnectionHealth, Symbol("connectionStatus"), connectionStatus)
        setfield!(o, Symbol("connectionStatus"), connectionStatus)
        validate_property(TunnelConnectionHealth, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        setfield!(o, Symbol("ingressBytesTransferred"), ingressBytesTransferred)
        validate_property(TunnelConnectionHealth, Symbol("egressBytesTransferred"), egressBytesTransferred)
        setfield!(o, Symbol("egressBytesTransferred"), egressBytesTransferred)
        validate_property(TunnelConnectionHealth, Symbol("lastConnectionEstablishedUtcTime"), lastConnectionEstablishedUtcTime)
        setfield!(o, Symbol("lastConnectionEstablishedUtcTime"), lastConnectionEstablishedUtcTime)
        o
    end
end # type TunnelConnectionHealth

const _property_map_TunnelConnectionHealth = Dict{Symbol,Symbol}(Symbol("tunnel")=>Symbol("tunnel"), Symbol("connectionStatus")=>Symbol("connectionStatus"), Symbol("ingressBytesTransferred")=>Symbol("ingressBytesTransferred"), Symbol("egressBytesTransferred")=>Symbol("egressBytesTransferred"), Symbol("lastConnectionEstablishedUtcTime")=>Symbol("lastConnectionEstablishedUtcTime"))
const _property_types_TunnelConnectionHealth = Dict{Symbol,String}(Symbol("tunnel")=>"String", Symbol("connectionStatus")=>"String", Symbol("ingressBytesTransferred")=>"Int64", Symbol("egressBytesTransferred")=>"Int64", Symbol("lastConnectionEstablishedUtcTime")=>"String")
Base.propertynames(::Type{ TunnelConnectionHealth }) = collect(keys(_property_map_TunnelConnectionHealth))
Swagger.property_type(::Type{ TunnelConnectionHealth }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_TunnelConnectionHealth[name]))}
Swagger.field_name(::Type{ TunnelConnectionHealth }, property_name::Symbol) =  _property_map_TunnelConnectionHealth[property_name]

const _allowed_TunnelConnectionHealth_connectionStatus = ["Unknown", "Connecting", "Connected", "NotConnected"]

function check_required(o::TunnelConnectionHealth)
    true
end

function validate_property(::Type{ TunnelConnectionHealth }, name::Symbol, val)
    if name === Symbol("connectionStatus")
        Swagger.validate_param(name, "TunnelConnectionHealth", :enum, val, _allowed_TunnelConnectionHealth_connectionStatus)
    end
end

ext: jl


text:
nax = [CartesianIndex()]
# include("./antenna.jl")
# include("./parser.jl")

"""
For Relative indexing like arr[1:-1]
"""
struct Slice
    start::Int
#    step::Int
    stop::Int
end

# const sslice = ⁝
⁝(start::Int=1, stop::Int=0) = Slice(start, stop)
#⁝(start::Int=1, step::Int=0, stop::Int=0) = Slice(start, 1, stop)

function Base.getindex(A::T, s::Slice) where {T<:AbstractArray}
    N = size(A)[1]
    a = s.start < 1 ? N + s.start : s.start
    b = s.stop < 1 ? N + s.stop : s.stop

    R = ndims(A)
    return if R == 3
        getindex(A, a:b, :, :)
    elseif R == 2
        getindex(A, a:b, :)
    elseif R == 1
        getindex(A, a:b)
    end
end

"""
coord : 3xA array
cell : 3x3  ase-like array
"""
function cart2scaled(coord::Array, cell::Array)
    return cell' \ coord
end

"""
coords : 3xAxN array
cell : 3x3  ase-like array
return : 3xAxN array
"""
function cart2scaled(coords::Array{Float64, 3}, cell::Array{Float64, 2})
    scaleds = zeros(Float64, size(coords)...)
    # cell_t = cell'
    cell⁻¹ = inv(cell')
    for i = 1:size(coords)[end]
        scaleds[:, :, i] =  cell⁻¹ * coords[:, :, i]
    end
    return scaleds
end

"""
scaled : 3xA array
cell : 3x3  ase-like array
"""
function scaled2cart(scaled::Array{Float64, 2}, cell::Array{Float64, 2})
    return cell' * scaled
end

"""
scaled : 3xAxN array
cell : 3x3  ase-like array
return : 3xAxN array
"""
function scaled2cart(scaleds::Array{Float64, 3}, cell::Array{Float64, 2})
    coords = zeros(Float64, size(scaleds)...)
    for i = 1:size(scaleds)[end]
        coords[:, :, i] = cell' * scaleds[:, :, i]
    end
    return coords
end

"""
coord : 3xA array
cell : 3x3  ase-like array
n : Int index of atom to be centered
return : 3xA array
"""
function center_wrap(coord::Array{Float64, 2}, cell::Array{Float64, 2}, n::Int)
    scaled = cart2scaled(coord, cell)
    scaled .-= ((scaled[:, n] .+ 0.5) .% 1.) .- 0.5
    return scaled2cart(scaled, cell)
end

"""
coords : 3xAxN array
cell : 3x3  ase-like array
n : Int; index of atom to be centered
return : 3xAxN array
"""
function center_wrap(coords::Array{Float64, 3}, cell::Array{Float64, 2}, n::Int)
    scaleds = cart2scaled(coords, cell)
    scaleds .-= ((scaleds[:, n, nax, :] .+ 0.5) .% 1.) .- 0.5
    return scaled2cart(scaleds, cell)
end

"""
period = [false, false, false]
periodic = supercell_periodic_help(period)
"""
function supercell_periodic_help(period)
    per = Set([])
    for x in [-1, 0, 1] for y in [-2, 0, 2] for z in [3, 0, -3]
        dir = [x, y, z][period]
        dir == [0, 0, 0][period] ? continue : nothing
        union!(per, Set([dir]))
    end end end
    hcat([p for p in per]...)
end

"""
Minimum image convention cell padding
"""
function getmiccell(coords, cell, period)
    N, A, _ = size(coords)
    periodic = supercell_periodic_help(period)

    pD, pN = periodic == [] ? (0, 0) : size(periodic)
    supercell = similar(coords, 3, A * (1+pN), N)
    positionss = permutedims(coords, (3, 2, 1))
    supercell[:, 1:A, :] = positionss
    for i=1:pN
        translate = zeros(3)
        for j=1:pD
            p = periodic[j, i]
            p == 0 ? continue : nothing
            s = sign(p)
            idx = abs(p)
            translate += s .* cell[idx, :]
        end
        scell = positionss .+ translate
        supercell[:, A*i+1:A*(i+1), :] = scell
    end
    return supercell
end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
"""
    mutable struct SolidStateDetector{T <: SSDFloat, CS} <: AbstractConfig{T}

CS: Coordinate System: -> :cartesian / :cylindrical
"""
mutable struct SolidStateDetector{T <: SSDFloat, CS} <: AbstractConfig{T}
    name::String  # optional
    inputunits::Dict{String, Unitful.Units}
    world::World{T, 3}

    config_dict::Dict

    medium::NamedTuple # this should become a struct at some point

    semiconductors::Vector{Semiconductor{T}}
    contacts::Vector{Contact{T}}
    passives::Vector{Passive{T}}

    SolidStateDetector{T, CS}() where {T <: SSDFloat, CS} = new{T, CS}()
end

get_precision_type(d::SolidStateDetector{T}) where {T} = T
get_coordinate_system(d::SolidStateDetector{T, CS}) where {T, CS} = CS

function construct_units(dict::Dict)::Dict{String,Unitful.Units}
    result_dict::Dict{String,Unitful.Units} = Dict()
    haskey(dict,"length") ? result_dict["length"] = unit_conversion[dict["length"]] : result_dict["length"] = u"mm"
    haskey(dict,"angle") ? result_dict["angle"] = unit_conversion[dict["angle"]] : result_dict["angle"] = u"rad"
    haskey(dict,"potential") ? result_dict["potential"] = unit_conversion[dict["potential"]] : result_dict["potential"] = u"V"
    haskey(dict,"temperature") ? result_dict["temperature"] = unit_conversion[dict["temperature"]] : result_dict["temperature"] = u"K"
    result_dict
end

function construct_semiconductor(T, sc::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Semiconductor{T}(sc, inputunit_dict)
end

function construct_passive(T, pass::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Passive{T}(pass, inputunit_dict)
end

function construct_contact(T, contact::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Contact{T}(contact, inputunit_dict)
end

function construct_objects(T, objects::Vector, semiconductors, contacts, passives, inputunit_dict)::Nothing
    for obj in objects
        if obj["type"] == "semiconductor"
            push!(semiconductors, construct_semiconductor(T, obj, inputunit_dict))
        elseif obj["type"] == "contact"
            push!(contacts, construct_contact(T, obj, inputunit_dict))
        elseif obj["type"] == "passive"
            push!(passives, construct_passive(T, obj, inputunit_dict))
        else
            @warn "please spcify the calss to bei either a \"semiconductor\", a \"contact\", or \"passive\""
        end
    end
    nothing
end

function SolidStateDetector{T}(config_file::Dict)::SolidStateDetector{T} where{T <: SSDFloat}
    grid_type = Symbol(config_file["setup"]["grid"]["coordinates"])
    c = SolidStateDetector{T, grid_type}()
    c.name = config_file["name"]
    c.config_dict = config_file
    c.inputunits = construct_units(config_file["setup"]["units"])
    c.world = World(T, config_file["setup"]["grid"], c.inputunits)
    c.medium = material_properties[materials[config_file["setup"]["medium"]]]
    c.semiconductors, c.contacts, c.passives = [], [], []
    construct_objects(T, config_file["setup"]["objects"], c.semiconductors, c.contacts, c.passives, c.inputunits)
    return c
end

function SolidStateDetector(parsed_dict::Dict)
    SolidStateDetector{Float32}(parsed_dict)
end

function Base.sort!(v::AbstractVector{<:AbstractGeometry})
    hierarchies::Vector{Int} = map(x->x.hierarchy,v)
    v_result::typeof(v) = []
    for idx in sort!(unique!(hierarchies))
        push!(v_result,filter(x->x.hierarchy == hierarchies[idx],v)...)
    end
    return v_result
end

function SolidStateDetector{T}(parsed_dict::Dict) where T
    SolidStateDetector{T}(parsed_dict)
end

function contains(c::SolidStateDetector, point::AbstractCoordinatePoint{T,3})::Bool where T
    for contact in c.contacts
        if point in contact
            return true
        end
    end
    for sc in c.semiconductors
        if point in sc
            return true
        end
    end
    return false
end

function println(io::IO, d::SolidStateDetector{T, CS}) where {T <: SSDFloat, CS}
    println("________"*d.name*"________\n")
    # println("Class: ",d.class)
    println("---General Properties---")
    println("-Environment Material: \t $(d.medium.name)")
    println("-Grid Type: \t $(CS)")
    println()
    println("# Semiconductors: $(length(d.semiconductors))")
    for (isc, sc)  in enumerate(d.semiconductors)
        println("\t_____Semiconductor $(isc)_____\n")
        println(sc)
    end
    println()
    println("# Contacts: $(length(d.contacts))")
    if length(d.contacts)<=5
        for c in d.contacts
            println(c)
        end
    end
    println()
    println("# Passives: $(length(d.passives))")
    if length(d.passives)<=5
        for p in d.passives
            # println(c)
        end
    end
end

function show(io::IO, d::SolidStateDetector{T}) where {T <: SSDFloat} println(d) end
function print(io::IO, d::SolidStateDetector{T}) where {T <: SSDFloat} println(d) end
function display(io::IO, d::SolidStateDetector{T} ) where {T <: SSDFloat} println(d) end
function show(io::IO,::MIME"text/plain", d::SolidStateDetector) where {T <: SSDFloat}
    show(io, d)
end

# ToDo: Test it
function generate_random_startpositions(d::SolidStateDetector{T}, n::Int, Volume::NamedTuple=bounding_box(d), rng::AbstractRNG = MersenneTwister(), min_dist_from_boundary = 0.0001) where T
    delta = T(min_dist_from_boundary)
    n_filled::Int = 0
    positions = Vector{CartesianPoint{T}}(undef,n)
    while n_filled < n
        sample=CylindricalPoint{T}(rand(rng,Volume[:r_range].left:0.00001:Volume[:r_range].right),rand(rng,Volume[:φ_range].left:0.00001:Volume[:φ_range].right),rand(rng,Volume[:z_range].left:0.00001:Volume[:z_range].right))
        if !(sample in d.contacts) && contains(d,sample) && contains(d,CylindricalPoint{T}(sample.r+delta,sample.φ,sample.z))&& contains(d,CylindricalPoint{T}(sample.r-delta,sample.φ,sample.z))&& contains(d,CylindricalPoint{T}(sample.r,sample.φ,sample.z+delta))&& contains(d,CylindricalPoint{T}(sample.r,sample.φ,sample.z-delta))
            n_filled += 1
            positions[n_filled]=CartesianPoint(sample)
        end
    end
    positions
end


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

@doc raw"""PersistentVolumeClaimStatus is the current status of a persistent volume claim.

    IoK8sApiCoreV1PersistentVolumeClaimStatus(;
        accessModes=nothing,
        capacity=nothing,
        conditions=nothing,
        phase=nothing,
    )

    - accessModes::Vector{String} : AccessModes contains the actual access modes the volume backing the PVC has. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#access-modes-1
    - capacity::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Represents the actual resources of the underlying volume.
    - conditions::Vector{IoK8sApiCoreV1PersistentVolumeClaimCondition} : Current Condition of persistent volume claim. If underlying persistent volume is being resized then the Condition will be set to 'ResizeStarted'.
    - phase::String : Phase represents the current phase of PersistentVolumeClaim.
"""
mutable struct IoK8sApiCoreV1PersistentVolumeClaimStatus <: SwaggerModel
    accessModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: accessModes
    capacity::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: capacity
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PersistentVolumeClaimCondition} } # spec name: conditions
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase

    function IoK8sApiCoreV1PersistentVolumeClaimStatus(;accessModes=nothing, capacity=nothing, conditions=nothing, phase=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("accessModes"), accessModes)
        setfield!(o, Symbol("accessModes"), accessModes)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaimStatus

const _property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus = Dict{Symbol,Symbol}(Symbol("accessModes")=>Symbol("accessModes"), Symbol("capacity")=>Symbol("capacity"), Symbol("conditions")=>Symbol("conditions"), Symbol("phase")=>Symbol("phase"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaimStatus = Dict{Symbol,String}(Symbol("accessModes")=>"Vector{String}", Symbol("capacity")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("conditions")=>"Vector{IoK8sApiCoreV1PersistentVolumeClaimCondition}", Symbol("phase")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaimStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaimStatus)
    true
end

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


text:
using IterTools
using Statistics
include("Solution.jl")
include("Params.jl")
include("Recombinations.jl")

""" Genetic Algorithm
Initialize population with random candidate solutions
Evaluate each candidate
Repeat until Termination condition is satisfied
    Select parents (best fitness)
    Recombine pairs of solutions (Crossover)
    Mutate the resulting offspring (Mutation)
    Evaluate new candidates
    Select individuals for the next generation
"""

# Genetic Algorithm
function runGA(params::Params, maxGenerations::Int, populationSize::Int;
               recombineOp::Function=order1cx, distanceFunc::Function=hamming, entropyReduceFunc::Function=mean,
               mutate::Any=false, selectionOp::Function=fitnessSelection, selectionOpts::Dict{String,Any}=Dict{String,Any}("fitnessThreshold"=>0.2))
    selectionOpts["n"] = populationSize
    selectionOpts["distanceFunc"] = hamming
    startTime = time()
    # Check if population size is an even number
    if populationSize % 2 != 0
        populationSize += 1
    end
    results = Dict{String,Any}()
    results["maxGenerations"] = maxGenerations
    results["populationSize"] = populationSize
    results["entropy"] = []
    results["numIndividuals"] = []
    results["cost"] = []
    results["generationTime"] = []

    # Initialize population with random candidate solutions
    population = initializePopulation(params, populationSize)
    bestSol = sort(population)[1]

    append!(results["numIndividuals"],length(population))
    append!(results["entropy"],entropy(population, distanceFunc, entropyReduceFunc))
    append!(results["cost"],bestSol.cost)
    append!(results["generationTime"],time() - startTime)

    generation = 0
    # Repeat until Termination condition is satisfied
    while ((generation < maxGenerations) && (time() - startTime < params.maxTime) && (bestSol.cost != params.solutionCost))
        # Select parents
        population = shuffle!(population)  # shuffle population
        parents = [(population[i], population[i+Int(floor(length(population)/2))])
                   for i=1:Int(floor(length(population)/2))]

        # Recombine pairs of solutions
        offspring = recombinationOp(parents, recombineOp)

        # Mutate the resulting offspring
        if mutate != false
            offspring = mutationOp(offspring, mutate)
        end

        # Evaluate new candidates

        # Select individuals for the next generation
        population = vcat(population,offspring)
        survivors = selectionOp(population, selectionOpts)
        population = population[survivors]
        population = removeClones(population) # Remove clones from the population

        sol = sort(population)[1]
        # record best observed solution
        if sol < bestSol
            bestSol = sol
        end

        append!(results["numIndividuals"],length(population))
        append!(results["entropy"],entropy(population, distanceFunc, entropyReduceFunc))
        append!(results["cost"],bestSol.cost)
        append!(results["generationTime"],time() - startTime)
        results["generations"] = generation
        generation += 1
    end
    # Results
    results["time"] = time() - startTime
    results["finalCost"] = bestSol.cost
    results["permutation"] = bestSol.permutation
    if isnan(params.solutionCost)
        results["gap"] = NaN
    else
        results["gap"] = (results["cost"] .- params.solutionCost) ./ params.solutionCost
    end
    return bestSol, results
end

function initializePopulation(params::Params, populationSize::Int)::Array{Solution}
    population = Array{Solution}(undef, populationSize)
    for k = 1:populationSize
        population[k] = Solution(params)
    end
    return population
end

##############################################################################
# Mutation operation

function mutationOp(population::Array{Solution}, mutate::Function=descentHeuristic)::Array{Solution}
    for k = 1:length(population)
        population[k] = mutate(population[k], max_iter=1.0)
    end
    return population
end

##############################################################################
# Selection operation

# Selection based on fitness. Selects n most fit individuals from the population.
function fitnessSelection(population::Array{Solution}, selectionOpts::Dict)::Array{Bool}
    costArray = [s1.cost for s1 in population]
    fitnessRank = getRanking(costArray)
    survivors = fitnessRank .<= selectionOpts["n"]
    return survivors
    #sort!(population)
    #population = population[1:n]
    #return population
end

# Selection based on biased fitness
function biasedFitnessSelection(population::Array{Solution}, selectionOpts::Dict)::Array{Bool}
    index, diversityScore, fitnessScore = fitnessDiversityScore(population, selectionOpts["distanceFunc"])
    biasedFitnessVal = biasedFitness(fitnessScore, diversityScore, selectionOpts["fitnessThreshold"])
    biasedFitnessRank = getRanking(biasedFitnessVal)
    #survivors = biasedFitnessVal .<= 1
    survivors = biasedFitnessRank .<= selectionOpts["n"]
    return survivors
    #population = population[survivors]
    #return population
end

function removeClones(population::Array{Solution})
    return unique(x -> x.permutation, population)
end

##############################################################################
# Distances
# Implemented distance functions:
# - Hamming distance
# - Kendall-Tau distance

# Hamming distance: number of positions with different values
function hamming(s1::Solution,s2::Solution)::Float64
    distance = sum(s1.permutation .!= s2.permutation)
    return distance
end

# Distance between two permutations: https://math.stackexchange.com/questions/2492954/distance-between-two-permutations

# Kendall-Tau distance:
# - https://en.wikipedia.org/wiki/Kendall_tau_distance
# - https://stats.stackexchange.com/questions/168602/whats-the-kendall-taus-distance-between-these-2-rankings/238852#238852
function kendallTau(s1::Solution,s2::Solution)::Float64
    A = s1.permutation
    B = s2.permutation
    pairs = subsets(1:length(A), 2)
    distance = 0
    for (x,y) in pairs
        a = A[x] - A[y]
        b = B[x] - B[y]
        # if discordant (different signs)
        if (a * b < 0)
            distance += 1
        end
    end
    return distance
end

# "Population entropy": maximum/average/median distance between any two solutions in the population
function entropy(population::Array{Solution}, distance::Function, reduce::Function)
    allpairs = subsets(population, 2)
    distances = zeros(Float64, length(allpairs)).-1
    for (i,pairs) in enumerate(allpairs)
        s1,s2 = pairs
        distances[i] = distance(s1,s2)
    end
    if ~all(distances .>= 0)
        throw("All distances should be non negative!")
    end
    entropy = reduce(distances)
    return entropy
end

##############################################################################
# Diversity management

function getRanking(input; rev=false)
    output = zeros(Int64,length(input))
    for (i,x) in enumerate(sort(1:length(input), by=y->input[y], rev=rev))
        output[x] = i
    end
    return output
end

# Individual contribution to diversity: mean distance from all other solutions
function diversityContrib(s1::Solution, population::Array{Solution}, distance::Function)
    dcontrib = mean([distance(s1,s2) for s2 in population if s1.permutation != s2.permutation])
    return dcontrib
end

# Ranks solutions by fFitness and diversity
function fitnessDiversityRank(population::Array{Solution}, distanceFunc::Function)
    """ Fitness Rank and Diversity Rank
    Fitness and Diversity Rank are values between 1 and the size of the population n
    with 1 meaning that the solution is the most fit/diverse
    and n meaning that the solution is the least fit/diverse
    """
    index = 1:length(population)
    costArray = [s1.cost for s1 in population]
    diversityArray = [diversityContrib(s1, population, distanceFunc) for s1 in population]
    fitnessRank = getRanking(costArray)
    diversityRank = getRanking(diversityArray,rev=true)
    # Check fitness ranking
    @assert population[fitnessRank .== 1][1].cost == minimum(costArray)
    @assert population[fitnessRank .== length(population)][1].cost == maximum(costArray)
    # Check diversity ranking
    @assert diversityArray[diversityRank .== 1][1] == maximum(diversityArray)
    @assert diversityArray[diversityRank .== length(population)][1] == minimum(diversityArray)
    return index, diversityRank, fitnessRank
end

function minMaxScale(x)
    x_norm = (x .- minimum(x)) ./ (maximum(x) - minimum(x))
    return x_norm
end

# Fitness and diversity score: fitness and diversity ranks normalized to values between 0 and 1
function fitnessDiversityScore(population::Array{Solution}, distanceFunc::Function)
    """ Fitness Score and Diversity Score
    Fitness and Diversity Score are values between 0 and 1
    with 0 meaning that the solution is the most fit/diverse
    and 1 meaning that the solution is the least fit/diverse
    """
    index, diversityRank, fitnessRank = fitnessDiversityRank(population, distanceFunc)
    diversityScore = minMaxScale(diversityRank)
    fitnessScore = minMaxScale(fitnessRank)
    return index, diversityScore, fitnessScore
end

# Biased fitness
function biasedFitness(fitnessScore::Array{Float64}, diversityScore::Array{Float64}, fitnessThreshold::Float64=0.2)
    biasedFitness = fitnessScore .+ (1 - fitnessThreshold).*diversityScore
    return biasedFitness
end

# TODO: function tests, especially for the recombination operations

# TODO: formulate and solve the Mathematical Programming model of the QAP. Compare time and objective values.

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
using DashBootstrapComponents, DashCoreComponents, DashHtmlComponents

accordion = html_div([
    dbc_accordion(
        [
            dbc_accordionitem(
                "This is the content of the first section",
                title = "Item 1",
                item_id = "item-1",
            ),
            dbc_accordionitem(
                "This is the content of the second section",
                title = "Item 2",
                item_id = "item-2",
            ),
            dbc_accordionitem(
                "This is the content of the third section",
                title = "Item 3",
                item_id = "item-3",
            ),
        ],
        id = "accordion",
        active_item = "item-1",
    ),
    html_div(id = "accordion-contents", className = "mt-3"),
])

callback!(
    app,
    Output("accordion-contents", "children"),
    Input("accordion", "active_item"),
) do item
    return "Item selected: $item"
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
""" Kalman Filter and Extended Kalman Filter Likelihoods """

function kf_likelihood(θ, param_kf::ParamKalmanFilter, state_beliefs::AbstractArray,
    action_history::AbstractArray, measurement_history::AbstractArray)
    # drop initial s0 belief
    @assert length(action_history) == length(measurement_history)
    N = length(measurement_history)
    # initialize log likelihood
    # l = o.R[1]
    l=0.0
    for (k, (s, y, u)) in enumerate(zip(state_beliefs, measurement_history, action_history))
        x_hat = param_kf.A(θ)*s.x + param_kf.B(θ)*u # predicted state prior
        P_hat = param_kf.A(θ)*s.P*param_kf.A(θ)' + param_kf.Q(θ) # a priori state covariance
        v = y - param_kf.H(θ)*x_hat # measurement pre fit residual
        S = param_kf.H(θ)*P_hat*param_kf.H(θ)' + param_kf.R(θ) # pre fit residual covariance
        l += 1/2*(v'*inv(S)*v + log(det(S)))
    end
    return l
end

function ekf_likelihood(θ, param_ekf::ExtendedParamKalmanFilter, state_beliefs::AbstractArray,
    action_history::AbstractArray, measurement_history::AbstractArray)
    # drop initial s0 belief
    @assert length(action_history) == length(measurement_history)
    N = length(measurement_history)
    # initialize log likelihood
    # l = o.R[1]
    l=0.0
    for (k, (s, y, u)) in enumerate(zip(state_beliefs, measurement_history, action_history))
        x_hat = param_ekf.f(θ, s.x, u) # predicted state prior
        F = ForwardDiff.jacobian(μ -> param_ekf.f(θ, μ, u), s.x)
        P_hat = F*s.P*F' + param_ekf.Q(θ) # a priori state covariance
        v = y - param_ekf.h(θ, x_hat) # measurement pre fit residual
        H = ForwardDiff.jacobian(μ -> param_ekf.h(θ, μ), s.x)
        S = H*P_hat*H' + param_ekf.R(θ) # pre fit residual covariance
        l += 1/2*(v'*inv(S)*v + log(det(S)))
    end
    return l
end


text:
using Test
using JuliaEuler

@testset "Examples 25 - 50" begin

    @testset "Ex25" begin
        @test ex25(3) == 12
        @test ex25(1000) == 4782

    end

    @testset "Ex26" begin
        @test ex26(10) == 7
        @test ex26(1000) == 983
    end

    @testset "Ex27" begin
        @test ex27(1000,1000) == -59231
        @test ex27(80,1602) == -126479
    end

    @testset "Ex28" begin
        @test ex28(3) == 25
        @test ex28(5) == 101
        @test ex28(7) == 261
        @test ex28(1001) == 669171001
    end

    @testset "Ex29" begin
        @test ex29(5,5) == 15
        @test ex29(100,100) == 9183
    end
end

max_stars_repo_name: binnisb/JuliaEuler

max_issues_repo_name: binnisb/JuliaEuler

max_forks_repo_name: binnisb/JuliaEuler


text:
import Base: length, start, next, done, eachindex, get, getindex, setindex!, show, display, keys, haskey, merge, similar, hash, ==, isequal, endof, merge

# Types
FrameProps = Dict{Symbol, Any}

immutable Frame <: Associative{Symbol, Any}
    props::FrameProps
end

AtomProps = Dict{Symbol, Array}

immutable Atoms <: Associative{Symbol, Array}
    indices::Vector{Int}
    props::AtomProps
end

length(atoms::Atoms) = length(atoms.indices)

immutable Atom <: Associative{Symbol, Array}
    index::Int
    props::AtomProps
end

length(atom::Atom) = 1

# Constructors
Frame(props::Pair...) = Frame(FrameProps(props))

Atoms(num_atoms::Int, props::AtomProps) = Atoms(collect(1:num_atoms), props)
Atoms(num_atoms::Int, props::Pair...) = Atoms(num_atoms, AtomProps(props))
Atoms{T <: AbstractVector{Int}}(indices::T, props::AtomProps) = Atoms(collect(indices), props)
Atoms(atoms1::Atoms, atoms2::Atoms) = Atoms(atoms1.indices, atoms2.props)
Atom(atom1::Atom, atoms2::Atoms) = Atom(atom1.index, atoms2.props)
function Atoms(atoms_list::Vector{Atom})
    indices = Int[]
    for item in atoms_list
        push!(indices, item.index)
    end

    # Assumes that all atoms_list.props is same props structure.  Find fast way to check this for sure.
    Atoms(indices, first(atoms_list).props)
end

# Frame accessors
length(f::Frame) = length(f.props)
haskey(f::Frame, key) = haskey(f.props, key)
start(f::Frame) = start(f.props)
next(f::Frame, args...) = next(f.props, args...)
done(f::Frame, args...) = done(f.props, args...)
get(f::Frame, args...) = get(f.props, args...)
setindex!(f::Frame, args...) = setindex!(f.props, args...)

# Atoms accessors
keys(atoms::Atoms) = keys(atoms.props)
haskey(atoms::Atoms, key) = haskey(atoms.props, key)
start(atoms::Atoms) = start(atoms.indices)
function next(atoms::Atoms, state)
    index, state = next(atoms.indices, state)
    Atom(index, atoms.props), state
end
done(atoms::Atoms, state) = done(atoms.indices, state)
eachindex(atoms::Atoms) = 1:length(atoms)
endof(atoms::Atoms) = endof(atoms.indices)
function getindex(atoms::Atoms, index::Int)
    @inbounds i = atoms.indices[index]
    Atom(i, atoms.props)
end
getindex{T <: AbstractVector{Int}}(atoms::Atoms, indices::T) = Atoms(atoms.indices[indices], atoms.props)

function getindex(atoms::Atoms, prop::Symbol)
    i = atoms.indices
    p = atoms.props[prop]

    if ndims(p) == 1
        @inbounds v = view(p, i)
        return v
    elseif ndims(p) == 2
        @inbounds v = view(p, :, i)
        return v
    elseif ndims(p) == 3
        @inbounds v = view(p, :, :, i)
    else
        return 0
    end
end

getindex(atoms::Atoms, atom::Atom) = Atom(atom, atoms)

function setindex!(atoms::Atoms, a::Vector, prop::Symbol)
    atoms.props[prop] = a
end

function merge(atoms::Atoms...)
    num_atoms = sum(length, atoms)
    prop_types = intersect(map(keys, atoms)...)
    props = Dict(
        prop => cat(
            ndims(first(atoms).props[prop]),
            map(a -> a[prop], atoms)...
        ) for prop in prop_types
    )

    Atoms(num_atoms, props)
end

function add_property!{T}(atoms::Atoms, prop::Symbol, default_value::T)
    num_atoms = last(size(first(values(atoms.props))))
    atoms.props[prop] = fill(default_value, num_atoms)
end

# Atom accessors
function getindex(atom::Atom, prop::Symbol)
    i = atom.index
    p = atom.props[prop]

    if ndims(p) == 1
        @inbounds v = p[i]
        return v
    elseif ndims(p) == 2
        @inbounds v = p[:, i]
        return v
    elseif ndims(p) == 3
        @inbounds v = p[:, :, i]
    else
        return 0
    end
end

function setindex!{T}(atom::Atom, value::T, prop::Symbol)
    @inbounds atom.props[prop][atom.index] = value
end

hash(atom::Atom) = hash(object_id(atom.props)) ⊻ hash(atom.index)
hash(atom::Atom, h::UInt) = hash(atom) ⊻ hash(h)
==(atom1::Atom, atom2::Atom) = hash(atom1) == hash(atom2)
isequal(atom1::Atom, atom2::Atom) = atom1 == atom2

# Display

function show(io::IO, frame::Frame)
    props = frame.props
    show_special = (:num_atoms, :num_bonds, :timestep, :atoms, :bonds)
    haskey(props, :atoms) && print(io, length(props[:atoms]), " atoms  ")
    haskey(props, :bonds) && print(io, size(props[:bonds], 2), " bonds  ")
    haskey(props, :timestep) && print(io, "t = ", props[:timestep], "  ")
    for (k, v) in props
        k ∈ show_special || print(io, k, " = ", v, "  ")
    end
end

display(frame::Frame) = show(frame)

function show(io::IO, atoms::Atoms)
    print(io, "$(length(atoms)) atoms")
    print(io, " of types $(tuple(unique(atoms[:type])...))")
end

function show(io::IO, atom::Atom)
    for (k, v) in atom.props
        println(io, k, ": ", atom[k])
    end
end

display(atoms::Atoms) = show(atoms)
display(atom::Atom) = show(atom)

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
# Utilities to filter occurrences of UMLS semantic type in pubmed articles
# Date: September 7, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using BioMedQuery.Entrez.DB

"""
umls_semantic_occurrences(db, umls_semantic_type)

Return a sparse matrix indicating the presence of MESH descriptors associated
with a given umls semantic type in all articles of the input database

###Output

* `des_ind_dict`: Dictionary matching row number to descriptor names
* `disease_occurances` : Sparse matrix. The columns correspond to a feature
vector, where each row is a MESH descriptor. There are as many
columns as articles. The occurance/abscense of a descriptor is labeled as 1/0
"""
function umls_semantic_occurrences(db, umls_semantic_type)

    #retrieve a list of filtered descriptors
    filtered_mesh = Set(filter_mesh_by_concept(db, umls_semantic_type))

    println("-------------------------------------------------------------")
    println("Found ", length(filtered_mesh), " MESH decriptor related to  ", umls_semantic_type)
    println(filtered_mesh)
    println("-------------------------------------------------------------")

    #create a map of filtered descriptor name to index to guarantee order
    des_ind_dict = Dict()

    for (i, fm) in enumerate(filtered_mesh)
        des_ind_dict[fm]= i
    end

    articles = all_pmids(db)

    #create the data-matrix
    disease_occurances = spzeros(length(filtered_mesh), length(articles))

    #Can this process be more efficient using database join/select?
    narticle = 0
    for (i, pmid) in enumerate(articles)

        #get all mesh descriptors associated with give article
        article_mesh = Set(Entrez.DB.get_article_mesh(db, pmid))

        #not all mesh are of the desired semantic type
        article_filtered_mesh = intersect(article_mesh, filtered_mesh)

        #skip if empty
        if isempty(article_filtered_mesh)
            continue
        end

        #otherwise form feature vector for this article
        indices = []
        for d in article_filtered_mesh
            push!(indices, des_ind_dict[d])
        end

        #TO DO: Not sure about the type. Should we choose bool to save space
        # or float to support opperations
        article_dis_feature  = zeros(Int, (length(filtered_mesh),1))
        article_dis_feature[indices] = 1

        #append to data matrix
        disease_occurances[:, i] = article_dis_feature
        narticle+=1
    end

    println("-------------------------------------------------------------")
    println("Found ", narticle, " articles with valid descriptors")
    println("-------------------------------------------------------------")
    return des_ind_dict, disease_occurances

end

# Retrieve all mesh descriptors associated with the given umls_concept
function filter_mesh_by_concept(db, umls_concept)

    uc = string("'", replace(umls_concept, "'", "''") , "'")
    query  = db_query(db, "SELECT mesh FROM mesh2umls
    WHERE umls LIKE $uc ")

    #return data array
    return get_value(query.columns[1])

end


text:
abstract type CoordinateTransform end

# Based on the pseudo-code on page 156 in Loschchilov's PhD Thesis from 2013.
mutable struct AdaptiveEncoding <: CoordinateTransform
  n::Int
  sqrtn::Float64
  initialized::Bool
  c_p::Float64
  c_mu::Float64
  c_1::Float64
  p::Array{Float64, 2}
  C::Array{Float64, 2}
  B::Array{Float64, 2}
  m::Array{Float64, 2}
  m_prev::Array{Float64, 2}

  AdaptiveEncoding(mu, n) = begin
    c_p = 1.0 / sqrt(n)
    c_mu = c_1 = 0.5 / n
    p = zeros(n, 1)
    C = Matrix{Float64}(I, n, n)
    B = Matrix{Float64}(I, n, n)
    new(n, sqrt(n), false, c_p, c_mu, c_1,
      p, C, B, zeros(n, 1), zeros(n, 1))
  end
end

# Update the coordinate transform given a new set of points (x) and their
# utilities/fitnesses (us). Based on
# http://en.wikipedia.org/wiki/CMA-ES#Example_code_in_MATLAB.2FOctave
function update_transform!(ae::AdaptiveEncoding, x::Array{Float64, 2}, us::Array{Float64,2})
  # Init m on the first call.
  if !ae.initialized
    ae.m = x * us
    ae.initialized = true
    # The other variables were already updated in the constructor...
    return ae.B
  end

  ae.m_prev = ae.m
  ae.m = x * us
  Binv = inv(ae.B)

  diff_m = ae.m - ae.m_prev
  z0 = (ae.sqrtn / norm(Binv * diff_m)) * diff_m

  # Update evolution path.
  ae.p = (1.0 - ae.c_p) * ae.p + sqrt(ae.c_p * (2 - ae.c_p)) * z0

  # Skip the rank-mu update when you need to save time.
  # Only do this randomly, at times, then, to save time?
  xis_minus_m_prev = broadcast(-, xi, m_prev)
  zis = (ae.sqrtn / norm(Binv * xis_minus_m_prev)) * xis_minus_m_prev
  C_mu = zeros(ae.n, ae.n)
  # Should be some simpler way to achieve this... Maybe C_mu = zis * diag(us) * zis'
  for i in 1:ae.n
    z = zis[:,i]
    C_mu += (us[i] * (z * z'))
  end

  # Now update the covar matrix. Skip the c_mu parts if not doing the rank-u update...
  ae.C = (1 - ae.c_1 - ae.c_mu) * ae.C + (ae.c_1 * ae.p * ae.p') + ae.c_mu * C_mu

  # Should we enforce symmetry in C?
  # tc = triu(ae.C)
  # ae.C = tc + tc' - Matrix{Float64}(I, ae.C, ae.C) * diag(ae.C)

  # Now B is the square root of the covar matrix, but faster to do eig decomposition.
  eigvalues, eigvectors = eig(ae.C)
  diagD = diag(eigvectors)
  ae.B = eigvalues .* repmat(sqrt(diagD)', ae.n, 1)

  return ae.B
end


text:
function Experiment(::Val{:JuliaRL}, ::Val{:Minimax}, ::Val{:OpenSpiel}, game;)
    env = OpenSpielEnv(string(game))
    agents = MultiAgentManager(
        NamedPolicy(0 => MinimaxPolicy()),
        NamedPolicy(1 => MinimaxPolicy()),
    )
    hooks = MultiAgentHook(0 => TotalRewardPerEpisode(), 1 => TotalRewardPerEpisode())
    description = """
      # Play `$game` in OpenSpiel with Minimax
      """
    Experiment(agents, env, StopAfterEpisode(1), hooks, description)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using UnregisteredWithProject
using Test
using Example

@test f(2) == 2

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using Test
using PotentialCalculation.IdenticalTools

@testset "IdenticalTools" begin
    ii = Identical()
    push!(ii, (2,3))
    @test ! areidentical(ii, (1,2))
    @test areidentical(ii, (3,2))

    jj = Identical(10)
    @test length(jj.identical) == 10

    push!(jj, (3,6,8))
    @test length(jj.identical) == 8

    @test ! areidentical(jj, (2,3))
    @test areidentical(jj, (3,6))
    @test areidentical(jj, (3,8))
    @test areidentical(jj, (3,6,8))
end


text:
using Distributions, MCMCDiagnostics

#NOTE: -----------Assumes IC50s are all in log10 domain-------------
#------------------------------MCMC functions------------------------------
#=
random variables:
latent neut distribution (could ignore this)
change mean
change std
measurement std
latent neut values ~ latent neut distribution
changes ~ change distribution
when a neut value is below threshold, P(value) = integral over implied gaussian for everything below the cutoff.
=#
function MCMC(founder_neuts, variant_neuts, censor_at; prop_bw = 0.02, maxiters = 10000000)

    log_censored_probability(mu,std) = logcdf(Normal(mu,std),censor_at)

    function point_log_likelihoods(points,mus,std; thresh = censor_at)
        LLs = logpdf(Normal(0.0,std),points .- mus)
        for i in 1:length(points)
            if points[i] <= thresh
                LLs[i] = log_censored_probability(mus[i],std)
            end
        end
        return sum(LLs)
    end

    function propose_some(v,std,prop)
        newv = copy(v)
        for i in 1:length(newv)
            if rand()<prop
                newv[i] = v[i] .+ randn()*std
            end
        end
        return newv
    end

    log_change_prior(mu,std) = logpdf(Normal(0.0,10.0),mu) + logpdf(Gamma(2.0,0.1),std)

    #Error STD expected to be fairly small.
    log_error_prior(error_std) = logpdf(Gamma(1.0,0.2),error_std)

    #Could intro hyperparams over this, but it won't make a diff.
    latent_neut_dist(latent_neuts) = sum(logpdf(Normal(2.5,5.0),latent_neuts))

    function log_lik(founder_neuts, variant_neuts, latent_neuts, changes, error_std)
        LL = point_log_likelihoods(founder_neuts,latent_neuts,error_std)
        LL += point_log_likelihoods(variant_neuts,latent_neuts .- changes,error_std)
        return LL
    end

    function log_prior(latent_neuts, changes, change_mean, change_std, error_std)
        LP = latent_neut_dist(latent_neuts)
        LP += sum(logpdf(Normal(change_mean,change_std),changes)) #P(changes|change dist)
        LP += log_change_prior(change_mean,change_std)
        LP += log_error_prior(error_std)
        return LP
    end

    #Slightly_guided_init
    error_std = 0.1
    change_mean = mean((founder_neuts .- variant_neuts)[variant_neuts .> censor_at])
    change_std = std((founder_neuts .- variant_neuts)[variant_neuts .> censor_at])
    latent_neuts = founder_neuts .+ randn()*0.1
    changes = (ones(N) .* change_mean) .+ randn(N)*0.3

    Lpost = log_lik(founder_neuts, variant_neuts, latent_neuts, changes, error_std) + log_prior(latent_neuts, changes, change_mean, change_std, error_std)
    chain = []

    for i in 1:maxiters
        new_error_std = error_std + randn()*prop_bw
        new_change_mean = change_mean + randn()*prop_bw
        new_change_std = change_std + randn()*prop_bw
        new_latent_neuts = propose_some(latent_neuts,prop_bw,0.3)
        new_changes = propose_some(changes,prop_bw,0.3)

        if  new_change_std > 0.05 && new_error_std > 0.05 #Truncating distribution because of stiffness
            new_Lprior = log_prior(new_latent_neuts, new_changes, new_change_mean, new_change_std, new_error_std)
            new_Lpost = new_Lprior + log_lik(founder_neuts, variant_neuts, new_latent_neuts, new_changes, new_error_std)
            if exp(new_Lpost - Lpost) > rand()
                error_std = new_error_std
                change_mean = new_change_mean
                change_std = new_change_std
                latent_neuts = new_latent_neuts
                changes = new_changes
                Lpost = new_Lpost
            end
        end
        if mod(i,2000)==0
            push!(chain,(Lpost,change_mean,change_std,error_std, latent_neuts, changes))
        end
    end
    return chain
end

max_stars_repo_name: MurrellGroup/CensoredNeuts

max_issues_repo_name: MurrellGroup/CensoredNeuts

max_forks_repo_name: MurrellGroup/CensoredNeuts


text:
using Pandas
using Base.Test

df = DataFrame(Dict(:name=>["a", "b"], :age=>[27, 30]))
age = values(df[:age])
age[2] = 31
@test loc(df)[1, "age"] == 31


text:
mutable struct Context
  cache::Union{IdDict{Any,Any},Nothing}
  globals::Union{Dict{GlobalRef,Any},Nothing}
end

Context() = Context(nothing, nothing)

cache(cx::Context) = cx.cache === nothing ? (cx.cache = IdDict()) : cx.cache
globals(cx::Context) = cx.globals === nothing ? (cx.globals = Dict{GlobalRef,Any}()) : cx.globals

struct Pullback{S,T}
  t::T
end

Pullback{S}(x) where S = Pullback{S,typeof(x)}(x)

struct CompileError
  T
  e
end

function Base.showerror(io::IO, e::CompileError)
  print(io, "Compiling $(e.T): ")
  showerror(io, e.e)
end

# interface2.jl

# Wrappers

_forward(f, args...) = _forward(Context(), f, args...)

tailmemaybe(::Nothing) = nothing
tailmemaybe(x::Tuple) = Base.tail(x)

function forward(f, args...)
  y, back = _forward(f, args...)
  y, Δ -> tailmemaybe(back(Δ))
end

function gradient(f, args...)
  y, back = forward(f, args...)
  y isa Real || error("Function output is not scalar")
  return back(Int8(1))
end

derivative(f::F, x) where F = gradient(f, x)[1]

Base.adjoint(f::Function) = x -> derivative(f, x)

# Param-style wrappers

# TODO store ids only
struct Params
  order::Vector{Any}
  params::IdSet{Any}
  Params() = new([], IdSet())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !(x in ps.params)
    push!(ps.order, x)
    push!(ps.params, x)
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)

Params(xs) = push!(Params(), xs...)

function Base.show(io::IO, ps::Params)
  print(io, "Params([")
  join(io, ps.order, ", ")
  print(io, "])")
end

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = print(io, "Grads(...)")

@forward Grads.grads Base.setindex!, Base.getindex, Base.haskey

function forward(f, ps::Params)
  cx = Context()
  y, back = _forward(cx, f)
  y, function (Δ)
    for p in ps
      cache(cx)[p] = nothing
    end
    back(Δ)
    Grads(cx.cache) # TODO make a copy
  end
end


text:
using PyPlot, PyCall

"""
    t = delaunay(p)

Delaunay triangulation `t` of N x 2 node array `p`.
"""
function delaunay(p)
    tri = pyimport("matplotlib.tri")
    t = tri[:Triangulation](p[:,1], p[:,2])
    return Int64.(t[:triangles] .+ 1)
end

"""
    edges, boundary_indices = all_edges(t)

Find all unique edges in the triangulation `t` (ne x 2 array)
Second output is indices to the boundary edges.
"""
function all_edges(t)
    edges = vcat(t[:,[1,2]], t[:,[2,3]], t[:,[3,1]])
    edges = sortslices(sort(edges, dims=2), dims=1)
    dup = all(edges[2:end,:] - edges[1:end-1,:] .== 0, dims=2)[:]
    keep = .![false;dup]
    edges = edges[keep,:]
    dup = [dup;false]
    dup = dup[keep]
    return edges, findall(.!dup)
end

"""
    e = boundary_nodes(t)

Find all boundary nodes in the triangulation `t`.
"""
function boundary_nodes(t)
    edges, boundary_indices = all_edges(t)
    return unique(edges[boundary_indices,:][:])
end

"""
    tplot(p, t, u=nothing)

If `u` == nothing: Plot triangular mesh with nodes `p` and triangles `t`.
If `u` == solution vector: Plot filled contour color plot of solution `u`.
"""
function tplot(p, t, u=nothing)
    clf()
    axis("equal")
    if u == nothing
        tripcolor(p[:,1], p[:,2], t .- 1, 0*t[:,1],
                  cmap="Set3", edgecolors="k", linewidth=1)
    else
        tricontourf(p[:,1], p[:,2], t .- 1, u, 20)
    end
    draw()
end

"""
    inside = inpolygon(p, pv)

Determine if each point in the N x 2 node array `p` is inside the polygon
described by the NE x 2 node array `pv`.
"""
function inpolygon(p::Array{Any,2}, pv::Array{Float64,2})
    path = pyimport("matplotlib.path")
    poly = path[:Path](pv)
    inside = [poly[:contains_point](p[ip,:]) for ip = 1:size(p,1)]
end

function div_poly(polys_xy0,hmax)
    x0=polys_xy0[1,1]
    y0=polys_xy0[1,2]
    ret = [x0 y0]
    for i = 2:size(polys_xy0)[1]
        #println(polys[i,:])
        x=polys_xy0[i,1]
        y=polys_xy0[i,2]
        len=sqrt((x-x0)^2+(y-y0)^2)
        n=ceil(len/hmax)
        #println("len=",len)
        #println("n=",n)
        #println([x0 y0])
        for j=1:n-1
          nx=(x0*(n-j)+x*j)/n
          ny=(y0*(n-j)+y*j)/n
          #println(nx,",",ny)
          ret = [ret; nx ny]
        end
        ret = [ret; x y]
        x0=x
        y0=y
    end
    return ret
end

function centroids(points, tri)
    ret = reshape([],0,2)
    for i = 1:size(tri)[1]
        nx = (points[tri[i,1],1]+points[tri[i,2],1]+points[tri[i,3],1])/3
        ny = (points[tri[i,1],2]+points[tri[i,2],2]+points[tri[i,3],2])/3
        ret = [ret; nx ny]
    end
    return ret
end

function area(points,tri)
   ret = []
   max = 0.0
   maxindex = -1
   for i = 1:size(tri)[1]
       na = (points[tri[i,1],1]*(points[tri[i,2],2]-points[tri[i,3],2])
           + points[tri[i,2],1]*(points[tri[i,3],2]-points[tri[i,1],2])
           + points[tri[i,3],1]*(points[tri[i,1],2]-points[tri[i,2],2]))/2
       ret = [ret; na]
       if na>max
           max=na
           maxindex=i
       end
   end
   return ret, maxindex, max
end

function circumcenter(p, t, it)
    ct = t[it,:]
    dp1 = p[ct[2], :] - p[ct[1], :]
    dp2 = p[ct[3], :] - p[ct[1], :]

    mid1 = (p[ct[2], :] + p[ct[1], :])/2
    mid2 = (p[ct[3], :] + p[ct[1], :])/2

    rhs = mid2-mid1
    s = [-dp1[2] dp2[2] ; dp1[1] -dp2[1]]\rhs
    cpc = mid1 + s[1] * [-dp1[2], dp1[1]]

    return cpc
end

function get_new_tri(inside,t)
    t_new=zeros(Int64,0,3)
    for i = 1:size(t)[1]
        if inside[i]
            t_new = [t_new;t[i,:]']
        end
    end
    return t_new
end

function pmesh(pv, hmax, nref)
    pvd = div_poly(pv,hmax)
    t = delaunay(pvd)

    centers = centroids(pvd,t)
    inside = inpolygon(centers,pv)
    t = get_new_tri(inside,t)
    areas, maxindex, max = area(pvd,t)

    while (max > 0.5*hmax.^2)
        added = circumcenter(pvd,t,maxindex)
        pvd = [pvd; added']
        t = delaunay(pvd)
        centers = centroids(pvd,t)

        inside = inpolygon(centers,pv)
        t = get_new_tri(inside,t)
        areas, maxindex, max = area(pvd,t)
    end

    for i = 1:nref
        edges = all_edges(t)[1]
        for n = 1:size(edges)[1]
            node = (pvd[edges[n,1],:] + pvd[edges[n,2],:])/2.0
            pvd = [pvd;node']
        end
        t = delaunay(pvd)
        centers = centroids(pvd,t)
        inside = inpolygon(centers,pv)
        t = get_new_tri(inside,t)
    end
    e = boundary_nodes(t)

    tplot(pvd,t)
    return pvd,t,e
end


text:
# script that wants to numerically test chatteris diaconis (misspelled with 99% prob)
# for the estimates of beta, fitness,ergm (many names..) in the DirBin1 case
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,KernelDensity,DynNets
using PyCall; pygui(); using PyPlot

Nsample = 100

# load my data
halfPeriod = false
## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

matY_T = YeMidWeekly_T
strIO_T = [sumSq(matY_T,2);sumSq(matY_T,1)]
N,~,T = size(matY_T)

i = 2
nnzStrIO = [ strIO_T[i,strIO_T[i,:].>0] for i =1:2N]
nnzLenghtIO = [length(nnzStrIO[i]) for i=1:2N]
indNnzStrIO = nnzLenghtIO.> 1
acorrNnzStrIO = [autocor(nnzStrIO[indNnzStrIO][i],[1])[1] for i=1:sum(indNnzStrIO)]

plot(nnzLenghtIO,".")
plot(strIO_T[i,strIO_T[i,:].>0])
plot(acorrNnzStrIO)

fracsI = Float64[]
 fracsO = Float64[]
 for t=1:T
    Y = YeMidWeekly_T[:,:,t]
    YfracIn = Y./sum(Y,2)
    YfracOut = Y./sum(Y,1)

    tmpO = YfracIn[.!isnan(YfracIn)][YfracIn[.!isnan(YfracIn)].>0]
    tmpI = YfracOut[.!isnan(YfracOut)][YfracOut[.!isnan(YfracOut)].>0]
    fracsI = [fracsI; tmpI]
    fracsO = [fracsO; tmpO]
end

xGrid = linspace(0,1,10000)
 plt[:hist](fracsI,100)

uniI = unique(fracsI )
countUniI = [sum(uniI[i] .== fracsI) for i=1:length(uniI)]
plot(uniI,countUniI,".")
tabI = sortrows([countUniI uniI ],rev = true)
sum(countUniI[countUniI.>10])./sum(countUniI)
semilogx(tabI[:,1])
prevDegEqFlag = (degsIO_T[:,1:end-1] .== degsIO_T[:,2:end])
prevDegEqFlagI,prevDegEqFlagO = prevDegEqFlag[1:N,:], prevDegEqFlag[1+1N:end,:]
tmp = (sum(prevDegEqFlagI,1) - N).*(sum(prevDegEqFlagO,1) - N)/(N*(N-1))
PyPlot.plot(tmp')

steps = 1
 prevDegDiff = abs.(degsIO_T[:,1:end-steps] - degsIO_T[:,steps+1:end])
 sum(prevDegDiff.<=1)
 degDiffMax = 20
 plt[:hist](prevDegDiff[prevDegDiff.<=degDiffMax],degDiffMax,density = true);grid()
 title("Absolute Values of One Step Degree Variation ")
 ylabel("Fraction of Occurences")
# plot(1:degDiffMax,cumsum(counts(prevDegDiff[:]))[1:degDiffMax]/sum(counts(prevDegDiff[:])))

prevStrDiff = (strIO_T[:,1:end-1] - strIO_T[:,2:end])[:]

prevStrDiffNnz = zeros(1); for i=1:2N  prevStrDiffNnz = vcat(prevStrDiffNnz , diff(strIO_T[i,:][strIO_T[i,:].!=0]) ); end; deleteat!(prevStrDiffNnz,1)
posDiffNnz = prevStrDiffNnz[prevStrDiffNnz.>0]
negDiffNnz = prevStrDiffNnz[prevStrDiffNnz.<0]

 diff(strIO_T[i,:][strIO_T[i,:].!=0])./strIO_T[i,:][strIO_T[i,:].!=0][2:end]

strIO_T

i =12
log(strIO_T[i,:][strIO_T[i,:].!=0][2:end]./strIO_T[i,:][strIO_T[i,:].!=0][1:end-1] )

prevStrRelDiffNnz = zeros(1); for i=1:2N  prevStrRelDiffNnz = vcat(prevStrRelDiffNnz , log(strIO_T[i,:][strIO_T[i,:].!=0][2:end]./strIO_T[i,:][strIO_T[i,:].!=0][1:end-1] )); end;
deleteat!(prevStrRelDiffNnz,1)
plot(prevStrRelDiffNnz[:])
using AverageShiftedHistograms,Plots
plotly()
Plots.plot(ash(prevStrRelDiffNnz),leg=false,size = (1200,700))
xlabel!("log( str_tm1/str_t)")
grid()
posDiff = prevStrDiff[prevStrDiff.>0]
negDiff = prevStrDiff[prevStrDiff.<0]
zeroDiff = prevStrDiff[prevStrDiff.==0]
posecdf = ecdf(posDiff)
negecdf = ecdf(abs.(negDiff))
posecdfNnz = ecdf(posDiffNnz)
negecdfNnz = ecdf(negDiffNnz)

Npoints = 100
tot = length(prevStrDiff)
x = exp(linspace(0,log(maximum(posDiffNnz)),Npoints))
 loglog(x,(1-posecdf(x))*(length(posDiffNnz)/tot))
 grid(which = :both ,alpha = 0.5)
 x = exp(linspace(0,log(maximum(abs.(negDiffNnz))),Npoints))
  loglog(x,(1-negecdf(x))*(length(negDiffNnz)/tot))
  grid(which = :both ,alpha = 0.5)

degDiffMax = 20
 plt[:hist](prevDegDiff[prevDegDiff.<=degDiffMax],degDiffMax,density = true);grid()
 title("Absolute Values of One Step Degree Variation ")
 ylabel("Fraction of Occurences")

prevDegEqFlagI,prevDegEqFlagO = prevDegEqFlag[1:N,:], prevDegEqFlag[1+1N:end,:]

tmp = (sum(prevDegEqFlagI,1) - N).*(sum(prevDegEqFlagO,1) - N)/(N*(N-1))
plot(tmp')

##
using Utilities
plot(mean(degsIO_T,2),autocor(degsIO_T',[1],demean=true)',".")

nnzDegs_T  = degsIO_T[sumSq(degsIO_T.>0,2).>T/10,:]

plot(sum(prevDegDiff.>0,2),autocor(degsIO_T',[20 ],demean=true)',".")
 grid()
autocor(randi(0,30,297),[1])
#
lagsVec = [1]
maxDeg = 20
 T = 297
 Nsample = 100
 sampleAutocorr = zeros(Nsample,2)
 for i=1:Nsample
    sampleAutocorr[i,:] = [autocor([randi(maxDeg) for i=1:T ],lagsVec) autocor(maxDeg*rand(T) ,lagsVec)]
 end
 plot(sampleAutocorr)

tmp = ones(100)
 tmp[1 ] = 12
 autocor(tmp,lagsVec)

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using MacroTools

struct Effect{E}
  effect::E
  cont
end

effect(e) = error("No handler available for effect $(e).")

cps(k, ::typeof(effect), e) = Effect(e, k)

function handle(f, T, h)
  r = cps(identity, f)
  while r isa Effect && r.effect isa T
    r = h(r.effect, x -> handle(() -> r.cont(x), T, h))
  end
  return r
end

macro effect(ex)
  @capture(ex, begin x__; ((e_::T_|e_), k_) -> h_ end) ||
    error("@effect begin ...; (e, k) -> ...; end")
  quote
    local f = () -> $(esc.(x)...)
    local h = ($(esc(e)), $(esc(k))) -> $(esc(h))
    handle(f, $(esc(T)), h, )
  end
end


text:
using Base.Test
using DataFrames, NLreg

const sd1 = readtable(Pkg.dir("NLreg","data","sd1.csv.gz"))
sd1[:ID] = pool(sd1[:ID])

# nl = gnfit(NonlinearLS(Logsd1(:(CONC ~ TIME), sd1)),true)

pl = fit([LogTr] * BolusSD1(CONC ~ TIME, sd1),true)

@test_approx_eq coef(pl) [1.1429558452268844,-1.4036989045425874]

pl = fit(AsympReg([2.:7],[18.6,22.6,25.1,27.2,29.1,30.1]),true)

const bod = DataFrame(Time=[1.:5,7],Demand=[8.3,10.3,19.0,16.0,15.6,19.8])

pl = fit(AsympOrig(Demand ~ Time,bod),true)


text:
@benchset "abstract_gp" begin
    for N in [10,]
        f = GP(CustomMean(sin), EQ(), GPC())
        x = randn(N)
        fx = f(x)
        y = rand(fx)
        @benchset "N=$N" begin
            create_benchmarks("logpdf EQ $N (x, y)", true, (x, y)->logpdf(f(x), y), x, y)
            create_benchmarks("rand EQ $N", true, x->rand(f(x)), x)
            create_benchmarks("cov(f(x)) $N", true, x->cov(f(x)), x)
        end
    end
end


text:
# domain.jl
# Definition of the abstract Domain type and its interface

spaceof(::Domain{T}) where {T} = spacetype(T)

eltype(::Type{Domain{T}}) where {T} = T
eltype(::Type{D}) where {D <: Domain} = eltype(supertype(D))

dimension(::Type{Domain{T}}) where {T} = dimension_type(T)
dimension(::Type{D}) where {D <: Domain} = dimension(supertype(D))
dimension(d::Domain) = dimension(typeof(d))
dimension_type(::Type{SVector{N,T}}) where {N,T} = N
dimension_type(::Type{T}) where {T <: Number} = 1

"""
If the type `T` is a container type, the elements of `T` may have a different
`subeltype`. If `T` is not a container, `subeltype` simply evaluates to `T`.
"""
subeltype(d::Domain) = subeltype(spaceof(d))
subeltype(::Type{T}) where {T} = subeltype(GSpace{T})

"A `EuclideanDomain` is any domain whose eltype is `SVector{N,T}`."
const EuclideanDomain{N,T} = Domain{SVector{N,T}}

# We implement the indicator function by overriding `in`.
# The implementation of `in` at the level of Domain converts the point into
# the element type of the domain using promotion, and calls `indomain`.
# Concrete subtypes should implement `indomain`, rather than `in`.
_in(x::T, d::Domain{T}) where {T} = indomain(x, d)

function _in(x, d::Domain)
   T = promote_type(typeof(x), eltype(d))
   T == Any && return false
   in(convert(T, x), convert(Domain{T}, d))
end

in(x, d::Domain) = _in(x, d)

# The user may supply a vector. We attempt to convert it to the right space.
in(x::Vector{T}, d::Domain) where {T} = in(convert(eltype(d), x), d)

"""
Return a suitable tolerance to use for verifying whether a point is close to
a domain. Typically, the tolerance is close to the precision limit of the numeric
type associated with the domain.
"""
default_tolerance(d::Domain) = default_tolerance(d, subeltype(d))
default_tolerance(d::Domain, ::Type{T}) where {T <: Real} = 100eps(T)
default_tolerance(d::Domain, ::Type{Complex{T}}) where {T <: Real} = 100eps(T)
# Default tolerance for integers and for any other type is zero
default_tolerance(d::Domain, ::Type{T}) where {T <: Integer} = zero(T)
default_tolerance(d::Domain, ::Type{T}) where {T} = zero(T)

approx_in(x::T, d::Domain{T}, tolerance = default_tolerance(d)) where {T} =
   approx_indomain(x, d, tolerance)

function approx_in(x, d::Domain, tolerance = default_tolerance(d))
   T = promote_type(typeof(x), eltype(d))
   T == Any && return false
   approx_in(convert(T, x), convert(Domain{T}, d), tolerance)
end

isapprox(d1::Domain, d2::Domain; kwds...) = d1 == d2

isreal(d::Domain) = isreal(spaceof(d))

infimum(d::Domain) = minimum(d)  # if the minimum exists, then it is also the infimum
supremum(d::Domain) = maximum(d)  # if the maximum exists, then it is also the supremum

# override minimum and maximum for closed sets
function boundary end
const ∂ = boundary


text:
@generated function struct_as_namedtuple(st)
  A = (Expr(:(=), n, :(st.$n)) for n in setdiff(fieldnames(st),(:kwargs,)))
  Expr(:tuple, A...)
end

Base.@pure remaker_of(prob::T) where {T} = parameterless_type(T){isinplace(prob)}

# Define `remaker_of` for the types that does not (make sense to)
# implement `isinplace` trait:
for T in [
    NoiseProblem,
    SplitFunction,  # TODO: use isinplace path for type-stability
    TwoPointBVPFunction,
    ]
  @eval remaker_of(::$T) = $T
end

"""
    remake(thing; <keyword arguments>)

Re-construct `thing` with new field values specified by the keyword
arguments.
"""
function remake(thing; kwargs...)
  T = remaker_of(thing)
  if :kwargs ∈ fieldnames(typeof(thing))
    T(; struct_as_namedtuple(thing)...,thing.kwargs...,kwargs...)
  else
    T(; struct_as_namedtuple(thing)...,kwargs...)
  end
end

isrecompile(prob::ODEProblem{iip}) where {iip} = (prob.f isa ODEFunction) ? !(typeof(prob.f.f) <: FunctionWrapper) : true

function remake(thing::ODEProblem; kwargs...)
  T = remaker_of(thing)
  tup = merge(merge(struct_as_namedtuple(thing),thing.kwargs),kwargs)
  if !isrecompile(thing)
    if isinplace(thing)
      f = wrapfun_iip(unwrap_fw(tup.f.f),(tup.u0,tup.u0,tup.p,tup.tspan[1]))
    else
      f = wrapfun_oop(unwrap_fw(tup.f.f),(tup.u0,tup.p,tup.tspan[1]))
    end
    tup2 = (f = convert(ODEFunction{isinplace(thing)},f),)
    tup = merge(tup, tup2)
  end
  T(; tup...)
end

function remake(thing::AbstractJumpProblem; kwargs...)
  parameterless_type(thing)(remake(thing.prob;kwargs...))
end


text:
export iauPr00
"""
Precession-rate part of the IAU 2000 precession-nutation models
(part of MHB2000).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2    double  TT as a 2-part Julian Date (Note 1)

Returned:
   dpsipr,depspr  double  precession corrections (Notes 2,3)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The precession adjustments are expressed as "nutation
   components", corrections in longitude and obliquity with respect
   to the J2000.0 equinox and ecliptic.

   3. Although the precession adjustments are stated to be with respect
   to Lieske et al. (1977), the MHB2000 model does not specify which
   set of Euler angles are to be used and how the adjustments are to
   be applied.  The most literal and straightforward procedure is to
   adopt the 4-rotation epsilon_0, psi_A, omega_A, xi_A option, and
   to add dpsipr to psi_A and depspr to both omega_A and eps_A.

   4. This is an implementation of one aspect of the IAU 2000A nutation
   model, formally adopted by the IAU General Assembly in 2000,
   namely MHB2000 (Mathews et al. 2002).

References:

   Lieske, J.H., Lederle, T., Fricke, W. & Morando, B., "Expressions
   for the precession quantities based upon the IAU (1976) System of
   Astronomical Constants", Astron.Astrophys., 58, 1-16 (1977)

   Mathews, P.M., Herring, T.A., Buffet, B.A., "Modeling of nutation
   and precession   New nutation series for nonrigid Earth and
   insights into the Earth's interior", J.Geophys.Res., 107, B4,
   2002.  The MHB2000 code itself was obtained on 9th September 2002
   from ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A.

   Wallace, P.T., "Software for Implementing the IAU 2000
   Resolutions", in IERS Workshop 5.1 (2002).

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPr00(double date1, double date2, double *dpsipr, double *depspr)

function iauPr00(date1::Real, date2::Real)
   ref_dpsipr = Ref{Float64}(0.0)
   ref_depspr = Ref{Float64}(0.0)

   ccall((:iauPr00, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
         convert(Float64, date1),
         convert(Float64, date2),
         ref_dpsipr, ref_depspr)

   return ref_dpsipr[], ref_depspr[]
end


text:
# to ensure that propertynames that aren't fieldnames are always viewed as
# "defined", in MLJType objects. See #115.
function _isdefined(object, name)
    pnames = propertynames(object)
    fnames = fieldnames(typeof(object))
    name in pnames && !(name in fnames) && return true
    isdefined(object, name)
end

function _equal_to_depth_one(x1, x2)
    names = propertynames(x1)
    names === propertynames(x2) || return false
    for name in names
        getproperty(x1, name) == getproperty(x2, name) || return false
    end
    return true
end

@doc """
    deep_properties(::Type{<:MLJType})

Given an `MLJType` subtype `M`, the value of this trait should be a
tuple of any properties of `M` to be regarded as "deep".

When two instances of type `M` are to be tested for equality, in the
sense of `==` or `is_same_except`, then the values of a "deep"
property (whose values are assumed to be of composite type) are deemed
to agree if all corresponding properties *of those property values*
are `==`.

Any property of `M` whose values are themselves of `MLJType` are
"deep" automatically, and should not be included in the trait return
value.

See also [`is_same_except`](@ref)

### Example

Consider an `MLJType` subtype `Foo`, with a single field of
type `Bar` which is *not* a subtype of `MLJType`:

```julia
mutable struct Bar
    x::Int
end

mutable struct Foo <: MLJType
    bar::Bar
end
```

Then the mutability of `Foo` implies `Foo(1) != Foo(1)` and so, by the
definition `==` for `MLJType` objects (see [`is_same_except`](@ref))
we have

```julia
Bar(Foo(1)) != Bar(Foo(1))
```

However after the declaration

```julia
MLJModelInterface.deep_properties(::Type{<:Foo}) = (:bar,)
```
We have

```julia
Bar(Foo(1)) == Bar(Foo(1))
```
"""
StatisticalTraits.deep_properties

"""
    is_same_except(m1, m2, exceptions::Symbol...; deep_properties=Symbol[])

If both `m1` and `m2` are of `MLJType`, return `true` if the
following conditions all hold, and `false` otherwise:

- `typeof(m1) === typeof(m2)`

- `propertynames(m1) === propertynames(m2)`

- with the exception of properties listed as `exceptions` or bound to
  an `AbstractRNG`, each pair of corresponding property values is
  either "equal" or both undefined. (If a property appears as a
  `propertyname` but not a `fieldname`, it is deemed as always defined.)

The meaining of "equal" depends on the type of the property value:

- values that are themselves of `MLJType` are "equal" if they are
equal in the sense of `is_same_except` with no exceptions.

- values that are not of `MLJType` are "equal" if they are `==`.

In the special case of a "deep" property, "equal" has a different
meaning; see [`deep_properties`](@ref)) for details.

If `m1` or `m2` are not `MLJType` objects, then return `==(m1, m2)`.

"""
is_same_except(x1, x2) = ==(x1, x2)

function is_same_except(m1::M1,
    m2::M2,
    exceptions::Symbol...
) where {M1<:MLJType, M2<:MLJType}
    typeof(m1) === typeof(m2) || return false
    names = propertynames(m1)
    propertynames(m2) === names || return false

    for name in names
        if !(name in exceptions)
            if !_isdefined(m1, name)
               !_isdefined(m2, name) || return false
            elseif _isdefined(m2, name)
                if name in deep_properties(M1)
                    _equal_to_depth_one(
                        getproperty(m1,name),
                        getproperty(m2, name)
                    ) || return false
                else
                    (
                        is_same_except(
                            getproperty(m1, name),
                            getproperty(m2, name)
                        ) ||
                        getproperty(m1, name) isa AbstractRNG ||
                        getproperty(m2, name) isa AbstractRNG
                    ) || return false
                end
            else
                return false
            end
        end
    end
    return true
end

==(m1::M1, m2::M2) where {M1<:MLJType, M2<:MLJType} = is_same_except(m1, m2)

# for using `replace` or `replace!` on collections of MLJType objects
# (eg, Model objects in a learning network) we need a stricter
# equality and a corresponding definition of `in`.
Base.isequal(m1::MLJType, m2::MLJType) = (m1 === m2)

# Note: To prevent julia crash, it seems we need to annotate the type
# of itr:
function special_in(x, itr)::Union{Bool,Missing}
    for y in itr
        ismissing(y) && return missing
        y === x && return true
    end
    return false
end

Base.in(x::MLJType, itr::Set) = special_in(x, itr)
Base.in(x::MLJType, itr::AbstractVector) = special_in(x, itr)
Base.in(x::MLJType, itr::Tuple) = special_in(x, itr)

# A version of `in` that actually uses `==`:

"""
    isrepresented(object::MLJType, objects)

Test if `object` has a representative in the iterable
`objects`. This is a weaker requirement than `object in objects`.

Here we say `m1` *respresents* `m2` if `is_same_except(m1, m2)` is
`true`.

"""
isrepresented(object::MLJType, ::Nothing) = false

function isrepresented(object::MLJType, itr)::Union{Bool,Missing}
    for m in itr
        ismissing(m) && return missing
        is_same_except(m, object) && return true
    end
    return false
end


text:
using PatModules
@once module subpackage
    include("one.jl")
    include("two.jl")

    import .one: say_hi
    import .two: say_bye
end


text:
"""
    eom_linear3d_nonlinearcoupling(x, p, n) -> Function

Equations of motion for a 3d linear system with nonlinear coupling [1].
The difference equations are

```math
\\begin{aligned}
x_1(t+1) &= a_1 x_1 (1-x_1(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{1}(t) \\
x_2(t+1) &= a_1 x_2 (1-x_2(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{2}(t) + b x_1 x_2 \\
x_3(t+1) &= a_3 x_3 (1-x_3(t))^2  e^{-x_3(t)^2} + 0.4 \\xi_{3}(t) + c x_{2}(t) \\
         &+ d x_{1}(t)^2.
\\end{aligned}
```

Here, ``\\xi_{1,2,3}(t)`` are independent normally distributed noise processes,
representing dynamical noise in the system, with zero mean and standard
deviations ``\\sigma_1``, ``\\sigma_2``, ``\\sigma_3``, respectively.

# References
1. Gourévitch, B., Le Bouquin-Jeannès, R., & Faucon, G. (2006). Linear and
nonlinear causality between signals: methods, examples and neurophysiological
applications. Biological Cybernetics, 95(4), 349–369.
"""
function eom_linear3d_nonlinearcoupling(x, p, n)
    x₁, x₂, x₃ = (x...,)
    a₁, a₂, a₃, b, c, d, σ₁, σ₂, σ₃ = (p...,)
    ξ₁ = rand(Normal(0, σ₁))
    ξ₂ = rand(Normal(0, σ₂))
    ξ₃ = rand(Normal(0, σ₃))

    dx₁ = a₁*x₁*(1-x₁)^2 * exp(-x₁^2) + 0.4*ξ₁
    dx₂ = a₂*x₂*(1-x₂)^2 * exp(-x₂^2) + 0.4*ξ₂ + b*x₁*x₂
    dx₃ = a₃*x₃*(1-x₃)^2 * exp(-x₃^2) + 0.4*ξ₃ + c*x₂ + d*x₁^2

    return SVector{3}(dx₁, dx₂, dx₃)
end

function linear3d_nonlinearcoupling(uᵢ, a₁, a₂, a₃, b, c, d, σ₁, σ₂, σ₃)
    p = [a₁, a₂, a₃, b, c, d, σ₁, σ₂, σ₃]
    s = DiscreteDynamicalSystem(eom_linear3d_nonlinearcoupling, uᵢ, p)
    return s
end

"""
    linear3d_nonlinearcoupling(;uᵢ = rand(3), σ₁ = 1.0, σ₂ = 1.0, σ₃ = 1.0,
        a₁ = 3.4, a₂ = 3.4, a₃ = 3.4, b = 0.5, c = 0.3, d = 0.5) -> DiscreteDynamicalSystem

A 3d linear system with nonlinear coupling [1]. The difference equations are

```math
\begin{aligned}
x_1(t+1) &= a_1 x_1 (1-x_1(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{1}(t) \\
x_2(t+1) &= a_1 x_2 (1-x_2(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{2}(t) + b x_1 x_2 \\
x_3(t+1) &= a_3 x_3 (1-x_3(t))^2  e^{-x_3(t)^2} + 0.4 \\xi_{3}(t) + c x_{2}(t) \\
         &+ d x_{1}(t)^2.
\end{aligned}
```

Here, ``\\xi_{1,2,3}(t)`` are independent normally distributed noise processes,
representing dynamical noise in the system, with zero mean and standard
deviations ``\\sigma_1``, ``\\sigma_2``, ``\\sigma_3``, respectively.

# References
1. Gourévitch, B., Le Bouquin-Jeannès, R., & Faucon, G. (2006). Linear and
nonlinear causality between signals: methods, examples and neurophysiological
applications. Biological Cybernetics, 95(4), 349–369.
"""
linear3d_nonlinearcoupling(;uᵢ = rand(3), σ₁ = 1.0, σ₂ = 1.0, σ₃ = 1.0,
        a₁ = 3.4, a₂ = 3.4, a₃ = 3.4, b = 0.5, c = 0.3, d = 0.5) =
    linear3d_nonlinearcoupling(uᵢ, a₁, a₂, a₃, b, c, d, σ₁, σ₂, σ₃)


text:
@testset "Include All Files Function Tests" begin

  @test isdefined(Julz, :include_all_files) == true

  cd("dummy") do

    Julz.generate("function", "baz")

    @test isdefined(:baz) != true

    Julz.include_all_files("$(pwd())/src/functions", is_revised=false)

    @test isdefined(:baz) == true

    Julz.destroy("function", "baz")

  end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using Pkg
using BinaryBuilder

name = "Libtask"
version = v"0.3.1"
commit_id = "fbe338053f402d76524d0a01f3796dd7da90b781"

# see https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/336
# ENV["CI_COMMIT_TAG"] = ENV["TRAVIS_TAG"] = "v" * string(version)

sources = [
    "https://github.com/TuringLang/Libtask.jl.git" => commit_id,
]

# Bash recipe for building across all platforms
script_file = joinpath(@__DIR__, "build_dylib.sh")
if !isfile(script_file) # when run generate_buildjl.jl
    script_file = joinpath(@__DIR__, "L/Libtask/build_dylib.sh")
end
script = read(script_file, String)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    # Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:aarch64),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libtask_v1_0", :libtask_v1_0)
    LibraryProduct("libtask_v1_1", :libtask_v1_1)
    LibraryProduct("libtask_v1_2", :libtask_v1_2)
    LibraryProduct("libtask_v1_3", :libtask_v1_3)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
# build_file = "products/build_$(name).v$(version).jl"
build_tarballs(ARGS, name, version, sources,
               script, platforms, products, dependencies)

max_stars_repo_name: jbrea/Yggdrasil

max_issues_repo_name: jbrea/Yggdrasil

max_forks_repo_name: jbrea/Yggdrasil


text:
module MCIntegrals
using Requires

include("rng.jl")
include("core.jl")
include("vegas.jl")
include("cuba.jl")
include("plots.jl")

function __init__()
    @require HCubature="19dc6840-f33b-545b-b366-655c7e3ffd49" include("hcubature.jl")
end

end # module


text:
padfac(::DFTAlg) = 1

# We don't pad a DFT since it is already correct
padimage(::DFTAlg, img) = img

function DFTAlg(u::AbstractArray, v::AbstractArray)
    @argcheck length(u) == length(v)
    uv = Matrix{eltype(u)}(undef, 2, length(u))
    uv[1,:] .= u
    uv[2,:] .= v
    return ObservedNUFT(DFTAlg(), uv)
end

function plan_nuft(alg::ObservedNUFT{<:DFTAlg}, img, dx, dy)
    uv = alg.uv
    xitr, yitr = imagepixels(img)
    dft = Array{Complex{eltype(uv)},3}(undef, size(uv,2), size(img)...)
    @fastmath for i in axes(img,2), j in axes(img,1), k in axes(uv,2)
        u = uv[1,k]
        v = uv[2,k]
        dft[k, j, i] = cispi(-2(u*xitr[i] + v*yitr[j]))
    end
    # reshape to a matrix so we can take advantage of an easy BLAS call
    return reshape(dft, size(uv,2), :)
end

function make_phases(alg::ObservedNUFT{<:DFTAlg}, img)
    u = @view alg.uv[1,:]
    v = @view alg.uv[2,:]
    dx, dy = pixelsizes(img)
    visibility.(Ref(img.pulse), u*dx, v*dy).*dx.*dy
end

@inline function create_cache(alg::ObservedNUFT{<:DFTAlg}, plan, phases, img)
    return NUFTCache(alg, plan, phases, reshape(img.im, :))
end


text:
################################################################################
# Control Input
################################################################################

function input_impulse!(joint::Translational{T},
    pbody::Node, cbody::Node,
    timestep::T, clear::Bool) where T

    xa, qa = current_configuration(pbody.state)
    xb, qb = current_configuration(cbody.state)

    input = joint.input
    Ta = impulse_transform(:parent, joint, xa, qa, xb, qb)
    Tb = impulse_transform(:child, joint, xa, qa, xb, qb)
    Faw = Ta[1:3,1:3] * input
    τaa = Ta[4:6,1:3] * input
    Fbw = Tb[1:3,1:3] * input
    τbb = Tb[4:6,1:3] * input

    pbody.state.F2 += Faw
    pbody.state.τ2 += τaa/2
    cbody.state.F2 += Fbw
    cbody.state.τ2 += τbb/2

    clear && (joint.input = szeros(T,3))
    return
end

################################################################################
# Control Jacobian
################################################################################

function input_jacobian_control(relative::Symbol,
    joint::Translational,
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion) where T

    Ta = impulse_transform(relative, joint, xa, qa, xb, qb)
    X = Ta[1:3,1:3]
    Q = 0.5 * Ta[4:6,1:3]
    return [X; Q]
end

function input_jacobian_configuration(relative::Symbol,
    joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion) where T

    # d[Faw;2τaa]/d[xa,qa]
    ∇aa = impulse_transform_jacobian(:parent, relative, joint, xa, qa, xb, qb, joint.input)
    FaXa = ∇aa[1:3, 1:3]
    FaQa = ∇aa[1:3, 4:6]
    τaXa = 0.5 * ∇aa[4:6, 1:3]
    τaQa = 0.5 * ∇aa[4:6, 4:6]

    # d[Fbw;2τbb]/d[xa,qa]
    ∇ba = impulse_transform_jacobian(:child, relative, joint, xa, qa, xb, qb, joint.input)
    FbXa = ∇ba[1:3,1:3]
    FbQa = ∇ba[1:3,4:6]
    τbXa = 0.5 * ∇ba[4:6,1:3]
    τbQa = 0.5 * ∇ba[4:6,4:6]

    return FaXa, FaQa, τaXa, τaQa, FbXa, FbQa, τbXa, τbQa
end


text:
#API for calling/creating services. Syntax is practically identical to rospy.
export Service, ServiceProxy, wait_for_service

"""
    ServiceProxy{T}(name; kwargs...)
    ServiceProxy(name, T; kwargs...)

Create a proxy object used to invoke a remote service. Use `srv_proxy(msg_request)` with the object
to invoke the service call. Keyword arguments are directly passed to rospy.
"""
struct ServiceProxy{SrvType <: AbstractService}
    o::PyObject

    function ServiceProxy{ST}(name::AbstractString; kwargs...) where ST <: AbstractService
        @debug("Creating <$ST> service proxy for '$name'")
        rospycls = _get_rospy_class(ST)
        new{ST}(__rospy__.ServiceProxy(ascii(name), rospycls; kwargs...))
    end
end

function ServiceProxy(name::AbstractString, srv::Type{ST}; kwargs...) where ST <: AbstractService
    ServiceProxy{ST}(ascii(name); kwargs...)
end

function (srv::ServiceProxy{ST})(req::AbstractSrv) where ST <: AbstractService
    if ! isa(req, _srv_reqtype(ST))
        throw(ArgumentError(
            string("Incorrect service request type: ", typeof(req),
                   ", expected: ", _srv_reqtype(ST))))
    end
    pyresp = pycall(srv.o, PyObject, convert(PyObject, req))
    resp = convert(_srv_resptype(ST), pyresp)
    resp
end

"""
    Service{T}(name, callback; kwargs...)
    Service(name, T, callback; kwargs...)

Create a service object that can receive requests and provide responses. The callback can be of
any callable type. Keyword arguments are directly passed to rospy.
"""
mutable struct Service{SrvType <: AbstractService}
    handler
    srv_obj::PyObject
    cb_interface::PyObject
    async_loop::Task

    function Service{ST}(name::AbstractString, handler; kwargs...) where ST <: AbstractService
        @debug("Providing <$ST> service at '$name'")
        rospycls = _get_rospy_class(ST)

        cond = Base.AsyncCondition()
        pysrv = _py_ros_callbacks."ServiceCallback"(CB_NOTIFY_PTR[], cond.handle)

        srvobj = try
            __rospy__.Service(ascii(name), rospycls, pysrv."srv_cb"; kwargs...)
        catch err
            if isa(err, PyCall.PyError)
                error("Problem during service creation: $(err.val.args[1])")
            else
                rethrow(err)
            end
        end

        rosobj = new{ST}(handler, srvobj, pysrv)

        cbloop = Task(() -> _callback_async_loop(rosobj, cond))
        schedule(cbloop)

        rosobj.async_loop = cbloop
        return rosobj
    end
end

function Service(name::AbstractString, srv::Type{ST}, handler; kwargs...) where ST <: AbstractService
    Service{ST}(ascii(name), handler; kwargs...)
end

"""
    wait_for_service(srv_name; kwargs...)

Block until the specified service is available. Keyword arguments are directly passed to rospy.
Throws an exception if the waiting timeout period is exceeded.
"""
function wait_for_service(service::AbstractString; kwargs...)
    try
        __rospy__.wait_for_service(ascii(service); kwargs...)
    catch ex
        error("Timeout exceeded waiting for service '$service'")
    end
end


text:
using Printf

Base.show(io::IO, b::SeawaterBuoyancy{FT}) where FT =
    println(io, "SeawaterBuoyancy{$FT}: g = $(b.gravitational_acceleration)", '\n',
                "└── equation of state: $(b.equation_of_state)")

Base.show(io::IO, eos::LinearEquationOfState{FT}) where FT =
    println(io, "LinearEquationOfState{$FT}: ", @sprintf("α = %.2e, β = %.2e", eos.α, eos.β))


text:
using Dates
import Base: isless

raw"""
    DatetimeRotatingFileLogger(dir, file_pattern; always_flush=true, rotation_callback=identity)
    DatetimeRotatingFileLogger(f::Function, dir, file_pattern; always_flush=true, rotation_callback=identity)

Construct a `DatetimeRotatingFileLogger` that rotates its file based on the current date.
The constructor takes a log output directory, `dir`, and a filename pattern.
The filename pattern given is interpreted through the `Dates.format()` string formatter,
allowing for yearly all the way down to minute-level log rotation.  Note that if you
wish to have a filename portion that is not interpreted as a format string, you may need
to escape portions of the filename, as shown in the example below.

It is possible to pass a formatter function as the first argument to control the output.
The formatting function should be of the form `f(io::IOContext, log_args::NamedTuple)`
where `log_args` has the following fields:
`(level, message, _module, group, id, file, line, kwargs)`.
See `?LoggingExtra.handle_message_args` for more information about what each field represents.

It is also possible to pass `rotation_callback::Function` as a keyword argument. This function
will be called every time a file rotation is happening. The function should accept one
argument which is the absolute path to the just-rotated file. The logger will block until
the callback function returns. Use `@async` if the callback is expensive.

# Examples

```julia
# Logger that logs to a new file every day
logger = DatetimeRotatingFileLogger(log_dir, raw"\a\c\c\e\s\s-yyyy-mm-dd.\l\o\g")

# Logger with a formatter function that rotates the log file hourly
logger = DatetimeRotatingFileLogger(log_dir, raw"yyyy-mm-dd-HH.\l\o\g") do io, args
    println(io, args.level, " | ", args.message)
end

# Example callback function to compress the recently-closed file
compressor(file) = run(`gzip $(file)`)
logger = DatetimeRotatingFileLogger(...; rotation_callback=compressor)
```
"""
mutable struct DatetimeRotatingFileLogger <: AbstractLogger
    logger::Union{SimpleLogger,FormatLogger}
    dir::String
    filename_pattern::DateFormat
    next_reopen_check::DateTime
    always_flush::Bool
    reopen_lock::ReentrantLock
    current_file::Union{String,Nothing}
    rotation_callback::Function
end

function DatetimeRotatingFileLogger(dir, filename_pattern; always_flush=true, rotation_callback=identity)
    DatetimeRotatingFileLogger(nothing, dir, filename_pattern; always_flush=always_flush, rotation_callback=rotation_callback)
end
function DatetimeRotatingFileLogger(f::Union{Function,Nothing}, dir, filename_pattern; always_flush=true, rotation_callback=identity)
    # Construct the backing logger with a temp IOBuffer that will be replaced
    # by the correct filestream in the call to reopen! below
    logger = if f === nothing
        SimpleLogger(IOBuffer(), BelowMinLevel)
    else # f isa Function
        FormatLogger(f, IOBuffer(); always_flush=false) # no need to flush twice
    end
    # abspath in case user constructs the logger with a relative path and later cd's.
    drfl = DatetimeRotatingFileLogger(logger, abspath(dir),
        DateFormat(filename_pattern), now(), always_flush, ReentrantLock(), nothing, rotation_callback)
    reopen!(drfl)
    return drfl
end

similar_logger(::SimpleLogger, io) = SimpleLogger(io, BelowMinLevel)
similar_logger(l::FormatLogger, io) = FormatLogger(l.f, io, l.always_flush)
function reopen!(drfl::DatetimeRotatingFileLogger)
    if drfl.current_file !== nothing
        # close the old IOStream and pass the file to the callback
        close(drfl.logger.stream)
        drfl.rotation_callback(drfl.current_file)
    end
    new_file = calc_logpath(drfl.dir, drfl.filename_pattern)
    drfl.current_file = new_file
    io = open(new_file, "a")
    drfl.logger = similar_logger(drfl.logger, io)
    drfl.next_reopen_check = next_datetime_transition(drfl.filename_pattern)
    return nothing
end

"""
    next_datetime_transition(fmt::DateFormat)

Given a DateFormat that is being applied to our filename, what is the next
time at which our filepath will need to change?
"""
function next_datetime_transition(fmt::DateFormat)
    extract_token(x::Dates.DatePart{T}) where {T} = T
    token_timescales = Dict(
        # Milliseconds is the smallest timescale
        's' => Millisecond(1),
        # Seconds
        'S' => Second(1),
        # Minutes
        'M' => Minute(1),
        # Hours
        'I' => Hour(1),
        'H' => Hour(1),
        # Days
        'd' => Day(1),
        'e' => Day(1),
        'E' => Day(1),
        # Month
        'm' => Month(1),
        'u' => Month(1),
        'U' => Month(1),
        # Year
        'y' => Year(1),
        'Y' => Year(1),
    )

    # Dates for some reason explicitly does not define equality between the smaller
    # timescales (Second, Minute, Day, etc..) and the larger, non-constant timescales
    # (Month, Year).  We do so explicitly here, without committing type piracy:
    custom_isless(x, y) = isless(x, y)
    custom_isless(x::Union{Millisecond,Second,Minute,Hour,Day}, y::Union{Month, Year}) = true
    custom_isless(x::Union{Month, Year}, y::Union{Millisecond,Second,Minute,Hour,Day}) = false

    tokens = filter(t -> isa(t, Dates.DatePart), collect(fmt.tokens))
    minimum_timescale = first(sort(map(t -> token_timescales[extract_token(t)], tokens), lt=custom_isless))
    if minimum_timescale < Minute(1)
        throw(ArgumentError("rotating the logger with sub-minute resolution not supported"))
    end
    return ceil(now(), minimum_timescale)
end

calc_logpath(dir, filename_pattern) = joinpath(dir, Dates.format(now(), filename_pattern))

function handle_message(drfl::DatetimeRotatingFileLogger, args...; kwargs...)
    lock(drfl.reopen_lock) do
        if now() >= drfl.next_reopen_check
            reopen!(drfl)
        end
    end
    handle_message(drfl.logger, args...; kwargs...)
    drfl.always_flush && flush(drfl.logger.stream)
end

shouldlog(drfl::DatetimeRotatingFileLogger, arg...) = true
min_enabled_level(drfl::DatetimeRotatingFileLogger) = BelowMinLevel
catch_exceptions(drfl::DatetimeRotatingFileLogger) = catch_exceptions(drfl.logger)