text:
using Test
using QuantumLattices.Prerequisites.Traits

import QuantumLattices.Prerequisites.Traits: parameternames, isparameterbound, contentnames, getcontent, contenttype, dissolve

abstract type FT{T} end
@inline parameternames(::Type{<:FT}) = (:content,)
@inline isparameterbound(::Type{<:FT}, ::Val{:content}, D) = !isconcretetype(D)
@inline contentnames(::Type{<:FT}) = (:content,)
@inline dissolve(m::FT, ::Val{:content}, f::Function, args::Tuple, kwargs::NamedTuple) = f(getcontent(m, :content), args...; kwargs...)

struct SameNameField <: FT{Vector{Int}}
    info::String
    content::Vector{Int}
end
@inline parameternames(::Type{SameNameField}) = ()
@inline contentnames(::Type{SameNameField}) = (:info, :content)

struct DiffNameField <: FT{Vector{Int}}
    info::String
    table::Vector{Int}
end
@inline parameternames(::Type{DiffNameField}) = ()
@inline contentnames(::Type{DiffNameField}) = (:info, :content)
@inline getcontent(m::DiffNameField, ::Val{:content}) = getfield(m, :table)
@inline contenttype(DF::Type{DiffNameField}, ::Val{:content}) = fieldtype(DF, :table)

@testset "commontype" begin
    @test commontype(Int) == Int
    @test commontype(Union{Int, Float64}) == Real
    @test commontype(Union{Int, Float64, Complex{Float64}}) == Number

    fx = x::Int->(x==1 ? 1 : x==2 ? 2.0 : 3.0im)
    @test commontype(fx, Tuple{Int}) == Number
end

@testset "abstracttypehelper" begin
    @test DataType(Int) == Int
    T = Vector{<:Integer}
    @test DataType(T) == T.body

    @test supertype(Vector, :Array) == Vector
    @test supertype(Vector, :AbstractArray) == AbstractVector

    @test parametercount(Vector) == 2
    @test parametertype(Vector{<:Real}, 1) == Real
    @test parametertype(Vector{<:Real}, 2) == 1
    @test parametertypes(Vector{<:Real}) == Tuple{Real, 1}
    @test isparameterbound(Vector, 1, Real) == false
    @test isparameterbound(Vector, 2, 1) == false
    @test isparameterbounds(Vector{<:Real}, Tuple{Real, 1}) == (false, false)
    @test reparameter(Vector, 1, Real, true) == Vector{<:Real}
    @test reparameter(Vector, 1, Real, false) ==  Vector{Real}
    @test reparameter(Vector{Int}, 2, 3, false) == Array{Int, 3}

    @test contentnames(Vector) == ()
    @test rawtype(Int) == Int
    @test rawtype(Vector{Int}) == Array
    @test rawtype(Vector) == Array
    @test rawtype(Vector{<:Real}) == Array
    @test fulltype(Array, Tuple{Int, 1}, (false, false)) == Array{Int, 1}
    @test fulltype(Array, Tuple{Real, 2}, (true, false)) == Array{<:Real, 2}

    @test promoteparameters(NamedTuple{(:a, :b), Tuple{Int, Float64}}, NamedTuple{(:b, :c), Tuple{Complex{Float64}, Int}}) == NamedTuple{(:a, :b, :c), Tuple{Int, Complex{Float64}, Int}}

    @test parametercount(FT) == 1
    @test parametername(FT, 1) == :content
    @test parameterorder(FT, :content) == 1
    @test parametertype(FT{Vector}, 1) == parametertype(FT{<:Vector}, :content) == Vector
    @test parametertype(FT{Vector{Int}}, 1) == parametertype(FT{Vector{Int}}, :content) == Vector{Int}
    @test parameterpair(FT{Vector}, :content) == parameterpair(FT{Vector}, 1) == Pair{:content, Vector}
    @test isparameterbound(FT, :nocotent, Vector) == false
    @test isparameterbound(FT, :content, Vector) == true
    @test isparameterbound(FT, :content, Vector{Int}) == false
    @test hasparameter(FT, :content) == true

    @test parameternames(FT) == (:content,)
    @test parametertypes(FT{Vector}) == Tuple{Vector}
    @test parameterpairs(FT{Vector}) == NamedTuple{(:content,), Tuple{Vector}}
    @test isparameterbounds(FT, NamedTuple{(:content,), Tuple{Vector}}) == (true,)
    @test isparameterbounds(FT, NamedTuple{(:content,), Tuple{Vector{Int}}}) == (false,)
    @test reparameter(FT{Vector{Int}}, :content, Vector) == FT{<:Vector}
    @test reparameter(FT{Vector{Int}}, :content, Vector{Float64}) == FT{Vector{Float64}}

    @test fulltype(FT, NamedTuple{(:content,), Tuple{Vector}}) == FT{<:Vector}
    @test fulltype(FT, NamedTuple{(:content,), Tuple{Vector{Int}}}) == FT{Vector{Int}}

    @test contentnames(FT) == (:content,)
    @test contentname(FT, 1) == :content
    @test contentorder(FT, :content) == 1
    @test contentcount(FT) == 1
    @test hascontent(FT, :content) == true

    @test contentnames(SameNameField) == (:info, :content)
    @test hasparameter(SameNameField, :content) == false
    @test contenttype(SameNameField, :info) == contenttype(SameNameField, Val(:info)) == String
    @test contenttype(SameNameField, :content) == contenttype(SameNameField, Val(:content)) == Vector{Int}
    @test contenttypes(SameNameField) == Tuple{String, Vector{Int}}

    a = SameNameField("abcd", [1, 2, 3, 4])
    @test getcontent(a, :content) == getcontent(a, 2) == getcontent(a, Val(:content)) == [1, 2, 3, 4]
    @test dissolve(a, getindex, (1,)) == ("abcd", 1)

    @test contentnames(DiffNameField) == (:info, :content)
    @test hasparameter(DiffNameField, :content) == false
    @test contenttype(DiffNameField, :info) == contenttype(DiffNameField, Val(:info)) == String
    @test contenttype(DiffNameField, :content) == contenttype(DiffNameField, Val(:content)) == Vector{Int}
    @test contenttypes(DiffNameField) == Tuple{String, Vector{Int}}

    b = DiffNameField("abcd", [1, 2, 3, 4])
    @test getcontent(b, :content) == getcontent(b, 2) == getcontent(b, Val(:content)) == [1, 2, 3, 4]
    @test dissolve(b, getindex, (2,)) == ("abcd", 2)
end

struct EFO{F1, F2, F3}
    f1::F1
    f2::F2
    f3::F3
end
Base.:(==)(fc1::EFO, fc2::EFO) = ==(efficientoperations, fc1, fc2)
Base.isequal(fc1::EFO, fc2::EFO) = isequal(efficientoperations, fc1, fc2)
Base.:<(fc1::EFO, fc2::EFO) = <(efficientoperations, fc1, fc2)
Base.isless(fc1::EFO, fc2::EFO) = isless(efficientoperations, fc1, fc2)
Base.isapprox(fc1::EFO, fc2::EFO; atol::Real=10^-5, rtol::Real=10^-5) = isapprox(efficientoperations, (:f1, :f2), fc1, fc2; atol=atol, rtol=rtol)
Base.replace(fc::EFO; kwargs...) = replace(efficientoperations, fc; kwargs...)

@testset "efficientoperations" begin
    @test ==(efficientoperations, (), ())
    @test isequal(efficientoperations, (), ())
    @test ==(efficientoperations, (1, 2), (1, 2, 3)) == false
    @test isequal(efficientoperations, (1, 2), (1, 2, 3)) == false
    @test isapprox(efficientoperations, (1.0, 2.0), (1, 2)) == true
    @test isapprox(efficientoperations, (), (), ()) == true
    @test isapprox(efficientoperations, (), (), (1,)) == false

    fc1, fc2 = EFO(1.0, 2, 3), EFO(1, 2, 3)
    @test fc1 == fc2
    @test isequal(fc1, fc2)
    fc1, fc2 = EFO(1.0, 2, 3), EFO(1, 2, 4.0)
    @test fc1 < fc2
    @test isless(fc1, fc2)

    fc1, fc2 = EFO(1.0+10^-6, 2, 3), EFO(1, 2, 3)
    @test fc1 ≈ fc2
    fc1, fc2 = EFO(1.0, 2-10^-6, 3), EFO(1, 2, 3)
    @test fc1 ≈ fc2
    fc1, fc2 = EFO(1.0, 2, 3-10^-6), EFO(1, 2, 3)
    @test fc1 ≉  fc2

    @test replace(EFO(1, 2, 3), f1='c') == EFO('c', 2, 3)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PORTA_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PORTA")
JLLWrappers.@generate_main_file("PORTA", UUID("c3fa2e09-48e0-5371-872a-ed3ac32dd1fc"))
end  # module PORTA_jll

ext: jl


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

const LineNum = Int

mutable struct InferenceState
    params::InferenceParams
    result::InferenceResult # remember where to put the result
    linfo::MethodInstance
    sptypes::Vector{Any}    # types of static parameter
    slottypes::Vector{Any}
    mod::Module
    currpc::LineNum

    # info on the state of inference and the linfo
    src::CodeInfo
    world::UInt
    valid_worlds::WorldRange
    nargs::Int
    stmt_types::Vector{Any}
    stmt_edges::Vector{Any}
    stmt_info::Vector{Any}
    # return type
    bestguess #::Type
    # current active instruction pointers
    ip::BitSet
    pc´´::LineNum
    nstmts::Int
    # current exception handler info
    cur_hand #::Union{Nothing, Pair{LineNum, prev_handler}}
    handler_at::Vector{Any}
    n_handlers::Int
    # ssavalue sparsity and restart info
    ssavalue_uses::Vector{BitSet}
    throw_blocks::BitSet

    cycle_backedges::Vector{Tuple{InferenceState, LineNum}} # call-graph backedges connecting from callee to caller
    callers_in_cycle::Vector{InferenceState}
    parent::Union{Nothing, InferenceState}

    # TODO: move these to InferenceResult / Params?
    cached::Bool
    limited::Bool
    inferred::Bool
    dont_work_on_me::Bool

    # The place to look up methods while working on this function.
    # In particular, we cache method lookup results for the same function to
    # fast path repeated queries.
    method_table::CachedMethodTable{InternalMethodTable}

    # The interpreter that created this inference state. Not looked at by
    # NativeInterpreter. But other interpreters may use this to detect cycles
    interp::AbstractInterpreter

    # src is assumed to be a newly-allocated CodeInfo, that can be modified in-place to contain intermediate results
    function InferenceState(result::InferenceResult, src::CodeInfo,
                            cached::Bool, interp::AbstractInterpreter)
        linfo = result.linfo
        code = src.code::Array{Any,1}
        toplevel = !isa(linfo.def, Method)

        sp = sptypes_from_meth_instance(linfo::MethodInstance)

        nssavalues = src.ssavaluetypes::Int
        src.ssavaluetypes = Any[ NOT_FOUND for i = 1:nssavalues ]
        stmt_info = Any[ nothing for i = 1:length(code) ]

        n = length(code)
        s_edges = Any[ nothing for i = 1:n ]
        s_types = Any[ nothing for i = 1:n ]

        # initial types
        nslots = length(src.slotflags)
        argtypes = result.argtypes
        nargs = length(argtypes)
        s_argtypes = VarTable(undef, nslots)
        slottypes = Vector{Any}(undef, nslots)
        for i in 1:nslots
            at = (i > nargs) ? Bottom : argtypes[i]
            s_argtypes[i] = VarState(at, i > nargs)
            slottypes[i] = at
        end
        s_types[1] = s_argtypes

        ssavalue_uses = find_ssavalue_uses(code, nssavalues)
        throw_blocks = find_throw_blocks(code)

        # exception handlers
        cur_hand = nothing
        handler_at = Any[ nothing for i=1:n ]
        n_handlers = 0

        W = BitSet()
        push!(W, 1) #initial pc to visit

        if !toplevel
            meth = linfo.def
            inmodule = meth.module
        else
            inmodule = linfo.def::Module
        end

        valid_worlds = WorldRange(src.min_world,
            src.max_world == typemax(UInt) ? get_world_counter() : src.max_world)
        frame = new(
            InferenceParams(interp), result, linfo,
            sp, slottypes, inmodule, 0,
            src, get_world_counter(interp), valid_worlds,
            nargs, s_types, s_edges, stmt_info,
            Union{}, W, 1, n,
            cur_hand, handler_at, n_handlers,
            ssavalue_uses, throw_blocks,
            Vector{Tuple{InferenceState,LineNum}}(), # cycle_backedges
            Vector{InferenceState}(), # callers_in_cycle
            #=parent=#nothing,
            cached, false, false, false,
            CachedMethodTable(method_table(interp)),
            interp)
        result.result = frame
        cached && push!(get_inference_cache(interp), result)
        return frame
    end
end

method_table(interp::AbstractInterpreter, sv::InferenceState) = sv.method_table

function InferenceState(result::InferenceResult, cached::Bool, interp::AbstractInterpreter)
    # prepare an InferenceState object for inferring lambda
    src = retrieve_code_info(result.linfo)
    src === nothing && return nothing
    validate_code_in_debug_mode(result.linfo, src, "lowered")
    return InferenceState(result, src, cached, interp)
end

function sptypes_from_meth_instance(linfo::MethodInstance)
    toplevel = !isa(linfo.def, Method)
    if !toplevel && isempty(linfo.sparam_vals) && isa(linfo.def.sig, UnionAll)
        # linfo is unspecialized
        sp = Any[]
        sig = linfo.def.sig
        while isa(sig, UnionAll)
            push!(sp, sig.var)
            sig = sig.body
        end
    else
        sp = collect(Any, linfo.sparam_vals)
    end
    for i = 1:length(sp)
        v = sp[i]
        if v isa TypeVar
            fromArg = 0
            # if this parameter came from arg::Type{T}, then `arg` is more precise than
            # Type{T} where lb<:T<:ub
            sig = linfo.def.sig
            temp = sig
            for j = 1:i-1
                temp = temp.body
            end
            Pi = temp.var
            while temp isa UnionAll
                temp = temp.body
            end
            sigtypes = temp.parameters
            for j = 1:length(sigtypes)
                tj = sigtypes[j]
                if isType(tj) && tj.parameters[1] === Pi
                    fromArg = j
                    break
                end
            end
            if fromArg > 0
                ty = fieldtype(linfo.specTypes, fromArg)
            else
                ub = v.ub
                while ub isa TypeVar
                    ub = ub.ub
                end
                if has_free_typevars(ub)
                    ub = Any
                end
                lb = v.lb
                while lb isa TypeVar
                    lb = lb.lb
                end
                if has_free_typevars(lb)
                    lb = Bottom
                end
                if Any <: ub && lb <: Bottom
                    ty = Any
                else
                    tv = TypeVar(v.name, lb, ub)
                    ty = UnionAll(tv, Type{tv})
                end
            end
        elseif isa(v, Core.TypeofVararg)
            ty = Int
        else
            ty = Const(v)
        end
        sp[i] = ty
    end
    return sp
end

_topmod(sv::InferenceState) = _topmod(sv.mod)

# work towards converging the valid age range for sv
function update_valid_age!(sv::InferenceState, worlds::WorldRange)
    sv.valid_worlds = intersect(worlds, sv.valid_worlds)
    @assert(sv.world in sv.valid_worlds, "invalid age range update")
    nothing
end

update_valid_age!(edge::InferenceState, sv::InferenceState) = update_valid_age!(sv, edge.valid_worlds)

function record_ssa_assign(ssa_id::Int, @nospecialize(new), frame::InferenceState)
    old = frame.src.ssavaluetypes[ssa_id]
    if old === NOT_FOUND || !(new ⊑ old)
        # typically, we expect that old ⊑ new (that output information only
        # gets less precise with worse input information), but to actually
        # guarantee convergence we need to use tmerge here to ensure that is true
        frame.src.ssavaluetypes[ssa_id] = old === NOT_FOUND ? new : tmerge(old, new)
        W = frame.ip
        s = frame.stmt_types
        for r in frame.ssavalue_uses[ssa_id]
            if s[r] !== nothing # s[r] === nothing => unreached statement
                if r < frame.pc´´
                    frame.pc´´ = r
                end
                push!(W, r)
            end
        end
    end
    nothing
end

function add_cycle_backedge!(frame::InferenceState, caller::InferenceState, currpc::Int)
    update_valid_age!(frame, caller)
    backedge = (caller, currpc)
    contains_is(frame.cycle_backedges, backedge) || push!(frame.cycle_backedges, backedge)
    add_backedge!(frame.linfo, caller)
    return frame
end

# temporarily accumulate our edges to later add as backedges in the callee
function add_backedge!(li::MethodInstance, caller::InferenceState)
    isa(caller.linfo.def, Method) || return # don't add backedges to toplevel exprs
    if caller.stmt_edges[caller.currpc] === nothing
        caller.stmt_edges[caller.currpc] = []
    end
    push!(caller.stmt_edges[caller.currpc], li)
    nothing
end

# used to temporarily accumulate our no method errors to later add as backedges in the callee method table
function add_mt_backedge!(mt::Core.MethodTable, @nospecialize(typ), caller::InferenceState)
    isa(caller.linfo.def, Method) || return # don't add backedges to toplevel exprs
    if caller.stmt_edges[caller.currpc] === nothing
        caller.stmt_edges[caller.currpc] = []
    end
    push!(caller.stmt_edges[caller.currpc], mt)
    push!(caller.stmt_edges[caller.currpc], typ)
    nothing
end

function poison_callstack(infstate::InferenceState, topmost::InferenceState, poison_topmost::Bool)
    poison_topmost && (topmost = topmost.parent)
    while !(infstate === topmost)
        if call_result_unused(infstate)
            # If we won't propagate the result any further (since it's typically unused),
            # it's OK that we keep and cache the "limited" result in the parents
            # (non-typically, this means that we lose the ability to detect a guaranteed StackOverflow in some cases)
            # TODO: we might be able to halt progress much more strongly here,
            # since now we know we won't be able to keep anything much that we learned.
            # We were mainly only here to compute the calling convention return type,
            # but in most situations now, we are unlikely to be able to use that information.
            break
        end
        infstate.limited = true
        for infstate_cycle in infstate.callers_in_cycle
            infstate_cycle.limited = true
        end
        infstate = infstate.parent
        infstate === nothing && return
    end
end

function print_callstack(sv::InferenceState)
    while sv !== nothing
        print(sv.linfo)
        sv.limited && print("  [limited]")
        !sv.cached && print("  [uncached]")
        println()
        for cycle in sv.callers_in_cycle
            print(' ', cycle.linfo)
            cycle.limited && print("  [limited]")
            println()
        end
        sv = sv.parent
    end
end


text:
"""
Represents the API to interact with, either an actual QBO instance, a sandbox,
or a mock API for testing purposes
"""
abstract type API end

struct WebAPI <: API
    oauth::HttpCommon.URI
    endpoint::HttpCommon.URI
end
api_uri(api::WebAPI, path) = HttpCommon.URI(api.endpoint, path = string(api.endpoint.path,path))

const DEFAULT_API = WebAPI(HttpCommon.URI("https://appcenter.intuit.com/connect/oauth2"),
                            HttpCommon.URI("https://sandbox-quickbooks.api.intuit.com/v3"))

using Base.Meta
"""
For a method taking an API argument, add a new method without the API argument
that just calls the method with DEFAULT_API.
"""
macro api_default(func)
    call = func.args[1]
    has_kwargs = isexpr(call.args[2], :parameters)
    newcall = Expr(:call, call.args[1], (has_kwargs ?
        [Expr(:parameters, Expr(:..., :kwargs)); call.args[4:end]] : call.args[3:end])...)
    argnames = map(has_kwargs ? call.args[4:end] : call.args[3:end]) do expr
        isexpr(expr, :kw) && (expr = expr.args[1])
        isexpr(expr, Symbol("::")) && return expr.args[1]
        @assert isa(expr, Symbol)
        return expr
    end
    esc(Expr(:toplevel, func,
        Expr(:function, newcall, Expr(:block,
            has_kwargs ? :($(call.args[1])(DEFAULT_API, $(argnames...); kwargs...)) :
                            :($(call.args[1])(DEFAULT_API, $(argnames...)))
        ))))
end

abstract type Authorization end
abstract type QBObject end

struct AuthCode <: Authorization
    code::String
    realm::Int
end

struct OAuth2 <: Authorization
    token::String
end

struct AnonymousAuth <: Authorization
end

function authenticate_headers!(headers, auth::OAuth2)
    headers["Authorization"] = "Bearer $(auth.token)"
    return headers
end


text:
@inline function getVelocityDelta(joint::Rotational0, body1::AbstractBody, body2::Body{T}, ω::SVector{3,T}) where T
    Δω = ω # in body1 frame
    return Δω
end

@inline function getPositionDelta(joint::Rotational0, body1::AbstractBody, body2::Body{T}, θ::Quaternion{T}) where T
    Δq = θ # in body1 frame
    return Δq
end

@inline function setForce!(joint::Rotational0, body1::Body, body2::Body{T}, τ::SVector{3,T}, No) where T
    τ1 = vrotate(-τ, body1.q[No] * joint.qoff)
    τ2 = -τ1

    body1.τ[No] = τ1
    body2.τ[No] = τ2
    return
end

@inline function setForce!(joint::Rotational0, body1::Origin, body2::Body{T}, τ::SVector{3,T}, No) where T
    body2.τ[No] = vrotate(τ, joint.qoff)
    return
end

@inline function minimalCoordinates(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, No)
    body2.q[No]
end

@inline g(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = g(joint)

@inline ∂g∂posa(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, No) = ∂g∂posa(joint)
@inline ∂g∂posb(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, No) = ∂g∂posb(joint)
@inline ∂g∂vela(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = ∂g∂vela(joint)
@inline ∂g∂velb(joint::Rotational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = ∂g∂velb(joint)


text:
using Distributed
using SharedArrays
using GenGlobal
using GenGlobal.testModule

# -----------------------------------------
# demo on global variables and GenGlobal
# -----------------------------------------

# GenGlobal.testModule defines a global "globalx"

# set a global variable's value in another module
set_globalx(3)

# globalx not available in Main
try
    @show globalx
catch e
    println(e)
end

# but it is through module
@show GenGlobal.testModule.globalx

# Can't set global normally
try
    GenGlobal.testModule.globalx = 2
catch e
    @show e
end

# so use the generated functions
set_globalx(2)
@show get_globalx()

# -----------------------------------------
# Parallel computing & GenGlobal
# -----------------------------------------

# initialize workers
pids = addprocs()

# initialize a shared array visible to each worker in pids
ncols = 20
bigsa = SharedArray{Float64}(10,ncols)

# tell globals about the modules
@everywhere begin
    using GenGlobal
    using GenGlobal.testModule
end

# set globals in each worker
# Can set up temp arrays on remotes this way
@eval @everywhere begin
    set_globalx(3)
    set_g_sa($bigsa)
    set_g_ysame2(zeros(2,2))
    set_g_ydiff2(zeros(2,2))
end

# matrix of zeros
@show bigsa

# Update large shared array slice by slice
# sa[i,j] = i
s = @sync @distributed for j in 1:ncols
    update_shared!(j)  # sets      g_sa[:,j] = 1:end
end
@show fetch(s)
@show bigsa

# Compute things based on the shared array
# For example, can compute likelihood this way, and possibly gradient
s = @sync @distributed (pplus) for j in 1:ncols
    compute_shared(j)  # for each j, compute logsumexp(g_sa[:,j]) and add to ysame2. Add myid() to ydiff2
end
@show s

# Show what the different matrices are
@everywhere @show get_y2()

# add them up!
# Example: have each worker compute part of likelihood + gradient, then add up partial gradients on each worker
remote_mapreduce(get_y2, pplus)


text:
using SemanticReleaseExample
using Test

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


text:
using HigherOrderClustering
using Base.Test
using Combinatorics

function test1()
    I = [1, 1, 1, 2, 2, 3]
    J = [2, 3, 4, 3, 4, 4]
    A = convert(SparseMatrixCSC{Int64,Int64},
                sparse(I, J, ones(length(I)), 4, 4))
    A = max.(A, A')

    @test kcliques(A, 2) == [3, 3, 3, 3]
    @test kcliques(A, 3) == [3, 3, 3, 3]
    @test kcliques(A, 4) == [1, 1, 1, 1]
    @test kcliques(A, 5) == [0, 0, 0, 0]

    ccf2 = higher_order_ccfs(A, 2)
    @test ccf2.global_hoccf == 1.0
    @test ccf2.avg_hoccf    == 1.0
    @test ccf2.avg_hoccf2   == 1.0

    ccf3 = higher_order_ccfs(A, 3)
    @test ccf3.global_hoccf == 1.0
    @test ccf3.avg_hoccf    == 1.0
    @test ccf3.avg_hoccf2   == 1.0

    push!(I, 3, 4, 4)
    push!(J, 7, 7, 4)
    A = convert(SparseMatrixCSC{Int64,Int64},
                sparse(I, J, ones(length(I)), 7, 7))
    A = max.(A, A')

    @test sum(kcliques(A, 2)) == 8 * 2
    @test sum(kcliques(A, 3)) == 5 * 3
    @test sum(kcliques(A, 4)) == 1 * 4
    @test sum(kcliques(A, 5)) == 0 * 5

    @test kcliques(A, 4) == [1, 1, 1, 1, 0, 0, 0]
end

function test2()
    # Test k-cliques on a small network
    A = load_example_data("celegans.txt")
    degs = vec(sum(A, 2))
    n = size(A, 2)
    deg_order = zeros(Int64, n)
    deg_order[sortperm(vec(sum(A, 1)))] = collect(1:n)

    function higher_neighbors(j::Int64)
        oj = deg_order[j]
        nbrs = find(A[:, j])
        return filter(x -> deg_order[x] > oj, nbrs)
    end

    counts = zeros(Int64, size(A, 2))
    for j = 1:n
        for w in higher_neighbors(j)
            counts[[j, w]] += 1
        end
    end
    counts == kcliques(A, 2)

    counts = zeros(Int64, size(A, 2))
    for j = 1:n
        for (w, x) in combinations(higher_neighbors(j), 2)
            if A[w, x] > 0
                counts[[j, w, x]] += 1
            end
        end
    end
    @test counts == kcliques(A, 3)

    counts = zeros(Int64, size(A, 2))
    for j = 1:n
        for (w, x, y) in combinations(higher_neighbors(j), 3)
            if A[w, x] > 0 && A[w, y] > 0 && A[x, y] > 0
                counts[[j, w, x, y]] += 1
            end
        end
    end
    @test counts == kcliques(A, 4)

    counts = zeros(Int64, size(A, 2))
    for j = 1:n
        for (w, x, y, z) in combinations(higher_neighbors(j), 4)
            if A[w, x] > 0 && A[w, y] > 0 && A[x, y] > 0 &&
                A[w, z] > 0 && A[x, z] > 0 && A[y, z] > 0
                counts[[j, w, x, y, z]] += 1
            end
        end
    end
    @test counts == kcliques(A, 5)
end

test1()
test2()


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

using Markdown
using InteractiveUtils

# ╔═╡ 1da5ad0e-828f-11eb-3a21-1734f2d7bb36
using TestImages,FileIO,Images,ImageView,MosaicViews

# ╔═╡ b0079f7e-82c4-11eb-0540-ed1b691768e5
using Plots,ImageContrastAdjustment

# ╔═╡ fed405b6-828d-11eb-2d86-977a78c747ee
md"""
# Basic Image Import,Show,Processing and Export
This example shows how to read an image into the workspace, adjust the contrast in the image, and then write the adjusted image to a file
"""

# ╔═╡ c5e83452-8283-11eb-21e5-4d3d301bf70a
begin
	# Read Image
img1 = testimage("chelsea.jpg"); #one way
img2 = load("assets/reference.jpg"); #second way
end

# ╔═╡ d80d7608-8290-11eb-0dcf-53ee481f7193
begin
	imshow(img1); #imported through ImageView.jl
	mosaicview(img1,img2;nrow=1) #imported through MosaiViews.jl
end

# ╔═╡ d3334628-82be-11eb-37b4-a730795b4a8c
begin
 # Apparently no proper documentation on how to find details on a image
	x=typeof(img1);
	y=size(img1);
	z= sizeof(img1);
end

# ╔═╡ 4556d368-82c2-11eb-271e-916cba11f094
x

# ╔═╡ 4b7c750e-82c2-11eb-3846-ff49451b80f5
y

# ╔═╡ 4e4014e4-82c2-11eb-3fc1-971c351e14bb
z

# ╔═╡ 708364d8-82c3-11eb-29ab-619106b39d5a
val=build_histogram(img1)

# ╔═╡ fbed91ae-82c5-11eb-089d-eb27ae7a0a8c
histogram(val[2][0:255])

# ╔═╡ a885205e-82c5-11eb-2a12-17cf529de2eb
img_equalized = adjust_histogram(img1, Equalization(nbins = 256))

# ╔═╡ 2c4f2c18-82c6-11eb-272b-d5c0169ed2c8
begin
	#method to save a image with a particular directory in a specific format froma a loaded image
	save("assets/cat.png",img1)
	save("assets/img_equalised.png",img_equalized)
end

# ╔═╡ Cell order:
# ╠═fed405b6-828d-11eb-2d86-977a78c747ee
# ╠═1da5ad0e-828f-11eb-3a21-1734f2d7bb36
# ╠═c5e83452-8283-11eb-21e5-4d3d301bf70a
# ╠═d80d7608-8290-11eb-0dcf-53ee481f7193
# ╠═d3334628-82be-11eb-37b4-a730795b4a8c
# ╠═4556d368-82c2-11eb-271e-916cba11f094
# ╠═4b7c750e-82c2-11eb-3846-ff49451b80f5
# ╠═4e4014e4-82c2-11eb-3fc1-971c351e14bb
# ╠═b0079f7e-82c4-11eb-0540-ed1b691768e5
# ╠═708364d8-82c3-11eb-29ab-619106b39d5a
# ╠═fbed91ae-82c5-11eb-089d-eb27ae7a0a8c
# ╠═a885205e-82c5-11eb-2a12-17cf529de2eb
# ╠═2c4f2c18-82c6-11eb-272b-d5c0169ed2c8

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
using Test, LinearAlgebra
# using Revise
using Plots, BinningAnalysis, BinningAnalysisPlots, Distributions
# pyplot()
gr()

function generate_correlated_data(n)
    n_thermalize = floor(Int, 0.03*n)
    A = Base.clamp.(0.3 .+ rand(n_thermalize,n_thermalize), 0, 1);
    A = 0.5*(A+A');
    issymmetric(A)
    A = A + n_thermalize*I;
    @assert isposdef(A)
    C = cholesky(A);
    y_corr = C.L'*rand(n_thermalize)

    i = floor(Int, n_thermalize*0.9)
    μ = mean(y_corr[i:end])
    σ = std(y_corr[i:end])
    y = vcat(y_corr, rand(Normal(μ, σ), n-n_thermalize))
end

n = 10_000
data = generate_correlated_data(n)
x = FullBinner(data)

plot(x)
savefig(joinpath(@__DIR__, "plot.png"))

histogram(x)
savefig(joinpath(@__DIR__, "histogram.png"))

binning(x)
savefig(joinpath(@__DIR__, "binning.png"))

corrplot(x)
savefig(joinpath(@__DIR__, "corrplot.png"))

tauplot(x)
savefig(joinpath(@__DIR__, "tauplot.png"))

x = LogBinner(data)

binning(x)

tauplot(x)

corrplot(x)


text:
export readeventlist, writeeventlist

"""
    readeventlist(fd::IO; <keyword args>)
    readeventlist(file::AbstractString; <keyword args>) -> SparseMatrixCSC{Events}, node list

Read list of events from file returns undirected dynamic network represented as sparse matrix.
An event is an interaction between two nodes from start time to stop time.

# Keyword arguments
- `symmetric=true` : If false, only fill the top right triangle, else make resulting matrix symmetric.
- `header=false` : If true, ignore first line.
- `sortby=nothing`: If not `nothing`, sort nodes w.r.t this function (`by` argument in `sort`).
- `format=:timefirst` : If `format=:timefirst`, each line has format (`start_time stop_time node1 node2`).
    If `format=:nodefirst`, each line has format (`node1 node2 start_time stop_time`).
"""
function readeventlist(fd::IO; symmetric=true, header=false,
                       sortby=nothing, format=:timefirst)
    nodes = Set{String}()
    events = Vector{Tuple{String,String,Float64,Float64}}()
    if header
        readline(fd)
    end
    while !eof(fd)
        line = readskipping(fd)
        if isempty(line); continue; end
        vals = split(strip(line))
        if length(vals)==4
            if format==:nodefirst
                vals = vals[[3,4,1,2]]
            end
            t_s = parse(Float64,vals[1])
            t_e = parse(Float64,vals[2])
            if t_s > t_e error("Event time start > end") end
            vals[3]==vals[4] && continue
            push!(events,(vals[3],vals[4],t_s,t_e))
            push!(nodes,vals[3])
            push!(nodes,vals[4])
        else
            error("Format error: $line")
        end
    end
    nodes = collect(nodes)
    if sortby!=nothing
        nodes = sort(nodes,by=sortby)
    end
    nodesnum = indexmap(nodes)
    edges = Matrix{Int}(length(events),2)
    timestamps = Vector{Events}(length(events))
    for i = 1:length(events)
        vals = events[i]
        n1 = nodesnum[vals[1]]
        n2 = nodesnum[vals[2]]
        edges[i,1] = min(n1,n2)
        edges[i,2] = max(n1,n2)
        timestamps[i] = Events([(vals[3],vals[4])])
    end
    events2dynet(edges[:,1], edges[:,2], timestamps,
                 length(nodes), nodes,
                 symmetric=symmetric)
end
readeventlist(file::AbstractString; args...) =
    open(x -> readeventlist(x;args...), file, "r")

"""
    writeeventlist(fd::IO, dy::DynamicNetwork; <keyword args>)
    writeeventlist(file::AbstractString, dy::DynamicNetwork; <keyword args>)

Write list of events to file from undirected dynamic network represented as sparse matrix.
An event is an interaction between two nodes from start time to stop time.

# Keyword arguments
- `format=:timefirst` : If `format=:timefirst`, each line has format (`start_time stop_time node1 node2`).
    If `format=:nodefirst`, each line has format (`node1 node2 start_time stop_time`).
"""
function writeeventlist(fd::IO, dy::DynamicNetwork; format=:timefirst)
    G = dy.G
    nodes = dy.nodes
    m,n = size(G)
    for col = 1:n
        for j = nzrange(G,col)
            row = G.rowval[j]
            val = G.nzval[j]
            for x in val.timestamps
                if format==:nodefirst
                    println(fd, nodes[row], " ", nodes[col], " ", x[1], " ", x[2])
                else
                    println(fd, x[1], " ", x[2], " ", nodes[row], " ", nodes[col])
                end
            end
        end
    end
end
writeeventlist(file::AbstractString, args...) =
    open(fd -> writeeventlist(fd,args...), file, "w")


text:
"""
    exactShift!(hh,qq, iq, qrows, qcols, neq)

Compute the exact shiftrights and store them in qq.
"""
function exactShift!(hh::Array{Float64,2}, qq::Array{Float64,2}, iq::Int64, qRows::Int64, qCols::Int64, neq::Int64)

    # total number of shifts
    nexact = 0

    # functions to seperate hh
    left = 1:qCols
    right = (qCols + 1):(qCols + neq)

    # get right most columns of hh
    zerorows = copy(hh[:, right]')

    # compute absolute value
    zerorows = abs.(zerorows)

    # take the sum of the rows (aka 1st dimenison in julia)
    zerorows = Compat.sum(zerorows; dims=1)

    # anon function returns index of the rows who sum is 0
    zerorows = LinearIndices(zerorows)[findall(row->(row == 0), zerorows)]

    # continues until all zerorow indexes have been processed
    while length(zerorows) != 0 && (iq <= qRows)
        nz = size(zerorows, 1)

        # insert only the indexes found with zerorows into qq
        qq[(iq + 1):(iq + nz), :] = hh[zerorows, left]

	# update by shifting right by $neq columns
        hh[zerorows,:] = shiftRight!(hh[zerorows, :], neq)

        # increment the variables the amount of zerorows found
        iq = iq + nz
        nexact = nexact + nz

        # update zerorows as before but with our new hh matrix
        zerorows = hh[:, right]'
        zerorows = abs.(zerorows)
        zerorows = Compat.sum(zerorows; dims=1)
	zerorows = LinearIndices(zerorows)[findall(row->(row == 0), zerorows)]
    end # while

    return (hh, qq, iq, nexact)

end # exactShift


text:
module ParallelValueIteration

using Printf
using POMDPs
using POMDPModelTools
using LocalFunctionApproximation
using GridInterpolations
using Distributed
using SharedArrays

import POMDPs: Solver, solve, Policy, action, value

export
    ValueIterationPolicy,
    ParallelValueIterationSolver,
    solve,
    action,
    value

include("parallelSolver.jl")

end # module

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
module TestDay09

import AOC2021.Day09.part1, AOC2021.Day09.part2, AOC2021.Day09.string_block_to_matrix
using Test

test_input = """
2199943210
3987894921
9856789892
8767896789
9899965678"""

@testset "Day 09" begin
    @testset "part 1" begin
        @test part1(test_input) == 15
    end
    @testset "part 2" begin
        @test part2(test_input) == 1134
    end
end

end # module


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "AppVeyor" begin
    @testset "Plugin creation" begin
        p = AppVeyor()
        @test isempty(p.gitignore)
        @test p.src == joinpath(PkgTemplates.DEFAULTS_DIR, "appveyor.yml")
        @test p.dest == ".appveyor.yml"
        @test p.badges == [
            Badge(
                "Build Status",
                "https://ci.appveyor.com/api/projects/status/github/{{USER}}/{{PKGNAME}}.jl?svg=true",
                "https://ci.appveyor.com/project/{{USER}}/{{PKGNAME}}-jl",
            ),
        ]
        @test isempty(p.view)
        p = AppVeyor(; config_file=nothing)
        @test p.src === nothing
        p = AppVeyor(; config_file=test_file)
        @test p.src == test_file
        @test_throws ArgumentError AppVeyor(; config_file=fake_path)
    end

    @testset "Badge generation" begin
        p = AppVeyor()
        @test badges(p, me, test_pkg) == ["[![Build Status](https://ci.appveyor.com/api/projects/status/github/$me/$test_pkg.jl?svg=true)](https://ci.appveyor.com/project/$me/$test_pkg-jl)"]
    end

    @testset "File generation" begin
        # Without a coverage plugin in the template, there should be no post-test step.
        p = AppVeyor()
        @test gen_plugin(p, t, test_pkg) == [".appveyor.yml"]
        @test isfile(joinpath(pkg_dir, ".appveyor.yml"))
        appveyor = read(joinpath(pkg_dir, ".appveyor.yml"), String)
        @test !occursin("on_success", appveyor)
        @test !occursin("%JL_CODECOV_SCRIPT%", appveyor)
        rm(joinpath(pkg_dir, ".appveyor.yml"))

        # Generating the plugin with Codecov in the template should create a post-test step.
        t.plugins[Codecov] = Codecov()
        gen_plugin(p, t, test_pkg)
        delete!(t.plugins, Codecov)
        appveyor = read(joinpath(pkg_dir, ".appveyor.yml"), String)
        @test occursin("on_success", appveyor)
        @test occursin("%JL_CODECOV_SCRIPT%", appveyor)
        rm(joinpath(pkg_dir, ".appveyor.yml"))

        # TODO: Add Coveralls tests when AppVeyor.jl supports it.

        p = AppVeyor(; config_file=nothing)
        @test isempty(gen_plugin(p, t, test_pkg))
        @test !isfile(joinpath(pkg_dir, ".appveyor.yml"))
    end
end

rm(pkg_dir; recursive=true)


text:
"""
Parameters for the sampling solver.
"""
mutable struct SimpleSamplerParams{Tv,Ti}
    startingSize::Ti    # the initial size of the linked list storage structure
    blockSize::Ti       # the size of each consecutive block of memory assigned to a certain element

    # debug parameters, we can get rid of them later
    verboseSS::Bool
    verbose::Bool
    returnCN::Bool
    CNTol::Tv
    order::Symbol    # options :min, :approx, :tree
end

mutable struct AugTreeHybridParams
    n0::Int64      # the number of edges at which to go direct
    frac::Float64  # add frac*n edges to the tree
    subIts::Float64  # number of its of sub solver
    subTol::Float64  # tol to spec for sub solver
    verbose::Bool
    sampParams::SimpleSamplerParams
end


text:
function predict_mu(gp::GPE, X::AbstractMatrix, cK::AbstractMatrix)
    n = size(X, 2)
    mu = mean(gp.mean, X) + cK*gp.alpha        # Predictive mean
    return mu
end

""" Copy of a Gaussian process that shares X and cK but allows
    modification of the outcomes Y. This is useful for bootstrapping.
"""
function modifiable(gp::GPE)
    gp_copy = GPE(gp.x, copy(gp.y), gp.mean, gp.kernel,
                  gp.logNoise, gp.covstrat, gp.data, gp.cK)
    return gp_copy
end
function update_alpha!(gp::GPE)
    m = mean(gp.mean,gp.x)
    gp.alpha = gp.cK \ (gp.y - m)
end
function mll(gp::GPE)
    μ = mean(gp.mean,gp.x)
    return -dot((gp.y - μ),gp.alpha)/2.0 - logdet(gp.cK)/2.0 - gp.nobs*log(2π)/2.0 # Marginal log-likelihood
end

function addcov!(cK::AbstractMatrix{Float64}, kernel::Kernel, X::AbstractMatrix{Float64}, data::KernelData)
    dim, nobs = size(X)
    @assert size(cK, 1) == nobs
    @assert size(cK, 2) == nobs
    @inbounds for j in 1:nobs
        for i in 1:j
            cK[i,j] += cov_ij(kernel, X, X, data, i, j, dim)
            cK[j,i] = cK[i,j]
        end
    end
end
function add_diag!(mat::AbstractMatrix, d::Real)
    n = size(mat, 1)
    if n!=size(mat, 2)
        throw("Matrix is not square")
    end
    @inbounds for i in 1:n
        mat[i,i] += d
    end
    return mat
end
function update_chol!(pd::PDMats.PDMat)
    Σbuffer = mat(pd)
    Σbuffer, chol = make_posdef!(Σbuffer, cholfactors(pd))
    new_pd = wrap_cK(pd, Σbuffer, chol)
    return new_pd
end


text:
# file for containing constants required for using MeteomaticsAPI.
# currently limited to user credentials (which can be overridden in a function call)

username = "python-community"
password = "Umivipawe179"

max_stars_repo_name: meteomatics/julia-connector-api

max_issues_repo_name: meteomatics/julia-connector-api

max_forks_repo_name: meteomatics/julia-connector-api


text:
include("Line.jl")
include("Edge.jl")
include("Ellipse.jl")


text:
import J4Org: tokenized, untokenize
import J4Org: identifier
import J4Org: skip
import J4Org: raw_string
import J4Org: extract_comment
import J4Org: extract_variable
import J4Org: extract_function
import J4Org: extract_struct
import J4Org: extract_export
import J4Org: extract_abstract
import J4Org: extract_enum
import J4Org: extract_macro

@testset "extract_comment" begin
    filename="$(dirname(@__FILE__))/code_examples/basic.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_comment(t,1)
    @test r!=nothing
    @test raw_string(r) == "# This is a test file\n"
    r = extract_comment(t,skip(r))
    @test r!=nothing
    @test raw_string(r) == "#+MyTag1,MyTag2\n#     \n# My documentedItem\n"
end;

@testset "extract_comment_tab" begin
    filename="$(dirname(@__FILE__))/code_examples/basic_tab.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_comment(t,1)
    @test r!=nothing
    @test raw_string(r)  == "# This is a test file\n"

    r = extract_comment(t,3)
    @test r!=nothing
    @test raw_string(r) == "#+MyTag1,MyTag2\n#     \n# My documentedItem\n"
end;

@testset "function_1" begin
    filename="$(dirname(@__FILE__))/code_examples/function_1.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "   function    myFunction (a;\n                           x::Float64=sin(1))"
end;

@testset "function_2" begin
    filename="$(dirname(@__FILE__))/code_examples/function_2.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "function myFunction (a{T};\n                     x::Float64=sin(1))::Bool"
end;

@testset "function_3" begin
    filename="$(dirname(@__FILE__))/code_examples/function_3.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "function myFunction (a{T};\n                     x::Float64=sin(1))::Bool where {T}"
end;

@testset "function_1_short" begin
    filename="$(dirname(@__FILE__))/code_examples/function_1_short.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "       myFunction (a;\n                   x::Float64=sin(1))"
end;

@testset "function_2_short" begin
    filename="$(dirname(@__FILE__))/code_examples/function_2_short.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "myFunction (a{T};\n            x::Float64=sin(1))::Bool"
end;

@testset "function_3_short" begin
    filename="$(dirname(@__FILE__))/code_examples/function_3_short.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_function(t,1)
    @test r!=nothing
    @test raw_string(r) == "myFunction (a{T};\n            x::Float64=sin(1))::Bool where {T}"
end;

@testset "struct_1" begin
    filename="$(dirname(@__FILE__))/code_examples/struct_1.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_struct(t,1)
    @test r!=nothing
    @test raw_string(r) == "mutable struct Bar"
end;

@testset "struct_2" begin
    filename="$(dirname(@__FILE__))/code_examples/struct_2.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_struct(t,1)
    @test r!=nothing
    @test raw_string(r) == "struct Bar"
end;

@testset "struct_3" begin
    filename="$(dirname(@__FILE__))/code_examples/struct_3.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_struct(t,1)
    @test r!=nothing
    @test raw_string(r) == "struct Point{T<:Real} <: Pointy{T}"
end;

@testset "export" begin
    filename="$(dirname(@__FILE__))/code_examples/export.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_export(t,1)
    @test r!=nothing
    @test raw_string(r) == "export A,B,   C"
end;

@testset "abstract" begin
    filename="$(dirname(@__FILE__))/code_examples/abstract.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_abstract(t,1)
    @test r!=nothing
    @test raw_string(r) == "abstract type UDWT_Filter_Biorthogonal{T<:Number} "
end;

@testset "enum" begin
    filename="$(dirname(@__FILE__))/code_examples/enum.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_enum(t,1)
    @test r!=nothing
    @test raw_string(r) == "@enum Alpha A  B=1 C"
    @test identifier(r) == "Alpha"
end;

@testset "variable" begin
    filename="$(dirname(@__FILE__))/code_examples/variable.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_variable(t,1)
    @test r!=nothing
    @test raw_string(r) == "A"
    @test identifier(r) == "A"
end;

@testset "variable_const" begin
    filename="$(dirname(@__FILE__))/code_examples/variable_const.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_variable(t,1)
    @test r!=nothing
    @test raw_string(r) == "const A"
    @test identifier(r) == "A"
end;

@testset "variable_const_global" begin
    filename="$(dirname(@__FILE__))/code_examples/variable_const_global.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_variable(t,1)
    @test r!=nothing
    @test raw_string(r) == "const global A"
    @test identifier(r) == "A"
end;

@testset "variable_global_const" begin
    filename="$(dirname(@__FILE__))/code_examples/variable_global_const.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_variable(t,1)
    @test r===nothing
end;

@testset "macro" begin
    filename="$(dirname(@__FILE__))/code_examples/macro.jl"
    code=readcode(filename)
    t=tokenized(code)

    r = extract_macro(t,1)
    @test r!=nothing
    @test raw_string(r) == "macro swap(x,y)"
    @test identifier(r) == "swap"
end;


text:
abstract type AbstractKalmanFilter <: AbstractFilter end

@with_kw struct KalmanFilter{AT,BT,CT,DT,R1T,R2T,R2DT,D0T,XT,RT} <: AbstractKalmanFilter
    A::AT
    B::BT
    C::CT
    D::DT
    R1::R1T
    R2::R2T
    R2d::R2DT
    d0::D0T
    x::XT
    R::RT
    t::Ref{Int} = Ref(1)
end

"""
    KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))

The matrices `A,B,C,D` define the dynamics
```
x' = Ax + Bu + w
y  = Cx + Du + e
```
where `w ~ N(0, R1)`, `e ~ N(0, R2)` and `x(0) ~ d0`
"""
function KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(Matrix(R1)))
    try
        cR1 = cond(R1)
        cR2 = cond(R2)
        (cond(cR1) > 1e8 || cond(cR2) > 1e8) && @warn("Covariance matrices are poorly conditioned")
    catch
        nothing
    end

    KalmanFilter(A,B,C,D,R1,R2,MvNormal(Matrix(R2)), d0, Vector(d0.μ), Matrix(d0.Σ), Ref(1))
end

function Base.propertynames(kf::KF, private::Bool=false) where KF <: AbstractKalmanFilter
    return fieldnames(KF)
end

sample_state(kf::AbstractKalmanFilter) = rand(kf.d0)
sample_state(kf::AbstractKalmanFilter, x, u, t) = kf.A*x .+ kf.B*u .+ rand(MvNormal(kf.R1))
sample_measurement(kf::AbstractKalmanFilter, x, u, t) = kf.C*x .+ kf.D*u .+ rand(MvNormal(kf.R2))
particletype(kf::AbstractKalmanFilter) = typeof(kf.x)
covtype(kf::AbstractKalmanFilter)      = typeof(kf.R)
state(kf::AbstractKalmanFilter)        = kf.x
covariance(kf::AbstractKalmanFilter)   = kf.R

function reset!(kf::AbstractKalmanFilter)
    kf.x .= Vector(kf.d0.μ)
    kf.R .= copy(Matrix(kf.d0.Σ))
    kf.t[] = 1
end

# UKF ==========================================================================

@with_kw struct UnscentedKalmanFilter{DT,MT,R1T,R2T,R2DT,D0T,VT,XT,RT} <: AbstractKalmanFilter
    dynamics::DT
    measurement::MT
    R1::R1T
    R2::R2T
    R2d::R2DT
    d0::D0T
    xs::Vector{VT}
    x::XT
    R::RT
    t::Ref{Int} = Ref(1)
end

"""
    UnscentedKalmanFilter(dynamics,measurement,R1,R2,d0=MvNormal(Matrix(R1)))
"""
function UnscentedKalmanFilter(dynamics,measurement,R1,R2,d0=MvNormal(Matrix(R1)))
    try
        cR1 = cond(R1)
        cR2 = cond(R2)
        (cond(cR1) > 1e8 || cond(cR2) > 1e8) && @warn("Covariance matrices are poorly conditioned")
    catch
        nothing
    end
    n = size(R1,1)
    p = size(R2,1)
    R1s = SMatrix{n,n}(R1)
    R2s = SMatrix{p,p}(R2)
    xs = sigmapoints(mean(d0), cov(d0))
    UnscentedKalmanFilter(dynamics,measurement,R1s,R2s,MvNormal(Matrix(R2s)), d0, xs, Vector(d0.μ), Matrix(d0.Σ), Ref(1))
end

sample_state(kf::UnscentedKalmanFilter) = rand(kf.d0)
sample_state(kf::UnscentedKalmanFilter, x, u, t) = kf.dynamics(x,u,t) .+ rand(MvNormal(Matrix(kf.R1)))
sample_measurement(kf::UnscentedKalmanFilter, x, u, t) = kf.measurement(x, u, t) .+ rand(MvNormal(Matrix(kf.R2)))

# function transform_moments!(S,X,m,L)
#     X .-= mean(X) # Normalize the sample
#     for i in eachindex(X)
#         S[i] = m .+ L*X[i]
#     end
#     S
# end

function sigmapoints(m, Σ)
    T = promote_type(eltype(m), eltype(Σ))
    n = max(length(m), size(Σ,1))
    xs = [@SVector zeros(T, n) for _ in 1:(2n+1)]
    sigmapoints!(xs,m,Σ)
end

function sigmapoints!(xs, m, Σ::AbstractMatrix)
    n = length(xs[1])
    X = sqrt(Symmetric(n*Σ))
    for i in 1:n
        xs[i] = X[:,i]
        xs[i+n] = -xs[i] .+ m
        xs[i] = xs[i] .+ m
    end
    xs[end] = m
    xs
end

function predict!(ukf::UnscentedKalmanFilter, u, t::Integer = index(ukf))
    @unpack dynamics,measurement,x,xs,R,R1 = ukf
    ns = length(xs)
    sigmapoints!(xs,x,R) # TODO: these are calculated in the update step
    r = copy(R)
    for i in eachindex(xs)
        xs[i] = dynamics(xs[i], u, t)
    end
    x .= mean(xs)
    # for i in eachindex(xs)
    #     d = xs[i]-x
    #     r .+= Symmetric(d*d')
    # end
    # R .= Symmetric(r)./n + R1
    R .= symmetrize(cov(xs)) + R1
    ukf.t[] += 1
end

correct!(ukf::UnscentedKalmanFilter, y, t::Integer = index(ukf)) = correct!(ukf::UnscentedKalmanFilter, y, 0, t)
function correct!(ukf::UnscentedKalmanFilter, u, y, t::Integer = index(ukf))
    @unpack measurement,x,xs,R,R1,R2,R2d = ukf
    n = size(R1,1)
    p = size(R2,1)
    ns = length(xs)
    sigmapoints!(xs,x,R) # Update sigmapoints here since untransformed points required
    C = @SMatrix zeros(n,p)
    ys = map(xs) do x
        measurement(x, u, t)
    end
    ym = mean(ys)
    @inbounds for i in eachindex(ys) # Cross cov between x and y
        d   = ys[i]-ym
        ca  = (xs[i]-x)*d'
        C  += ca
    end
    e   = y .- ym
    S   = symmetrize(cov(ys)) + R2 # cov of y
    Sᵪ = cholesky(S)
    K   = (C./ns)/Sᵪ # ns normalization to make it a covariance matrix
    x .+= K*e
    # mul!(x, K, e, 1, 1) # K and e will be SVectors if ukf correctly initialized
    RmKSKT!(R, K, S)
    ll = logpdf(MvNormal(PDMat(S,Sᵪ)), e) #- 1/2*logdet(S) # logdet is included in logpdf
    ll, e
end

@inline function RmKSKT!(R, K, S)
    R .-= K*S*K'
    symmetrize(R)
    nothing
end

# SigmaFilter ==========================================================================

struct SigmaFilter{DT,MT,MLT,D0T,DFT,VT} <: AbstractParticleFilter
    dynamics::DT
    measurement::MT
    measurement_likelihood::MLT
    w::Vector{Float64}
    we::Vector{Float64}
    initial_density::D0T
    dynamics_density::DFT
    x::Vector{VT}
    xprev::Vector{VT}
    xm::Vector{Float64}
    R::Matrix{Float64}
    t::Ref{Int}
end

"""
    SigmaFilter(N,dynamics,measurement,measurement_likelihood,df,d0)
"""
function SigmaFilter(N,dynamics,measurement,measurement_likelihood,df,d0)
    @show n = length(d0)
    xs = rand(d0,N)
    xs = reinterpret(SVector{n,Float64}, xs)[:]
    SigmaFilter(dynamics,measurement,measurement_likelihood,
        zeros(length(xs)), # w
        zeros(length(xs)), # we
        d0, df,
        xs,                # x
        copy(xs),          # xprev
        Vector(d0.μ),      # xm
        Matrix(d0.Σ),      # R
        Ref(1))
end

sample_state(sf::SigmaFilter) = rand(sf.initial_density)
sample_state(sf::SigmaFilter, x, u, t) = sf.dynamics(x,u,t,true)
sample_measurement(sf::SigmaFilter, x, u, t) = sf.measurement(x,u,t,true)
num_particles(sf::SigmaFilter) = length(sf.x)
particles(sf::SigmaFilter) = sf.x
weights(sf::SigmaFilter) = sf.w
expweights(sf::SigmaFilter) = sf.we
state(sf::SigmaFilter) = sf
rng(sf::SigmaFilter) = Random.GLOBAL_RNG

function predict!(sf::SigmaFilter, u, t::Integer = index(sf))
    @unpack dynamics,measurement,x,xprev,xm,R,w,we = sf
    N = length(x)
    n = length(x[1])
    xm .= sum(x->x[1]*x[2], zip(x,we)) # Must update based on correction
    R .= 0
    for i in eachindex(x)
        δ = x[i]-xm
        R .+= (δ*δ')*we[i]
    end
    R ./= (1-sum(abs2,we))
    R .= (R + R')/2 + 1e-16I
    d = MvNormal(xm,R)
    noisevec = zeros(length(d))
    for i in eachindex(x)
        xi = SVector{n,Float64}(rand!(d, noisevec))
        x[i] = dynamics(xi, u, t, true)
    end
    w  .= -log(N)
    we .= 1/N
    xm .= mean(x)
    R  .= cov(x)
    sf.t[] += 1
end

Base.@propagate_inbounds propagate_particles!(sf::SigmaFilter, u, j::Vector{Int}, t::Int, noise=true) = propagate_particles!(sf, u, t, noise)

# Base.@propagate_inbounds function measurement_equation!(sf::SigmaFilter, y, t, w = weights(sf))
#     g = measurement_likelihood(sf)
#     any(ismissing.(y)) && return w
#     x = particles(sf)
#     @inbounds for i = 1:num_particles(sf)
#         w[i] += g(x[i],y,t)
#     end
#     w
# end


text:
cd(joinpath(@__DIR__,"..","test")) do
    # clean out the artifacts directory
    artifacts = joinpath(@__DIR__,"..","test","artifacts")
    isdir(artifacts) && rm(artifacts,recursive=true)
    mkdir(artifacts)

    include(joinpath(@__DIR__,"..","test","noise_analysis.jl"))
    include(joinpath(@__DIR__,"..","test","ljh.jl"))
    include(joinpath(@__DIR__,"..","test","ljh3.jl"))
end


text:
# As of 8/23/2017, this program does not fully run.
# Since the "fixed" option is not used in the current nearly-neutral paper, this option has not been kept up to date.
# Front end for src/nn_poplist.jl
# For the case where selection coefficients are drawn from a gamma distribution or a double-sided gamma distributions
#include("../src/NearlyNeutral.jl")
#include("../src/InfAlleles.jl")
include("../src/InfSites.jl")
if length(ARGS) == 0
  simname = "../experiments/examples/fi_example2"
else
  simname = ARGS[1]
end
include("$(simname).jl")
println("simname: ",simname)
stream = open("$(simname).csv","w")
println("stream: ",stream)

current_dir = pwd()
date_string = "../data/"*Dates.format(now(),"mm_dd_yy")*"/"
#try mkdir(date_string) catch end   # create today's directory with no error if it already exists
println("date string: ",date_string)
if !isdefined(:mu_list_flag)
  mu_list_flag=false
end

@doc """ type trial_result
  A universal trial result for all run_trials functions
  Constructors for specific cases are defined below.
"""
type trial_result
  nn_simtype::Int64
  n::Int64    # sample size, must be <= N
  N::Int64    # popsize
  N_mu::Float64  # N*mu, population mutation rate
  ngens::Int64
  burn_in::Float64
  fix_minimum::Float64
  dfe::Function
  dfe_s::Float64
  dfe_modulus::Int64
  expected_richness::Float64  # sum_{i=0}^{n-1} theta/(theta+i) where theta = 2*N_mu.
  average_richness::Float64   # Average number of traits in populations
  average_richness_stddev::Float64   # Average number of traits in populations
  expected_w_homoz::Float64
  w_homoz::Float64
  w_homoz_stddev::Float64
  IQV::Float64
  IQV_stddev::Float64
  number_active::Int64
  number_extinctions::Int64
  number_fixations::Int64
  extinction_time::Float64
  fixation_time::Float64
end

function add_expected_richness( tr::trial_result )
  sum = 0.0
  theta = 2*tr.N_mu
  for i = 0:(tr.N-1)
    sum += 1.0/(theta+i)
  end
  tr.expected_richness = theta*sum
end

# Constructor that sets the parameters
function trial_result( nn_symtype::Int64, n::Int64, N::Int64, N_mu::Float64, ngens::Int64,
    fix_minimum::Float64=0.5, burn_in::Float64=2.0, dfe::Function=dfe_neutral,
    dfe_s::Float64=0.0, dfe_modulus::Int64=10 )
  #tr = trial_result( nn_simtype, n, N, N_mu, ngens, burn_in, fix_minimum, dfe, dfe_s, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0 )
  tr = trial_result( nn_simtype, n, N, N_mu, ngens, burn_in, fix_minimum, dfe, dfe_s, dfe_modulus, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0.0, 0.0 )
  add_expected_richness( tr )
  tr
end

function print_trial_result( tr::trial_result )
  if tr.nn_simtype == 1
    println("\ninfinite alleleles model")
  elseif tr.nn_simtype == 2
    println("\ninfinite sites model")
  else
    println("nn_simtype: ", tr.nn_simtype)
  end
  println("n: ", tr.n)
  println("N: ", tr.N)
  println("N_mu: ", tr.N_mu)
  println("mu: ", tr.N_mu/tr.N)
  println("ngens: ", tr.ngens)
  println("burn_in: ", tr.burn_in)
  println("fix_minimum: ", tr.fix_minimum)
  println("dfe: ", tr.dfe)
  println("dfe_s: ", tr.dfe_s)
  println("dfe_modulus: ", tr.dfe_modulus)
  println("expected richness: ",tr.expected_richness)
  println("average richness: ",tr.average_richness)
  println("average richness_stddev: ",tr.average_richness_stddev)
  println("expected_w_homoz: ",tr.expected_w_homoz)
  println("w_homoz: ",tr.w_homoz)
  println("w_homoz_stddev: ",tr.w_homoz_stddev)
  println("IQV: ",tr.IQV)
  println("IQV_stddev: ",tr.IQV_stddev)
  println("number_active: ", tr.number_active)
  println("number_extinctions: ", tr.number_extinctions)
  println("number_fixations: ", tr.number_fixations)
  println("extinction_time: ", tr.extinction_time)
  println("fixation_time: ", tr.fixation_time)
end

# Convert an innovation_oollection to a trial_result
# The innovation_connection and the trial_result must agree in fields n, N,
function convert_to_trial_result( ic::innovation_collection, tr::trial_result )
  tr.number_active = length(ic.active)
  tr.number_extinctions = length(ic.extinct)
  tr.number_fixations = length(ic.fixed)
  tr.fix_minimum = ic.fix_minimum
  tr.extinction_time = average_time_to_extinction( ic )
  tr.fixation_time = average_time_to_fixation( ic )
  #tr.average_fitness_fixed = average_fitness_fixed( ic::innovation_collection )
  #tr.average_fitness_extinct = average_fitness_extinct( ic::innovation_collection )
  #tr.average_fitness_all = average_fitness_all( ic::innovation_collection )
end

function run_trials(popsize_multiplier_list::Vector{Int64}=[1]; mu_list_flag::Bool=false)
  println("stream: ",stream)
  trial = 1
  N = N_list[1]
  n = Int(floor(N*(1//popsize_multiplier_list[1])))
  if !mu_list_flag
    tr = trial_result( nn_simtype, n, N, N_mu_list[1], ngens, fix_minimum, burn_in, InfSites.dfe_mod, dfe_s_list[1], dfe_modulus )
    writeheader(stream, popsize_multiplier_list, N_list, N_mu_list, tr )
    for dfe_s in dfe_s_list
      dfe = x->dfe_mod( x, modulus=dfe_modulus, fit_inc=1.0+dfe_s )
      for N_mu in N_mu_list
        for N in N_list
          for psize_m in popsize_multiplier_list
            n = Int(floor(N*(1//psize_m)))
            tr = trial_result( nn_simtype, n, N, N_mu, ngens, fix_minimum, burn_in, dfe, dfe_s, dfe_modulus )
            run_trial( tr )
            writerow(stream, trial, tr )
            trial += 1
          end
        end
      end
    end
  else  # if mu_list_flag
    tr = trial_result( nn_simtype, n, N, N*mu_list[1], ngens, fix_minimum, burn_in, dfe_mod, dfe_s_list[1], dfe_modulus )
    writeheader(stream, popsize_multiplier_list, N_list, mu_list, tr, mu_list_flag=mu_list_flag )
    for dfe_s in dfe_s_list
      dfe = x->dfe_mod( x, modulus=dfe_mod, fit_inc=dfe_s )
      for mu in mu_list
        for N in N_list
          for psize_m in popsize_multiplier_list
            n = Int(floor(N*(1//psize_m)))
            N_mu = N*mu
            tr = trial_result( nn_simtype, n, N, N_mu, ngens, fix_minimum, burn_in, dfe, dfe_s, dfe_modulus )
            run_trial( tr )
            writerow(stream, trial, tr, mu_list_flag=mu_list_flag )
            trial += 1
          end
        end
      end
    end
  end  # if !mu_list_flag
end

function run_trials_mu(popsize_multiplier_list::Vector{Int64}=[1])
  println("stream: ",stream)
  trial = 1
  N = N_list[1]
  n = Int(floor(N*(1//popsize_multiplier_list[1])))
  tr = trial_result( nn_simtype, n, N, N*mu_list[1], ngens, fix_minimum, burn_in, dfe_mod, dfe_s_list[1], dfe_modulus )
  writeheader(stream, popsize_multiplier_list, N_list, mu_list, tr )
  for dfe_s in dfe_s_list
    dfe = x->dfe_mod( x, modulus=dfe_mod, fit_inc=dfe_s )
    for mu in mu_list
      for N in N_list
        for psize_m in popsize_multiplier_list
          n = Int(floor(N*(1//psize_m)))
          N_mu = N*mu
          tr = trial_result( nn_simtype, n, N, N_mu, ngens, fix_minimum, burn_in, dfe, dfe_s, dfe_modulus )
          run_trial( tr )
          writerow(stream, trial, tr )
          trial += 1
        end
      end
    end
  end
end

function run_trial( tr::trial_result )
  if tr.nn_simtype == 1
    println("dfe(tr.dfe_modulus): ",tr.dfe(tr.dfe_modulus))
    ic = innovation_collection( tr.N, tr.fix_minimum )
    poplist = nn_poplist(tr.N,tr.N_mu,tr.ngens,tr.dfe,combine=false,ic=ic)
    convert_to_trial_result( ic, tr )
    if tr.n < tr.N
      poplist = map(x->sample_population(x,tr.n),poplist)
    end
    add_stats_to_trial_result!( tr, poplist )
    print_trial_result( tr )
    return tr
  else
    println("nn_simtype ",nn_simtype," not implemented")
  end
end

function add_stats_to_trial_result!( tr::trial_result, poplist::Vector{Population} )
  #println("poplist: ",poplist)
  #println("poplist fitnesses: ",map(y->map(tr.dfe,y)),poplist)
  mu = tr.N_mu/tr.N
  filtered_poplist = map(y->filter(x->(x%tr.dfe_modulus==0),y),poplist)
  filtered_N = map(length,filtered_poplist)  # vector of number of individuals whose modulus is dfe_modulus
  #println("filtered_poplist: ",filtered_poplist)
  filtered_pcounts = map(pop_counts64,filtered_poplist)
  richness = map(length,filtered_pcounts)  # vector of richnesses over generations
  w_homoz = map(watterson_homozygosity,filtered_pcounts)
  IQV_lst = map(IQV,filtered_pcounts)
  expected_filtered_N = tr.N/tr.dfe_modulus   # Expected number of indivs whose modulus is dfe_modulus
  println("length(w_homoz): ",length(w_homoz),"  length(IQV_lst): ",length(IQV_lst),"  mean_filtered_N: ",mean(filtered_N))
  println("w_homoz[1]: ",w_homoz[1],"  sum(w_homoz): ",sum(w_homoz))
  tr.average_richness = mean(richness)
  tr.average_richness_stddev = std(richness)
  #tr.expected_w_homoz = 1.0/(1.0+2.0*tr.N_mu)
  tr.expected_w_homoz = 1.0/(1.0+2.0*expected_filtered_N*mu)
  tr.w_homoz = mean(w_homoz)
  println("tr.w_homoz: ",tr.w_homoz)
  tr.w_homoz_stddev = std(w_homoz)
  tr.IQV = mean(IQV_lst)
  tr.IQV_stddev = std(IQV_lst)
end

@doc """ function writeheader()
    burn_in::Float64, slat_reps::Int64=100000 )
Write header line for output CSV file that corresponds to stream.
See the comments for run_simulation for a description of the parameters.
"""
function writeheader(stream::IO, popsize_multiplier_list::Vector{Int64}, N_list::Vector{Int64},
    mu_list::Vector{Float64}, tr::trial_result; mu_list_flag::Bool=false )
  dfe_params = [:dfe_adv_prob, :dfe_adv_alpha, :dfe_adv_theta, :dfe_disadv_prob, :dfe_disadv_alpha, :dfe_disadv_theta]
  param_strings = [
    "# $(string(Dates.today()))",
    "# $((nn_simtype==1)?"infinite alleles model":"infinite_sites_model")",
    "# popsize_multiplier_list=$(popsize_multiplier_list)",
    "# N_list=$(N_list)",
    mu_list_flag ? "# mu_list=$(mu_list)": # N_mu_list=$(N_mu_list)",
    "# dfe_modulus=$(tr.dfe_modulus)",
    "# ngens=$(tr.ngens)",
    "# burn_in=$(tr.burn_in)",
    "# fix_minimum=$(tr.fix_minimum)",
    #"# dfe=$(tr.dfe)",
    #"# dfe_s=$(tr.dfe_s)"
    ]
  for dp in dfe_params
    if isdefined(dp)
      Base.push!(param_strings,"# $(string(dp))=$(eval(dp))")
    end
  end
  write(stream, join(param_strings, "\n"), "\n")
  first_heads = [
    "type",
    #  "n",
    "N",
    (mu_list_flag ? "mu":"N_mu"),
    "s"
  ]
  last_heads =
  [ "expected_richness",
    "average_richness",
    "average_richness_stddev",
    "exp_neut_w_heteroz",
    "w_heteroz",
    "w_heteroz_stddev",
    "IQV",
    "IQV_stddev",
    "num_extinct",
    "num_fixed",
    "fraction_fixed",
    "ave_extinct_time",
    "ave_fixed_time"
  ]
  line = join(vcat( first_heads, last_heads), ",")
  write(stream, line, "\n")
end

function writerow(stream::IO, trial::Int64, tr::trial_result; mu_list_flag::Bool=false   )
  first = Any[
    (isdefined(:type_str) ? type_str :trial),
    #tr.n,           # sample size
    tr.N,           # popsize
    mu_list_flag ? tr.N_mu/tr.N:tr.N_mu,
    tr.dfe_s,
  ]
  if tr.nn_simtype == 0
    mid = Any[]
  elseif tr.nn_simtype == 1
    mid = Any[
      tr.expected_richness,
      tr.average_richness,
      tr.average_richness_stddev,
      1.0-tr.expected_w_homoz,
      1.0-tr.w_homoz,
      tr.w_homoz_stddev,
      tr.IQV,
      tr.IQV_stddev,
      tr.number_extinctions,
      tr.number_fixations,
      Float64(tr.number_fixations)/(tr.number_extinctions+tr.number_fixations),
      tr.extinction_time,
      tr.fixation_time
    ]
  end
  line = join( vcat( first, mid ), "," )
  write(stream, line, "\n")
end

if isdefined(:popsize_multiplier_list)
  run_trials(popsize_multiplier_list, mu_list_flag=mu_list_flag )
else
  run_trials( mu_list_flag=mu_list_flag)
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using BenchmarkTools: leaves
using PkgBenchmark: benchmarkgroup
using TransducersBenchmarksReports: loadresult
using DataFrames
using Statistics
using VegaLite

fanout(; kwargs...) =
    x -> NamedTuple{keys(kwargs)}(Ref(x) .|> Tuple(values(kwargs)))

macro fanout(names...)
    kwargs = [:($n=$n) for n in names]
    esc(:($fanout(; $(kwargs...))))
end

function trialasrow((key, trial),)
    # stats = @fanout(mean, median, var, minimum, maximum)(trial.times)
    stats = @fanout(mean, median, minimum)(trial.times)
    if length(key) == 3
        impl, simd, n = key
    else
        impl, n = key
        simd = "none"
    end
    return (impl=impl, simd=simd, n=n, stats...)
end

group = benchmarkgroup(loadresult())["gemm"]
data = trialasrow.(leaves(group["mul"]))

df = melt(DataFrame(data), [:impl, :simd, :n])

dfnorm = join(
    df[df[:impl] .!= "linalg", :],
    deletecols!(df[df[:impl] .== "linalg", :], [:impl, :simd]),
    on = [:n, :variable],
    kind = :left,
    makeunique = true,
)
dfnorm[:normalized] = dfnorm[:value] ./ dfnorm[:value_1]
dfnorm[:fullname] = string.(dfnorm[:impl], ", ", dfnorm[:simd])

plt = dfnorm |> @vlplot(
    :text,
    x = {
        field = :n,
        typ = :quantitative,
        scale = {typ = :log},
    },
    y = {
        field = :normalized,
        typ = :quantitative,
        scale = {typ = :log},
    },
    column = :variable,
    # color = :fullname,
    color = :simd,
    text = :impl,
    resolve = {
        scale = {y = :independent},
    }
)


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function basic_checkbox(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.Checkbox(store, s->get_state(s).checkboxes["basic_checkbox"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end
end


text:
"""
    format(path::String=pwd())

Run formatter at `path`, the `.JuliaFormatter.toml`
should be used to config the formatter.

# Args

- `path`: path to the folder/file to format.
"""
@cast function format(path::String=pwd())
    JuliaFormatter.format(path)
    return
end


text:
using Plots
plotly()

function Φ(z, λ, w0)
    zR = pi*w0^2/λ
    return atan(z/zR)
end

function w(z, λ, w0)
    zR = pi*w0^2/λ
    return w0*sqrt(1+(z/zR)^2)
end

function R(z, λ, w0)
    zR = pi*w0^2/λ
    return z+zR^2/z
end

function q(z, λ, w0)
    zR = pi*w0^2/λ
    return z + zR*1im
end

function ψ(z, p, m, λ, w0)
    zR = pi*w0^2/λ
    N = abs(m) + 2p
    return (N+1)*atan(z/zR)
end

δ(k::Integer,j::Integer) = k == j ? 1 : 0

function gLaguerre(x, p, m)
    if p == 0
        return 1
    elseif p == 1
        return 1+m-x
    else
        k = p-1
        return ((2k+1+m-x)*gLaguerre(x, k, m)-(k+m)*gLaguerre(x, k-1, m))/(k+1)
    end
end

function H(x,n)
    if n==0
        return 1
    elseif n==1
        return 2x
    else
        return 2x*H(x, n-1)-2*(n-1)*H(x, n-2)
    end
end

function r_and_θ(x, y)
    r = sqrt.(x.^2 .+ y'.^2)
    θ = atan.(y', x)
    return r, θ
end

#Laguerre-Gaussian modes
function LG(r, θ, z, p, m; λ=800e-9, w0=2e-2)
    k = 2pi/λ
    _1 = @. sqrt(2*factorial(p)/(pi*factorial(abs(m)+p)))
    _2 = @. w0/w(z, λ, w0)*(sqrt(2)*r/w(z, λ, w0))^abs(m)
    _3 = @. exp(-r^2/w(z, λ, w0)^2)
    _4 = @. gLaguerre(2r^2/w(z, λ, w0)^2,p,abs(m))
    _5 = @. exp(-1im*k*r^2/(2*R(z, λ, w0))-1im*θ*m)
    _6 = @. exp(1im*ψ(z, p, m, λ, w0))
    return  @. _1*_2*_3*_4*_5*_6
    # return _3*_4
end

#Hermite-Gaussian modes
function HG(x, y, z, l, m; λ=800e-9, w0=2e-2)
    k = 2pi/λ
    r, θ = r_and_θ(x, y)
    # outer product
    _1 =  w0/w(z, λ, w0).*H.(sqrt(2)*x/w(z, λ, w0), l).*H.(sqrt(2)*y/w(z, λ, w0), m)'
    _2 = @. exp(-(r^2)/w(z, λ, w0)^2)
    _3 = @. exp(-1im*k*(r^2)/(2*R(z, λ, w0)))
    _4 = @. exp(1im*ψ(z, l, m, λ, w0))
    return @. _1*_2*_3*_4
end

# x = LinRange{Float64}(-2,10, 100)
x = Vector(LinRange{Float64}(-5e-2,5e-2, 500))
y = Vector(LinRange{Float64}(-5e-2,5e-2, 500))
r, θ = r_and_θ(x, y)

# heatmap(x,y,real.(LG(r,θ,0,1,3)).^2)
# plot!(x,gLaguerre.(x, 2, 2))

# x = LinRange{Float64}(-2,3, 100)
# plot!(x, H.(x,4))
heatmap(x,y, real.(HG(x, y, 0, 2, 2)).^2)

# heatmap(r)

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
using ForwardDiff

import NLPModels.obj, NLPModels.grad, NLPModels.grad!, NLPModels.cons, NLPModels.cons!
import NLPModels.jac, NLPModels.jprod, NLPModels.jprod!, NLPModels.jtprod
import NLPModels.jtprod!, NLPModels.hess, NLPModels.hprod, NLPModels.hprod!
import NLPModels.increment!

export HOADNLPModel,
    obj,
    grad,
    grad!,
    cons,
    cons!,
    jac,
    jprod,
    jprod!,
    jtprod,
    jtprod!,
    hess,
    hprod,
    hprod!,
    ∇f³xuv,
    ∇f³xuv2

"""HOADNLPModel is an AbstractNLPModel using ForwardDiff to compute the
derivatives. HOHOADNLPModel uses ForwardDiff to compute 3rd order information.
This interface is very similar to HOADNLPModel with few additional functions which
are documented here.
For further documentation, refer to the HOADNLPModel.
"""
mutable struct HOADNLPModel <: AbstractNLPModel
    meta::NLPModelMeta

    counters::Counters

    # Functions
    f::Function
    c::Function

    # AD Packages
    AD::Any
end

function HOADNLPModel(
    f::Function,
    x0::AbstractVector;
    y0::AbstractVector = eltype(x0)[],
    lvar::AbstractVector = eltype(x0)[],
    uvar::AbstractVector = eltype(x0)[],
    lcon::AbstractVector = eltype(x0)[],
    ucon::AbstractVector = eltype(x0)[],
    c::Function = (args...) -> throw(NotImplementedError("cons")),
    name::String = "Generic",
    lin::AbstractVector{Int} = Int[],
    AD::Any = ForwardDiff,
)

    nvar = length(x0)
    length(lvar) == 0 && (lvar = -Inf * ones(nvar))
    length(uvar) == 0 && (uvar = Inf * ones(nvar))
    ncon = maximum([length(lcon); length(ucon); length(y0)])

    A = AD.hessian(f, x0)
    for i = 1:ncon
        A += AD.hessian(x -> c(x)[i], x0) * (-1)^i
    end
    nnzh = nvar * (nvar + 1) / 2
    nnzj = 0

    if ncon > 0
        length(lcon) == 0 && (lcon = -Inf * ones(ncon))
        length(ucon) == 0 && (ucon = Inf * ones(ncon))
        length(y0) == 0 && (y0 = zeros(ncon))
        nnzj = nvar * ncon
    end
    nln = setdiff(1:ncon, lin)

    meta = NLPModelMeta(
        nvar,
        x0 = x0,
        lvar = lvar,
        uvar = uvar,
        ncon = ncon,
        y0 = y0,
        lcon = lcon,
        ucon = ucon,
        nnzj = nnzj,
        nnzh = nnzh,
        lin = lin,
        nln = nln,
        minimize = true,
        islp = false,
        name = name,
    )

    return HOADNLPModel(meta, Counters(), f, c, AD)
end

function obj(nlp::HOADNLPModel, x::AbstractVector)
    # increment!(nlp, :neval_obj)
    return nlp.f(x)
end

function grad(nlp::HOADNLPModel, x::AbstractVector)
    # increment!(nlp, :neval_grad)
    return ForwardDiff.gradient(nlp.f, x)
end

function grad!(nlp::HOADNLPModel, x::AbstractVector, g::AbstractVector)
    # increment!(nlp, :neval_grad)
    ForwardDiff.gradient!(view(g, 1:length(x)), nlp.f, x)
    return g
end

function cons(nlp::HOADNLPModel, x::AbstractVector)
    # increment!(nlp, :neval_cons)
    return nlp.c(x)
end

function cons!(nlp::HOADNLPModel, x::AbstractVector, c::AbstractVector)
    increment!(nlp, :neval_cons)
    c[1:nlp.meta.ncon] = nlp.c(x)
    return c
end

function jac(nlp::HOADNLPModel, x::AbstractVector)
    # increment!(nlp, :neval_jac)
    return ForwardDiff.jacobian(nlp.c, x)
end

function jac_coord(nlp::HOADNLPModel, x::AbstractVector)
    Jx = jac(nlp, x)
    m, n = nlp.meta.ncon, nlp.meta.nvar
    if VERSION < v"1.0"
        I = [(i, j) for i = 1:m, j = 1:n]
    else
        I = ((i, j) for i = 1:m, j = 1:n)
    end
    return (getindex.(I, 1)[:], getindex.(I, 2)[:], Jx[:])
end

function jprod(nlp::HOADNLPModel, x::AbstractVector, v::AbstractVector)
    # increment!(nlp, :neval_jprod)
    return ForwardDiff.jacobian(nlp.c, x) * v
end

function jprod!(nlp::HOADNLPModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
    # increment!(nlp, :neval_jprod)
    Jv[1:nlp.meta.ncon] = ForwardDiff.jacobian(nlp.c, x) * v
    return Jv
end

function jtprod(nlp::HOADNLPModel, x::AbstractVector, v::AbstractVector)
    # increment!(nlp, :neval_jtprod)
    return ForwardDiff.jacobian(nlp.c, x)' * v
end

function jtprod!(
    nlp::HOADNLPModel,
    x::AbstractVector,
    v::AbstractVector,
    Jtv::AbstractVector,
)
    # increment!(nlp, :neval_jtprod)
    Jtv[1:nlp.meta.nvar] = ForwardDiff.jacobian(nlp.c, x)' * v
    return Jtv
end

function hess(
    nlp::HOADNLPModel,
    x::AbstractVector;
    obj_weight::Real = one(eltype(x)),
    y::AbstractVector = eltype(x)[],
)
    increment!(nlp, :neval_hess)
    Hx =
        obj_weight == 0.0 ? spzeros(nlp.meta.nvar, nlp.meta.nvar) :
        ForwardDiff.hessian(nlp.f, x) * obj_weight
    for i = 1:min(length(y), nlp.meta.ncon)
        if y[i] != 0.0
            Hx += ForwardDiff.hessian(x -> nlp.c(x)[i], x) * y[i]
        end
    end
    return tril(Hx)
end

function hess_coord(
    nlp::HOADNLPModel,
    x::AbstractVector;
    obj_weight::Real = one(eltype(x)),
    y::AbstractVector = eltype(x)[],
)
    Hx = hess(nlp, x, obj_weight = obj_weight, y = y)
    n = nlp.meta.nvar
    if VERSION < v"1.0"
        I = [(i, j, Hx[i, j]) for i = 1:n, j = 1:n if i ≥ j]
    else
        I = ((i, j, Hx[i, j]) for i = 1:n, j = 1:n if i ≥ j)
    end
    return (getindex.(I, 1), getindex.(I, 2), getindex.(I, 3))
end

function hprod(
    nlp::HOADNLPModel,
    x::AbstractVector,
    v::AbstractVector;
    obj_weight::Real = one(eltype(x)),
    y::AbstractVector = eltype(x)[],
)
    Hv = zeros(nlp.meta.nvar)
    return hprod!(nlp, x, v, Hv, obj_weight = obj_weight, y = y)
end

function hprod!(
    nlp::HOADNLPModel,
    x::AbstractVector,
    v::AbstractVector,
    Hv::AbstractVector;
    obj_weight::Real = one(eltype(x)),
    y::AbstractVector = eltype(x)[],
)
    increment!(nlp, :neval_hprod)
    n = nlp.meta.nvar
    Hv[1:n] =
        obj_weight == 0.0 ? zeros(nlp.meta.nvar) :
        ForwardDiff.hessian(nlp.f, x) * v * obj_weight
    for i = 1:min(length(y), nlp.meta.ncon)
        if y[i] != 0.0
            Hv[1:n] += ForwardDiff.hessian(x -> nlp.c(x)[i], x) * v * y[i]
        end
    end
    return Hv
end

function ∇f³xuv(nlp::HOADNLPModel, x::AbstractVector, u::AbstractVector, v::AbstractVector)
    ∇f(x) = ForwardDiff.gradient(nlp.f, x)
    # println("On  a ∇f")
    ϕᵤ(x) = ∇f(x)' * u
    # println("On  a ϕᵤ")
    ∇ϕᵤ(x) = ForwardDiff.gradient(ϕᵤ, x)
    # println("On  a ∇ϕᵤ")
    ψᵥ(x) = ∇ϕᵤ(x)' * v
    # println("On  a ψᵥ")
    ∇ψᵥ(x) = ForwardDiff.gradient(ψᵥ, x)
    # println("On  a ∇ψᵥ")
    # increment!(nlp, :neval_Tuvprod)
    return ∇ψᵥ(x)
end

function ∇f³xuv2(nlp::HOADNLPModel, x::AbstractVector, u::AbstractVector, v::AbstractVector)
    ∇²fu(x) = ForwardDiff.hessian(nlp.f, x) * u
    ∇²fuv(x) = ∇²fu(x)' * v
    return ForwardDiff.gradient(∇²fuv, x)
end


text:
# All functions needed for the gradient/variable environment experiments.
# Randomly change habitat precipitation and temperature.

"""
    changetemp!(world, sdtemp)

Change the temperature of all patches according to a normal distribution.
"""
function changetemp!(world::Array{Patch,1}, sdtemp::Float64)
    sdtemp == 0 && return
    deltaval = rand(Normal(0.0, sdtemp))
    for patch in world
        patch.temp += deltaval
    end
    markthem!(world)
end

"""
    changeprec!(world, sdprec)

Change the precipitation of all patches according to a normal distribution.
"""
function changeprec!(world::Array{Patch,1}, sdprec::Float64)
    sdprec == 0 && return
    deltaval = rand(Normal(0.0, sdprec))
    for patch in world
        patch.prec += deltaval
    end
    markthem!(world)
end

"""
    changehabitat!(world)

Carry out 'global change' on all patches.
"""
function changehabitat!(world::Array{Patch,1})
    # TODO: record trajectory? input trajectory?
    changetemp!(world, setting("sdtemp"))
    changeprec!(world, setting("sdprec"))
end

"""
    markthem!(community)

Set each individual in the community (= array of individuals) as "marked".
"""
function markthem!(community::Array{Individual, 1})
    for ind in community
        ind.marked = true
    end
end

"""
    markthem!(habitat)

Set each individual in the given patch as "marked".
"""
function markthem!(habitat::Patch)
    markthem!(habitat.community)
end

"""
    markthem!(world)

Set every individual in the world as "marked".
"""
function markthem!(world::Array{Patch, 1})
    for habitat in world
        markthem!(habitat)
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
#module HT

if debug
const LIBHTS = "/haplox/users/guo/Github/htslib/libhts.so"
const LIBKSTR = joinpath(dirname(@__FILE__), "c/libkstr.so")
end

macro fncall(fn, return_type)
    f = eval(fn)
    quote
        ccall(($(Meta.quot(f)), $LIBHTS), $return_type,
              ())
    end
end

macro fncall(fn, return_type, argtypes, args...)
    f = eval(fn)
    args = map(esc, args)
    quote
        ccall(($(Meta.quot(f)), $LIBHTS), $return_type,
              $argtypes, $(args...))
    end
end

macro kstrcall(fn, return_type)
    f = eval(fn)
    quote
        ccall(($(Meta.quot(f)), $LIBKSTR), $return_type,
              ())
    end
end

include("type.jl")
include("fn.jl")
include("io.jl")
include("utils.jl")

#end


text:
function GaussianRBM(θv::AbstractArray, γv::AbstractArray, θh::AbstractArray, γh::AbstractArray, w::AbstractArray)
    @assert size(θv) == size(γv)
    @assert size(θh) == size(γh)
    @assert size(w) == (size(θv)..., size(θh)...)
    return RBM(Gaussian(θv, γv), Gaussian(θh, γh), w)
end

function GaussianRBM(θv::AbstractArray, γv::AbstractArray, w::AbstractArray)
    @assert size(θv) == size(γv) == size(w)[1:ndims(θv)]
    θh = fill!(similar(θv, eltype(θv), size(w)[(ndims(θv) + 1):end]),  0)
    γh = fill!(similar(γv, eltype(γv), size(w)[(ndims(γv) + 1):end]),  1)
    return GaussianRBM(θv, γv, θh, γh, w)
end

# For a Gaussian-Gaussian RBM we can use the analytical expression
function log_partition(rbm::RBM{<:Gaussian, <:Gaussian})
    θ = [vec(visible(rbm).θ); vec(hidden(rbm).θ)]
    γv = vec(abs.(visible(rbm).γ))
    γh = vec(abs.(hidden(rbm).γ))
    w = flat_w(rbm)

    lA = block_matrix_logdet(
        Diagonal(γv), -w,
        -w', Diagonal(γh)
    )

    iA = block_matrix_invert(
        Diagonal(γv), -w,
        -w', Diagonal(γh)
    )

    return length(θ)/2 * log(2π) + dot(θ, iA, θ)/2 - lA/2
end


text:
abstract type pseudoGPUArray <: SolidStateDetectors.GPUArrays.AbstractGPUArray{T, N where N} end
SolidStateDetectors.get_device(::Type{pseudoGPUArray}) = SolidStateDetectors.KernelAbstractions.CPU()

@testset "Simulate example detector: Inverted Coax" begin
    sim = Simulation{T}(SSD_examples[:InvertedCoax])

    calculate_electric_potential!(
        sim,
        device_array_type = pseudoGPUArray,
        convergence_limit = 0, # For this we have to probably add an addition kernel function
        max_n_iterations = 2000, # For now just perform always `max_n_iterations` iterations
        refinement_limits = [0.2, 0.1],
        depletion_handling = true
    )

    W = SolidStateDetectors.calculate_stored_energy(sim)
    @test isapprox(ustrip(W), 2.41e-4, atol = 1e-5)
end

@testset "Simulate example detector: Inverted Coax" begin
    sim = Simulation{T}(SSD_examples[:CGD])

    calculate_electric_potential!(
        sim,
        device_array_type = pseudoGPUArray,
        convergence_limit = 0, # For this we have to probably add an addition kernel function
        max_n_iterations = 2000, # For now just perform always `max_n_iterations` iterations
        refinement_limits = [0.2],
        depletion_handling = true
    )

    W = SolidStateDetectors.calculate_stored_energy(sim)
    @test isapprox(ustrip(W), 4.39e-6, atol = 1e-5)
end


text:
using SafeTestsets

@safetestset "Gradients                                                                       " begin include("gradient_tests.jl") end
@safetestset "Jacobians                                                                       " begin include("jacobian_tests.jl") end
@safetestset "Hessians                                                                        " begin include("hessian_tests.jl") end
@safetestset "Objectives                                                                      " begin include("objectives_tests.jl") end
@safetestset "Linear Solvers                                                                  " begin include("linear_solvers_tests.jl") end
@safetestset "Line Searches                                                                   " begin include("line_searches_tests.jl") end
@safetestset "Nonlinear Solvers                                                               " begin include("nonlinear_solvers_tests.jl") end
@safetestset "Optimizers                                                                      " begin include("optimizers_tests.jl") end


text:
"Square the sum of the first `n` positive integers"
function square_of_sum(n)
  return sum(1:n) * sum(1:n)
end

"Sum the squares of the first `n` positive integers"
function sum_of_squares(n)
  return sum(map(x -> x * x, 1:n))
end

"Subtract the sum of squares from square of the sum of the first `n` positive ints"
function difference(n)
  return square_of_sum(n) - sum_of_squares(n)
end


text:
# Accessor functions to parse CellML xml(doc) files
# naming convention:
#   list_*:         returns a list of 'Ezxml(doc).Node's
#   get_*:          returns a single 'Ezxml(doc).Node'
#   list_*_names:   returns a list of String

EzXML.root(doc::Document) = root(doc.xmls[1])

cellml_ns(doc::Document) = cellml_ns(doc.xmls[1])
cellml_ns(comp::Component) = cellml_ns(comp.node)

nodeof(xml::EzXML.Document) = root(xml)
nodeof(node::EzXML.Node) = node
nodeof(comp::Component) = comp.node

"""
    list_components returns the list of CellML <Component>s
"""
list_components(xml) = findall("//x:model/x:component", root(xml), ["x"=>cellml_ns(xml)])

"""
    get_model returns the single <model> element of a CellML file
"""
get_model(xml) = findfirst("//x:model", root(xml), ["x"=>cellml_ns(xml)])

"""
    get_component_variables returns the list of the variables of a component
    comp is an Ezxml(doc) as returned by list_components
"""
list_component_variables(comp) = findall("./x:variable", nodeof(comp), ["x"=>cellml_ns(comp)])

"""
    list_initiated_variables returns all variables that have an initial_value
"""
# list_initiated_variables(xml::EzXML.Document) = findall("//x:component/x:variable[@initial_value]", root(xml), ["x"=>cellml_ns(xml)])

list_initiated_variables(comp::Component) = findall("./x:variable[@initial_value]", nodeof(comp), ["x"=>cellml_ns(comp)])

"""
    list_connections returns the list of <connection> nodes in the CellML document
"""
list_connections(xml) = findall("//x:connection", root(xml), ["x"=>cellml_ns(xml)])

"""
    get_connection_variables returns the pair of components for the given connection
    conn is an Ezxml(doc) node as returned by list_connections
"""
get_connection_component(conn) = findfirst("./x:map_components", nodeof(conn), ["x"=>cellml_ns(conn)])

"""
    components_of converts the output of get_connection_component to a pair
    of strings (the names of the components)
"""
components_of(x) = (x["component_1"], x["component_2"])

"""
    list_connection_variables returns a list composed of pairs of variable
    for the given connection
    conn is an Ezxml(doc) node as returned by list_connections
"""
list_connection_variables(conn) = findall("./x:map_variables", nodeof(conn), ["x"=>cellml_ns(conn)])

"""
    components_of converts the output of list_connection_variables to a pair
    of strings (the names of the variables)
"""
variables_of(x) = (x["variable_1"], x["variable_2"])

"""
    list_state_names returns the names of variables that occur on the left-hand-side
    of an ODE; hence are *state* variables
    comp is an Ezxml(doc) as returned by list_components
"""
function find_state_names(comp)
    nodes = findall("./y:math/y:apply", nodeof(comp), ["y"=>mathml_ns])
    names = String[]
    for n in nodes
        e = elements(n)
        if e[1].name == "eq" && e[2].name == "apply"
            h = elements(e[2])
            if h[1].name == "diff" && h[3].name == "ci"
                push!(names, strip(nodecontent(h[3])))
            end
        end
    end
    return names
end

"""
    list_alg_names returns the names of variables that occur on the left-hand-side
    of a non-ODE equation; hence are *algebraic* variables
    comp is an Ezxml(doc) as returned by list_components
"""
function find_alg_names(comp)
    nodes = findall("./y:math/y:apply", nodeof(comp), ["y"=>mathml_ns])
    names = String[]
    for n in nodes
        e = elements(n)
        if e[1].name == "eq" && e[2].name == "ci"
            push!(names, strip(nodecontent(e[2])))
        end
    end
    return names
end

"""
    list_component_math returns a list of the math elements in the
    given component
"""
list_component_math(comp) = findall("./y:math", nodeof(comp), ["y"=>mathml_ns])

list_component_bvar(comp) = findall(".//y:math//y:bvar/y:ci", nodeof(comp), ["y"=>mathml_ns])

"""
    list_imports returns the list of <import> nodes in the CellML document
"""
list_imports(xml) = findall("//x:import", get_model(xml), ["x"=>cellml_ns(xml)])

"""
    list_import_components returns the list of component nodes of the given
    import element
    node: an import node as returned by list_imports
"""
list_import_components(node) = findall("./x:component", nodeof(node), ["x"=>cellml_ns(node)])

function list_encapsulation(doc, comp)
    name = string(nameof(comp))
    ns = cellml_ns(doc)
    groups = parentnode.(findall("//x:group/x:relationship_ref[@relationship='encapsulation']", root(doc), ["x"=>ns]))
    if isempty(groups)  # CellML ver 2.0
        return findall("//x:encapsulation//x:component_ref[@component='$name']//x:component_ref", root(doc), ["x"=>ns])
    else                # CellML ver 1.0 and 1.1
        for g in groups
            nodes = findall(".//x:component_ref[@component='$name']//x:component_ref", g, ["x"=>ns])
            if !isempty(nodes)
                return nodes
            end
        end
        return EzXML.Node[]
    end
end


text:
function assess_pheno_damage(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::POOLConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,maximum(damage_times)+199))

    evaluate_individual_view!(damage_times,damage_percentage,p_float,p_int,cell_config,exp_config,grn_config,results)

    o =  evaluate_fitness(results[:,:,evo_config.time_steps],evo_config)
    d = evaluate_fitness(results[:,:,maximum(damage_times)],evo_config)
    r = evaluate_fitness(results[:,:,maximum(damage_times) + 49],evo_config)

    return o,d,r

end

function assess_pheno_damage(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::EMAPConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,maximum(damage_times)+199))

    evaluate_individual_view!(damage_times,damage_percentage,p_float,p_int,cell_config,exp_config,grn_config,results)

    o =  evaluate_fitness(results[:,:,evo_config.time_steps],evo_config)[1]
    d = evaluate_fitness(results[:,:,maximum(damage_times)],evo_config)[1]
    r = evaluate_fitness(results[:,:,maximum(damage_times) + 49],evo_config)[1]

    return o,d,r

end

function assess_genome_damage(dp::Int,p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::POOLConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    bi = copy(p_int)
    bf = copy(p_float)

    for _ in 1:dp
        bi,bf = deletion_mutation(bf,bi,evo_config.genome_preserve,1.)
    end

    original = evaluate_individual_fitness(p_float,p_int,evo_config,cell_config,exp_config,grn_config)

    damage = evaluate_individual_fitness(bf,bi,evo_config,cell_config,exp_config,grn_config)

    return original,damage
end

function assess_genome_damage(dp::Int,p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::EMAPConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    bi = copy(p_int)
    bf = copy(p_float)

    for _ in 1:dp
        bi,bf = deletion_mutation(bf,bi,evo_config.genome_preserve,1.)
    end

    original = evaluate_individual_fitness(p_float,p_int,evo_config,cell_config,exp_config,grn_config)[1]

    damage = evaluate_individual_fitness(bf,bi,evo_config,cell_config,exp_config,grn_config)[1]

    return original,damage
end

function assess_stability(times::Vector{Int64},p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::POOLConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,1500))

    evaluate_individual_view!(p_float,p_int,cell_config,exp_config,grn_config,results)

    r = zeros(size(times,1) + 1)

    idx = 1

    r[idx] = evaluate_fitness(results[:,:,800],evo_config)

    for t in times
        idx+=1
        r[idx] =  evaluate_fitness(results[:,:,t],evo_config)
    end

    return r[1],r[2],r[3],r[4]

end

function assess_stability(times::Vector{Int64},p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::EMAPConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,1500))

    evaluate_individual_view!(p_float,p_int,cell_config,exp_config,grn_config,results)

    r = zeros(size(times,1) + 1)

    idx = 1

    r[idx] = evaluate_fitness(results[:,:,800],evo_config)[1]

    for t in times
        idx+=1
        r[idx] =  evaluate_fitness(results[:,:,t],evo_config)[1]
    end

    return r[1],r[2],r[3],r[4]

end

function compute_system_energy(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::POOLConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    energies = zeros(Float64,evo_config.time_steps + 50);
    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,evo_config.time_steps + 50));

    evaluate_individual_view!(damage_times,damage_percentage,p_float,p_int,cell_config,exp_config,grn_config,energies,results)
    fitness = zeros(evo_config.time_steps + 50)

    for step in 1:size(results,3)
        fitness[step] = evaluate_fitness(results[:,:,step],evo_config)
    end

    return energies,fitness
end

function compute_system_energy(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,evo_config::EMAPConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    energies = zeros(Float64,evo_config.time_steps + 50);
    results = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,evo_config.time_steps + 50));

    evaluate_individual_view!(damage_times,damage_percentage,p_float,p_int,cell_config,exp_config,grn_config,energies,results)
    fitness = zeros(evo_config.time_steps + 50)

    for step in 1:size(results,3)
        fitness[step] = evaluate_fitness(results[:,:,step],evo_config)[1]
    end

    return energies,fitness
end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
module FinOptions

export

        Black76,
        BlackScholes,
        GBlackScholes,
        KirkSpreadOption,
        GeoAveRate,
        GDelta,
        Delta,
        Elasticity,
        Gamma,
        DGammaDVol,
        GammaP,
        DDeltaDVol,
        Vega,
        Vomma,
        VegaP,
        DVegaDVol,
        Theta,
        Rho,
        RhoFut,
        Rho2,
        Carry,
        Speed,
        StrikeDelta,
        RNDensity,
        PutCallParity,
        VannaVolga,
        SwingOptionPricer

include("Black76.jl")
include("BlackScholes.jl")
include("CND.jl")
#include("CNDEV.jl"),
include("GBlackScholes.jl")
include("GDelta.jl")
include("Greeks.jl")
include("PutCallParity.jl")
include("VannaVolga.jl")
include("SwingOptionPricer.jl")

# package code goes here

end # module


text:
module SimpleMesh

include("SimpleMeshFunction2.jl")
include("../../../SummationByParts/src/SummationByParts.jl")

using SummationByParts

export simpleMesh

@doc """
### SimpleMesh

simpleMesh is a mesh type which generates a rectangular mesh in 2D with triangular elements. It has the following fields

**Fields**

*  `lengthx` : Length of domain on the X-axis
*  `lengthy` : Length of domain on the Y-axis
*  `nedx`    : Number of element edges along X-axis
*  `nedy`    : Number of element edges along Y-axis
*  `nnpe`    : Number of nodes per element edge (Same for all edges)
*  `order`   : Order of elements (linear, quadratic, cubic etc.)
*  `numNodes`  : Number of nodal points
*  `numVert`   : Total number of vertices in the mesh
*  `numEdge    : Total number of edges in the mesh
*  `numEl`     : Number of elements mesh
*  `numDof`    : Total degrees of freedom in the system
*  `numDofPerNode` : number of DOFs per node. Can be thought of number of equations being solved at a node
*  `numNodesPerElement` : Number of nodes per element
*  `numBoundaryEdges`   : Number of element edges on the boundary
*  `dofs`     : Element information matrix.
*  `vtxCoord` : Location coordinates of the vertices
*  `elemEdgeLengthx` : Length of element edge along the X-axis
*  `elemEdgeLengthy` : Length of element edge along the Y-axis
*  `bndryfaces` : An array of boundaries lying on the grid. It is of the type Boundary
*  `interfaces` : An array of all interior edges on the grid. It is of the type Interface

"""->

mutable struct simpleMesh{T <: FloatingPoint}

  lengthx::FloatingPoint
  lengthy::FloatingPoint
  nedx::Int
  nedy::Int
  nnpe::Int
  order::Int
  numNodes::Int
  numVert::Int
  numEdge::Int
  numEl::Int
  numDof::Int
  numDofPerNode::Int
  numNodesPerElement::Int
  numBoundaryEdges::Int
  dofs::Array{Int}
  vtxCoord::Array{T}
  elemEdgeLengthx::FloatingPoint
  elemEdgeLengthy::FloatingPoint
  bndryfaces::Array{Boundary}
  interfaces::Array{Interface}

  # dofs::Array{Int,3} # Store dof numbers of solution array for speedy assembly

  # numBC::Int  # number of boundary conditions

  function simpleMesh{T <: FloatingPoint}(lengthx,lengthy,nedx,nedy,nnpe,numDofPerNode) where T <: FloatingPoint
  order = nnpe-1 # Calculate the mesh order

  # calculating the members of the type using functions
    (dofs,numNodes,numEl,numNodesPerElement) = createMesh(lengthx,lengthy,nedx,nedy,nnpe,numDofPerNode)

  numVert = (nedx+1)*(nedy+1)
  numDof = numNodes * numDofPerNode
  numBoundaryEdges = 2*(nedx+nedy)

  (vtxCoord, elemEdgeLengthx, elemEdgeLengthy) = vtxLocation(lengthx, lengthy, nedx, nedy, numEl)

    (bndryfaces, numEdge) = boundaryEdges(nedx, nedy, numEl)
    interfaces = interiorEdgeinfo(nedx,nedy,numEl)

  new(lengthx,lengthy,nedx,nedy,nnpe,order,numNodes,numVert,numEdge,numEl,
    numDof,numDofPerNode,numNodesPerElement,numBoundaryEdges,dofs,vtxCoord,
    elemEdgeLengthx,elemEdgeLengthy,bndryfaces,interfaces)

  end # Ends the function simpleMesh
end # Ends the type simpleMesh

end # Ends the module


text:
using Base.Test
using SpecialFunctions, NaNMath
using DiffRules

srand(1)

function finitediff(f, x)
    ϵ = cbrt(eps(typeof(x))) * max(one(typeof(x)), abs(x))
    return (f(x + ϵ) - f(x - ϵ)) / (ϵ + ϵ)
end

non_numeric_arg_functions = [(:Base, :rem2pi, 2)]

for (M, f, arity) in DiffRules.diffrules()
    (M, f, arity) ∈ non_numeric_arg_functions && continue
    if arity == 1
        @test DiffRules.hasdiffrule(M, f, 1)
        deriv = DiffRules.diffrule(M, f, :x)
        modifier = in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)) ? 1 : 0
        @eval begin
            x = rand() + $modifier
            @test isapprox($deriv, finitediff($M.$f, x), rtol=0.05)
        end
    elseif arity == 2
        @test DiffRules.hasdiffrule(M, f, 2)
        derivs = DiffRules.diffrule(M, f, :x, :y)
        @eval begin
            x, y = rand(1:10), rand()
            dx, dy = $(derivs[1]), $(derivs[2])
            if !(isnan(dx))
                @test isapprox(dx, finitediff(z -> $M.$f(z, y), float(x)), rtol=0.05)
            end
            if !(isnan(dy))
                @test isapprox(dy, finitediff(z -> $M.$f(x, z), y), rtol=0.05)
            end
        end
    end
end

# Treat rem2pi separately because of its non-numeric second argument:
derivs = DiffRules.diffrule(:Base, :rem2pi, :x, :y)
for xtype in [:Float64, :BigFloat, :Int64]
    for mode in [:RoundUp, :RoundDown, :RoundToZero, :RoundNearest]
        @eval begin
            x = $xtype(rand(1 : 10))
            y = $mode
            dx, dy = $(derivs[1]), $(derivs[2])
            @test isapprox(dx, finitediff(z -> rem2pi(z, y), float(x)), rtol=0.05)
            @test isnan(dy)
        end
    end
end