text:
module EvoTrees

export init_evotree, grow_evotree!, grow_tree, fit_evotree, predict,
    EvoTreeRegressor, EvoTreeCount, EvoTreeClassifier, EvoTreeGaussian,
    EvoTreeRModels, importance

using Statistics
using Base.Threads: @threads
using StatsBase: sample, quantile
using Random: seed!
using StaticArrays
using Distributions
using CategoricalArrays
import MLJModelInterface
import MLJModelInterface: fit, update
import MLJModelInterface: predict

include("models.jl")
include("structs.jl")
include("loss.jl")
include("eval.jl")
include("predict.jl")
include("find_split.jl")
include("fit.jl")
include("importance.jl")
include("MLJ.jl")

end # module


text:
function MaximalOrder(K::NfAbsNS; discriminant::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[])
  try
    c = _get_maximal_order(K)::NfAbsOrd{NfAbsNS, NfAbsNSElem}
    return c
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    O = maximal_order_from_components(K)
    O.ismaximal = 1
    _set_maximal_order(K, O)
    return O
  end
end

###############################################################################
#
#  Generic code for orders
#
###############################################################################

function new_maximal_order(O::NfAbsOrd{S, T}; index_divisors::Vector{fmpz} = fmpz[], disc::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[]) where {S, T}
  return maximal_order_round_four(O, index_divisors= index_divisors, disc = disc, ramified_primes = ramified_primes)
end

function maximal_order_round_four(O::NfAbsOrd; index_divisors::Vector{fmpz} = fmpz[], disc::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[])
  OO = O
  M = trace_matrix(O)
  l = divisors(M, discriminant(O))
  if !isempty(index_divisors)
    push!(l, index_divisors)
  end
  if !isempty(ramified_primes)
    push!(l, ramified_primes)
  end
  l = coprime_base(l)
  for s in l
    if disc != -1
      u = divexact(discriminant(OO), disc)
      if isone(gcd(u, s))
        continue
      end
    end
    @vtime :NfOrd fac = factor(s)
    for (p, j) in fac
      @vprint :NfOrd 1 "Computing p-maximal overorder for $p ..."
      O1 = pmaximal_overorder(O, p)
      if valuation(discriminant(O1), p) < valuation(discriminant(OO),p)
        OO += O1
      end
      @vprint :NfOrd 1 "done\n"
    end
  end
  OO.ismaximal = 1
  return OO
end

function maximal_order_from_components(L::NfAbsNS; disc::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[])
  BKs, lp = maximal_order_of_components(L)
  B = product_basis(BKs)
  OO = Order(L, B, check = false, cached = false, isbasis = true)
  if disc != -1 && discriminant(OO) == disc
    return OO
  end
  if !isempty(ramified_primes)
    append!(lp, ramified_primes)
  end
  lp = coprime_base(lp)
  for p in lp
    if isprime(p)
      OO = pmaximal_overorder(OO, p)
    else
      fac = factor(p)
      for (k, v) in fac
        OO = pmaximal_overorder(OO, k)
      end
    end
  end
  OO.ismaximal = 1
  return OO

end

function product_basis(l::Vector{Vector{T}}) where T <: Union{NfAbsOrdElem, nf_elem, NfAbsNSElem}
  nelems = 1
  for i = 1:length(l)
    nelems *= length(l[i])
  end
  B = Vector{T}(undef, nelems)
  ind = length(l[1])
  for i = 1:ind
    B[i] = l[1][i]
  end
  for jj = 2:length(l)
    new_deg = length(l[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = B[j] * l[jj][i]
      end
    end
    ind *= new_deg
  end
  return B
end

function product_basis(l1::Vector{T}, l2::Vector{T}) where T <: Union{NfAbsOrdElem, nf_elem, NfAbsNSElem}
  B = Vector{T}(undef, length(l1)*length(l2))
  for i = 1:length(l1)
    B[i] = l1[i]
  end
  for i = 1:length(l2)
    for j = 1:length(l1)
      B[(i-1)* length(l1) + j] = B[j] * l2[i]
    end
  end
  return B
end

function maximal_order_of_components(L::NfAbsNS) where {S, T}
  Qx, x = PolynomialRing(FlintQQ, "x")
  fields = Vector{AnticNumberField}(undef, length(L.pol))
  for i = 1:length(L.pol)
    f = Qx(L.pol[i])
    K = NumberField(f, cached = false)[1];
    OK = maximal_order(K)
    fields[i] = K
  end
  mvpolring = parent(L.pol[1])
  gpols = gens(mvpolring)
  #Now, bring the maximal order of every component in L
  B = Vector{Vector{NfAbsNSElem}}(undef, length(fields))
  for i = 1:length(fields)
    OK = maximal_order(fields[i])
    BOK = OK.basis_nf
    BK = Vector{NfAbsNSElem}(undef, degree(OK))
    for j = 1:length(BK)
      polel = Qx(BOK[j])
      polm = evaluate(polel, gpols[i])
      BK[j] = L(polm)
    end
    B[i] = BK
  end
  lp = fmpz[]
  for i = 1:length(fields)
    for j = i+1:length(fields)
      push!(lp, gcd(discriminant(maximal_order(fields[i])), discriminant(maximal_order(fields[j]))))
    end
  end
  return B, lp
end

ext: jl


text:
module signalobject

include("utils/utils.jl")
include("config/ecgconfig.jl")

import Base

import .ecgconfig as cfg
using .utils

export SignalObject
mutable struct SignalObject{T<:Real}
    signals::Vector{T}
    signal_mean_absolute::Union{Vector{T}, UndefInitializer}
    signal_second_derivative::Union{Vector{T}, UndefInitializer}
    signal_enhance::Union{Vector{T}, UndefInitializer}
    SignalObject(signals::Vector{T}) where T = new{T}(signals)
end

function signal_mean_absolute(signals::Vector{<:Real}, windows_msec::Real = 20)

    windows = msec_to_sample(windows_msec, cfg.FS)
    return  utils.ma(abs.(signals), windows)
end

function signal_second_derivative(signals::Vector{T}) where T
    second_dif_signals = zeros(T, length(signals))
    second_dif_signals[2:end-1] = diff(diff(signals))
    return second_dif_signals
end

function Base.getproperty(sig::SignalObject, s::Symbol)
    if s == :signal_mean_absolute
        if isdefined(sig, s)
            return getfield(sig, s)
        else
            setfield!(sig, s, signal_mean_absolute(sig.signals))
            return getfield(sig, s)
        end
    elseif s == :signal_enhance
        if isdefined(sig, s)
            return getfield(sig, s)
        else
            return sig.signals
        end
    elseif s == :signal_second_derivative
        if isdefined(sig, s)
            return getfield(sig, s)
        else
            setfield!(sig, s, signal_second_derivative(sig.signals))
            return getfield(sig, s)
        end
    else
        return getfield(sig, s)
    end
end

Base.length(sig::SignalObject) = length(sig.signals)

end


text:
#=
Copyright (c) 2016, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module Checkpointing

import ParallelAccelerator

using CompilerTools
import CompilerTools.DebugMsg
DebugMsg.init()

import HPAT

# information about AST gathered and used in Checkpointing
type DomainState
    LambdaVarInfo :: CompilerTools.LambdaHandling.LambdaVarInfo
end

single_node_mttf = 4000.0   # default to 4000 hours MTTF of a single node
single_node_faults_per_million = 1000000.0 / single_node_mttf
checkpoint_time = 1.0 / 3600.0  # 1 minute checkpoint time converted to hours

@noinline function hpat_get_checkpoint_time(checkpoint_location)
    1.0
end

@noinline function hpat_get_sec_since_epoch()
    convert(Float64,1)
end

# Starts a checkpoint session.
# Return a Int64 checkpoint handle that should be passed to the other checkpoint routines.
@noinline function hpat_start_checkpoint(checkpoint_location)
    convert(Int32,1)
end

@noinline function hpat_value_checkpoint(checkpoint_handle :: Int32, value)
    convert(Int32,1)
end

# Finishes a checkpointing session.
# Takes the handle of the checkpoint session to terminate.
@noinline function hpat_end_checkpoint(checkpoint_handle :: Int32)
    convert(Int32,1)
end

@noinline function hpat_finish_checkpoint_region(checkpoint_location)
    convert(Int32,1)
end

# Start of checkpointing restore functions.
# Prevent hoist will be removed in the future once ParallelIR reordering issue is fixed.
@noinline function hpat_checkpoint_restore_start(checkpoint_location, prevent_hoist :: ANY)
    convert(Int32,1)
end

@noinline function hpat_checkpoint_restore_value(checkpoint_handle :: Int32, value)
    convert(Int32,1)
end

@noinline function hpat_checkpoint_restore_end(checkpoint_handle :: Int32)
    convert(Int32,1)
end
# End of checkpointing restore functions.

@noinline function hpat_dist_node_id()
    convert(Int32,1)
end

# Set the delay between checkpoints (in seconds).
# "0" means use Young's formula.
checkpoint_debug = 0
function setCheckpointDebug(checkpoint_delay)
    global checkpoint_debug = checkpoint_delay
end

# ENTRY to checkpointing
function from_root(function_name, ast :: Expr, with_restart :: Bool)
    @assert ast.head == :lambda "Input to Checkpointing should be :lambda Expr"
    @dprintln(1,"Starting main Checkpointing.from_root.  function = ", function_name, " ast = ", ast, " with_restart = ", with_restart)

    state::DomainState = DomainState(CompilerTools.LambdaHandling.lambdaExprToLambdaVarInfo(ast))

    body = CompilerTools.LambdaHandling.getBody(ast)

    lives = ParallelAccelerator.ParallelIR.computeLiveness(ast, state.LambdaVarInfo)
    @dprintln(3,"lives = ", lives)
    loop_info = CompilerTools.Loops.compute_dom_loops(lives.cfg)
    @dprintln(3,"loop_info = ", loop_info)

    loops = loop_info.loops

    if length(loops) != 1
        @dprintln(0,"Checkpointing.from_root currently only supports functions with exactly one loop.  ", function_name, " has ", length(loops), " loops.")
        return ast
    end
    for loop_index = 1:length(loops)
        the_loop = loops[loop_index]

        loop_entry = the_loop.head
        loop_entry_bb = lives.cfg.basic_blocks[loop_entry]
        @dprintln(3,"loop_entry = ", loop_entry)
        liveness_loop_entry = CompilerTools.LivenessAnalysis.getBasicBlockFromBlockNumber(loop_entry, lives)
        @dprintln(3,"liveness_loop_entry = ", liveness_loop_entry)
        bb_loop_members = map(x -> CompilerTools.LivenessAnalysis.getBasicBlockFromBlockNumber(x, lives), the_loop.members)
        @dprintln(3,"bb_loop_members = ", bb_loop_members)
        loop_live_in = liveness_loop_entry.live_in
        @dprintln(3,"loop_live_in = ", loop_live_in)
        loop_def = reduce((x,y) -> union(x,y), Set{LHSVar}(), map(x -> x.def, bb_loop_members))
        @dprintln(3,"loop_def = ", loop_def)
        live_in_and_def = intersect(loop_live_in, loop_def)
        @dprintln(3,"live_in_and_def = ", live_in_and_def)
        assert(!isempty(live_in_and_def))

        pre_loop_stmts = Any[]
        checkpoint_timer_sn = ParallelAccelerator.ParallelIR.createStateVar(state, "__hpat_checkpoint_timer", Float64, ParallelAccelerator.ParallelIR.ISASSIGNED)

        if with_restart
           restore_handle = ParallelAccelerator.ParallelIR.createStateVar(state, string("__hpat_restore_handle_", loop_index), Int32, ParallelAccelerator.ParallelIR.ISASSIGNED)
           liad_array = [live_in_and_def...]
           push!(pre_loop_stmts, ParallelAccelerator.ParallelIR.mk_assignment_expr(restore_handle, ParallelAccelerator.ParallelIR.TypedExpr(Int32, :call, GlobalRef(HPAT.Checkpointing,:hpat_checkpoint_restore_start), loop_index, liad_array[1])))
           for i = 1:length(liad_array)
               push!(pre_loop_stmts, ParallelAccelerator.ParallelIR.TypedExpr(Int32, :call, GlobalRef(HPAT.Checkpointing,:hpat_checkpoint_restore_value), restore_handle, liad_array[i]))

               @dprintln(3,"liad_array[i] = ", liad_array[i], " type = ", typeof(liad_array[i]))
               latyp = CompilerTools.LambdaHandling.getType(liad_array[i], state.LambdaVarInfo)
               @dprintln(3,"latyp = ", latyp)
               if CompilerTools.Helper.isArrayType(latyp)
                   @dprintln(3,"isArrayType")
                   push!(pre_loop_stmts, Expr(:call, GlobalRef(HPAT,:__hpat_dist_broadcast), liad_array[i], ParallelAccelerator.ParallelIR.TypedExpr(Int64, :call, GlobalRef(Base,:arraylen), liad_array[i])))
               else
                   @dprintln(3,"!isArrayType")
                   push!(pre_loop_stmts, Expr(:call, GlobalRef(HPAT,:__hpat_dist_broadcast), liad_array[i], 1))
               end
           end
           push!(pre_loop_stmts, ParallelAccelerator.ParallelIR.TypedExpr(Int32, :call, GlobalRef(HPAT.Checkpointing,:hpat_checkpoint_restore_end), restore_handle))
        end
        push!(pre_loop_stmts, ParallelAccelerator.ParallelIR.mk_assignment_expr(checkpoint_timer_sn, ParallelAccelerator.ParallelIR.TypedExpr(Float64, :call, GlobalRef(HPAT.Checkpointing,:hpat_get_sec_since_epoch))))

        CompilerTools.Loops.insertNewBlockBeforeLoop(the_loop, lives.cfg, pre_loop_stmts)
        @dprintln(3,"CFG after insert = ", lives.cfg)

        # Add a call after the loop to clean up any checkpoint file.
        postSet = CompilerTools.Loops.getPostLoopBlock(the_loop, lives.cfg)
        assert(length(postSet) == 1)
        postLoopBlockIndex = first(postSet)
        @dprintln(3,"postLoopBlockIndex = ", postLoopBlockIndex)
        @dprintln(3,"cfg = ", lives.cfg.basic_blocks)
        postLoopBB = lives.cfg.basic_blocks[postLoopBlockIndex]
        if checkpoint_debug == 0
          CompilerTools.CFGs.insertStatementBeginningOfBlock(lives.cfg, postLoopBB, ParallelAccelerator.ParallelIR.TypedExpr(Int32, :call, GlobalRef(HPAT.Checkpointing,:hpat_finish_checkpoint_region), loop_index))
        end

        # Create the checkpoint function as a string and then parse/eval to force it into existence.
        # The function takes the last checkpoint time.  If enough time has expired then do the checkpoint and return the current time.
        # If checkpoint time has not arrived then return the last checkpoint time.
        # The time between checkpoints is sqrt( 2 * time_to_checkpoint * system_mttf).
        # We get system_mttf by assuming some reasonable single_node_mttf, converting to failures per million hours, multiple by the
        # number of nodes in the system to get full system failures per million hours and then convert back to full system mttf (in hours).
        checkpoint_func_name = string("__hpat_checkpoint_func_", loop_index)

        # Creates the first line of the function with these characteristics.
        #   1. The function name is __hpat_checkpoint_func_ followed by a unique number to make the function name unique.
        #   2. The first parameter is the last checkpoint time.
        #   3. Other parameters are the elements that need to go in the checkpoint file.
        argument_names = [ string("arg",i) for i = 1:length(live_in_and_def) ]
        checkpoint_func_str = string(                     "function ", checkpoint_func_name, "(start_time, num_pes, ", foldl((a,b) -> "$a, $b", argument_names), ")\n")
        checkpoint_func_str = string(checkpoint_func_str, "    node_id = HPAT.Checkpointing.hpat_dist_node_id()\n")
        checkpoint_func_str = string(checkpoint_func_str, "    if node_id != 0\n")
        checkpoint_func_str = string(checkpoint_func_str, "        return start_time\n")
        checkpoint_func_str = string(checkpoint_func_str, "    end\n")
        checkpoint_func_str = string(checkpoint_func_str, "    system_faults_per_million_hours = num_pes * ", single_node_faults_per_million, "\n")
        checkpoint_func_str = string(checkpoint_func_str, "    system_mttf::Float64 = 1000000.0 / system_faults_per_million_hours\n")
        checkpoint_func_str = string(checkpoint_func_str, "    cur_time = HPAT.Checkpointing.hpat_get_sec_since_epoch()\n")
        if checkpoint_debug != 0
        checkpoint_func_str = string(checkpoint_func_str, "    if ((cur_time - start_time) > ", checkpoint_debug, ")\n")
        else
        checkpoint_func_str = string(checkpoint_func_str, "    if ((cur_time - start_time) / 3600.0) > sqrt(2 * system_mttf * HPAT.Checkpointing.hpat_get_checkpoint_time(", loop_index, "))\n")
        end
        # This num_pes < 1 ? loop_index : loop_index is a hack to get the input to start_checkpoint multiple defined
        # so that ParallelIR doesn't incorrectly hoist the start of the checkpoint before the conditional.
        checkpoint_func_str = string(checkpoint_func_str, "        checkpoint_handle = HPAT.Checkpointing.hpat_start_checkpoint(num_pes < 1 ? ", loop_index, ":", loop_index, ")\n")
        for i = 1:length(argument_names)
        checkpoint_func_str = string(checkpoint_func_str, "        HPAT.Checkpointing.hpat_value_checkpoint(checkpoint_handle, ", argument_names[i], ")\n")
        end
        checkpoint_func_str = string(checkpoint_func_str, "        HPAT.Checkpointing.hpat_end_checkpoint(checkpoint_handle)\n")
        checkpoint_func_str = string(checkpoint_func_str, "        return cur_time\n")
        checkpoint_func_str = string(checkpoint_func_str, "    end\n")
        checkpoint_func_str = string(checkpoint_func_str, "    return start_time\n")
        checkpoint_func_str = string(checkpoint_func_str, "end\n")
        @dprintln(3,"checkpoint_func_str = \n", checkpoint_func_str)
        Main.eval(parse(checkpoint_func_str))   # Force the new checkpoint function into existence.

        assert(!isempty(loop_entry_bb.statements))
        first_loop_entry_stmt = loop_entry_bb.statements[1]
        # We can use this version with the __hpat_num_pes variable once we are able to run the checkpointing pass after the distributed pass.
        #call_checkpoint_expr  = ParallelAccelerator.ParallelIR.mk_assignment_expr(checkpoint_timer_sn, ParallelAccelerator.ParallelIR.TypedExpr(Uint64, :call, TopNode(symbol(checkpoint_func_name)), checkpoint_timer_sn, :__hpat_num_pes, live_in_and_def...))
        pes_expr = ParallelAccelerator.ParallelIR.TypedExpr(Int32, :call, TopNode(:hpat_dist_num_pes))
        #pes_expr = 8
        call_checkpoint_expr  = ParallelAccelerator.ParallelIR.mk_assignment_expr(checkpoint_timer_sn, ParallelAccelerator.ParallelIR.TypedExpr(Float64, :call, GlobalRef(Main,symbol(checkpoint_func_name)), checkpoint_timer_sn, pes_expr, live_in_and_def...))
        CompilerTools.CFGs.insertStatementBefore(lives.cfg, loop_entry_bb, first_loop_entry_stmt.index, call_checkpoint_expr)
    end

    # Reconstitute the body with the new basic block.
    body.args = CompilerTools.CFGs.createFunctionBody(lives.cfg)
    # Re-create the lambda.
    new_ast = CompilerTools.LambdaHandling.LambdaVarInfoToLambdaExpr(state.LambdaVarInfo, body)

    @dprintln(1,"Checkpointing.from_root returns function = ", function_name, " ast = ", new_ast)
    return new_ast
end

end # module


text:
function mm_spec_dev1_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(4, 2, 3, 2, 3, 2, 1)
    return dimspec
end

function mm_spec_dev1_defn(x, pr, ps)

    f = exp(-x[1]+x[2]+2x[3]-5x[4]) + x[1]^2*x[4]^3*pr[1]^2 + x[2]*pr[2]^2 + x[3]*pr[1]*pr[2]*ps[1]

    ce = Vector{Num}(undef, 2)
    ce[1] = (x[1]-x[2])^2 - pr[1]*x[2]*x[3]^3
    ce[2] = (x[4]-x[1])^3 - pr[1]^2*x[4]

    ci = Vector{Num}(undef, 3)
    ci[1] = sin(0.5*pi*x[1])*pr[2]^3 + cos(x[1]^2*x[2]^3*x[3]^3*x[4]^4)
    ci[2] = pr[1]*(x[3]-2x[2])^4
    ci[3] = -x[4]^3*pr[2]^2 + x[1]^2

    F = Matrix{Num}(undef, 3, 2)
    F[1,1] = x[1]^3*x[2]^3*pr[1]^2
    F[2,1] = x[2]
    F[3,1] = x[3]
    F[1,2] = (x[1]+x[4])^5
    F[2,2] = pr[1]^2*x[3]^2
    F[3,2] = pr[1]*pr[2]*(x[1]-x[2])^2

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_dev1_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_dev1_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_dev1_knownsols()

    return [ ]
end

function mm_spec_dev1_parameterisations(t)

    return Vector{MPCCParameterisationDefn}( [ ] )
end


text:
"""
    LookingGlass.jl

This package contains a collection of code reflection and investigation utilities that can
be (potentially) useful for understanding, analyzing, and investigating julia code.

Most of the reflection functions are named hierarchically by what they reflect over. For
example:
- `func_specializations(f)` lists all specializations for a given function.
- `module_functions(m::Module)` lists all the functions in a given module.
- `module_globals_names(m::Module)` lists the names of all global variables in a given module.

There are also other non-reflection-focused utilities, such as `@quot`.
"""
module LookingGlass
import Base.Iterators
using InteractiveUtils: subtypes

export supertypes, print_typetree

"""
    @quot 2 + 3

Macro for debugging during metaprogramming; it simply returns the exact expression that a
macro gets as input. This can be helpful if you're wondering what inputs a macro can expect.
"""
macro quot(e) QuoteNode(e) end

# ---------------------------------------------------------------------------
# -- Utilities for Functions
# ---------------------------------------------------------------------------

"""
    func_specializations(f[, tt]) -> Dict(MethodInstance => Method)

Return all the specializations (MethodInstances) of each method of the given function (or
callable object), if there are any.

Specializations (also called MethodInstances) are the compiler-generated code for a method,
generated by specializing the method for specific argument types.

```julia-repl
julia> foo(x) = 2x
foo (generic function with 1 method)

julia> foo(2) + foo(2.0)
8.0

julia> keys(LookingGlass.func_specializations(foo))
Base.KeySet for a Dict{Core.MethodInstance,Method} with 2 entries. Keys:
  MethodInstance for foo(::Float64)
  MethodInstance for foo(::Int64)
```
"""
function func_specializations end

func_specializations(f) =
    Dict(mi => m
        for m in methods(f).ms
        for mi in _method_specializations(m)
    )
func_specializations(f, tt) =
    Dict(mi => m
        for m in methods(f, tt).ms
        for mi in _method_specializations(m)
    )

method_specializations(m) = collect(_method_specializations(m))

function _method_specializations(m::Method)
    s = m.specializations
    if VERSION > v"1.5-"
        return (s[i] for i in 1:length(s) if isassigned(s, i))
    else
        # m.specializations returns a custom iterate, which we collect here via `s.next`
        return Channel() do ch
            while s isa Core.TypeMapEntry  # TODO: Can there be other things in here besides TypeMapEntry?
                put!(ch, s.func)
                s=s.next
            end
        end
    end
end

"""
    func_backedges(f) -> Dict((method, signature) => [backedges...])

Return all the backedges on all specializations of a function, if there are any.
Backedges are a compiler-internal mechanism used for recursively invalidating the cached
compilations of a method -- its MethodInstances -- if the inlined code for any functions it
calls changes. This is the mechanism that fixes the infamous Julia #265 Issue.

Since these are _back_ edges, not _forward_ edges, if a function `foo` has a backedge listed
to function `bar`, this means `bar` has _inlined `foo`_ when it was compiled, so if the code
for `foo` is redefined by the user, `bar` needs to be recompiled the next time it's called.

```julia-repl
julia> foo(x) = 2x
foo (generic function with 1 method)

julia> bar(x) = foo(x) + 1
bar (generic function with 1 method)

julia> bar(1)
3

julia> LookingGlass.func_backedges(foo)
Dict{Any,Array{Any,1}} with 2 entries:
  (foo(x) in Main at none:1, Tuple{typeof(foo),Int64}) => Any[MethodInstance for bar(::Int64)]
  :MethodTable                                         => Any[]
```

Backedges on a function's MethodTable means changes to _any methods_ (including adding
a method) will trigger recompilation of the target function.
"""
function func_backedges(f)
    out = Dict{Any, Vector{Any}}(
        # Method backedges
        mi => try mi.backedges catch ; [] end
        for (mi,m) in func_specializations(f)
    )
    # MethodTable backedges
    for (typ, method) in
            # MethodTable edges are pairs of (type, method_instance) (but stored as a flat array)
            try Tuple.(Iterators.partition(methods(f).mt.backedges, 2)) catch ; [] end
        methods = get!(out, (:MethodTable, typ), [])
        push!(methods, method)
    end
    out
end

# ---------------------------------------------------------------------------
# -- Utilities for Modules
# ---------------------------------------------------------------------------

"""
    module_submodules(m::Module; recursive=true, base=false) -> Vector{Module}
Return a list of all Modules that are submodules of `m`. If `recursive=true`, this returns
all recursive submodules. If `recursive=false`, it returns only direct children of `m`.

By default, it skips Julia's `Base` and `Core` modules, but you can enable those with
`base=true`.
NOTE: currently, `module_submodules(Main, recursive=true, base=true)` will trigger a
StackOverflowError, since recursively listing submodules of `Core` or `Base` infinite loops.
"""
module_submodules(m::Module; recursive=true, base=false, kwargs...) =
    if recursive
        _module_recursive_submodules(m; base=base, seen=Set{Module}(), kwargs...)
    else
        _module_direct_submodules(m; base=base, kwargs...)
    end

function _module_recursive_submodules(m; base, seen, kwargs...)
    submodules = _module_direct_submodules(m; base=base, seen=seen, kwargs...)
    for m in submodules
        push!(seen, m)
    end
    modules = collect(Iterators.flatten(
        [x, _module_recursive_submodules(x; base=base, seen=seen, kwargs...)...]
            for x in submodules
        ))
    return modules
end
_module_direct_submodules(m; base, seen, kwargs...) =
    Module[submodule
        for x in filter(x->name_is_submodule(m,x), names(m; all=true, kwargs...))
        for submodule in (Core.eval(m, x),)  # assign to temporary variable (comprehensions are weird)
        if submodule ∉ seen]

name_is_submodule(m::Module, s::Symbol) =
    isdefined(m, s) && isa(Core.eval(m,s), Module) && nameof(m) != s
#module_is_submodule(m::Module, s::Module) = issubmodule(m, nameof(s))  # BROKEN: might just be a name collision. This is not the right way to check this.

"""
    module_objects(m::Module)
Return a list of the names of _all reachable objects_ in Module `m`.

There are more-specific helper functions as well:
 - `module_submodules(m, n)`
 - `module_functions_names(m, n)`
 - `module_globals(m, n)`
 - ...
"""
module_objects(m::Module) = [
    Core.eval(m, n)
    for n in names(m, all=true)
    if isdefined(m, n)]

"""
    module_functions_names(m::Module) -> Vector{Symbol}
Return a list of the names of all functions defined in Module `m`.
"""
module_functions_names(m::Module) =
    [name
     for name in names(m, all=true)
     if module_name_isfunction(m, name)]

function module_name_isfunction(m::Module, name::Symbol)
    # TODO: Why is this isdefined() check needed? Why can `names(all=true)` return names
    # that aren't defined?
    return isdefined(m, name) &&
        name ∉ (:include, :eval) && Core.eval(m, name) isa Function
end

"""
    module_functions(m::Module) -> Dict{Symbol, Function}
Return a list of all the Function objects defined in Module `m`.
"""
module_functions(m::Module) = [Core.eval(m, n) for n in module_functions_names(m)]

"""
    module_globals_names(m; constness=:all, mutability=:all) -> Vector{Symbol}
Return a list of the names of all global variables in Module `m`.

To return only const globals or only non-const globals, set the `constness=` keyword
argument to one of `constness=:const` or `constness=:nonconst`.
To return only mutable globals or only immutable globals, set the `mutability=` keyword
argument to one of `mutability=:mutable` or `mutability=:immutable`.
"""
module_globals_names(m::Module; constness=:all, mutability=:all, kwargs...) =
    [n for n in names(m, all=true; kwargs...)
        if module_name_isglobal(m, n; constness=constness, mutability=mutability)]
function module_name_isglobal(m::Module, n::Symbol; constness, mutability)
    @assert constness ∈ (:all, :const, :nonconst)
    @assert mutability ∈ (:all, :mutable, :immutable)
    try
        if String(n)[1] == '#' return false end
        v = Core.eval(m, n)
        return !isa(v, Union{DataType, UnionAll, Function, Module}) &&
            (constness == :all || (constness == :const && _isconst_global(m, n) ||
                                   constness == :nonconst && !_isconst_global(m, n))) &&
            (mutability == :all || (mutability == :mutable && !isimmutable(v) ||
                                    mutability == :immutable && isimmutable(v)))
    catch
        false
    end
end
"""
    module_globals(m::Module) -> Dict{Symbol, Any}
Return a Dict mapping the name to the value of all global variables in Module `m`.
"""
module_globals(m::Module) = Dict(n => Core.eval(m, n) for n in module_globals_names(m))

"""
    module_recursive_globals(m::Module) -> Dict{Tuple{Module,Symbol}, Any}
Return a Dict mapping the fully qualified name to the value of all global variables in
Module `m` and all its recursive submodules.
"""
function module_recursive_globals(m::Module; kwargs...)
    return Dict((mod,n) => Core.eval(mod, n)
                for (mod,names) in module_recursive_globals_names(m; kwargs...)
                for n in names
                )
end

function _isconst_global(m::Module, s::Symbol)
    b = _getbinding(m,s)
    b != nothing && b.constp
end

#--- Compiler-Internal Julia Binding struct ------------------------------------------------
# Stores information about a variable binding, which we can interrogate for information
# such as whether a variable is a global `const`.
# NOTE: This is very specific to a given version of Julia.
struct _Binding_t
    name::Ptr{Nothing}
    value::Ptr{Nothing}
    globalref::Ptr{Nothing}
    owner::Ptr{Nothing}
    constp::Bool
end
function _getbinding(m, s)
    p = ccall(:jl_get_binding, Ptr{_Binding_t}, (Any, Any), m, s)
    p == Ptr{Nothing}(0) ? nothing : unsafe_load(p)
end
#-------------------------------------------------------------------------------------------

"""
    module_recursive_globals_names(m; constness=:all, mutability=:all) -> Dict(submodule => name)
Return a list of the names of all global variables for each submodule in Module `m`.

Defaults to all globals, can toggle only const-globals or nonconst-globals via keyword args.
See Also:
- [`module_globals`](@ref)
- [`module_recursive_globals`](@ref)
"""
module_recursive_globals_names(m::Module; constness=:all, mutability=:all, kwargs...) =
    merge!(
        Dict(m => module_globals_names(m; constness=constness, mutability=mutability, kwargs...)),
        Dict(
            sm => names
            for sm in module_submodules(m; recursive=true, kwargs...)
            for names in (module_globals_names(sm; constness=constness, mutability=mutability, kwargs...),)
            #if !isempty(names)
        ))

"""
    print_typetree(t::Type)

Prints the type hierarchy rooted at the specified type. Note that this is intended to be used interactively, so the current implementation works by printing to `stdout`, and returns `nothing`.

```julia-repl
julia> print_typetree(Integer)
Integer
    Bool
    Signed
        BigInt
        Int128
        Int16
        Int32
        Int64
        Int8
    Unsigned
        UInt128
        UInt16
        UInt32
        UInt64
        UInt8
```
"""
function print_typetree(t::Type, depth=0)
    println(repeat("    ", depth), t)

    for s in subtypes(t)
        print_typetree(s, depth + 1)
    end

    return nothing
end

"""
    supertypes(x) :: Vector{Type}

`supertypes(x::Type)` returns a `Vector{Type}` of all supertypes of `x`.

If `x` is not a `Type`, `supertypes(x::T)` returns a `Vector{Type}` containing `T` and all supertypes of `T`.

```julia-repl
julia> supertypes(2)
6-element Array{Type,1}:
 Int64
 Signed
 Integer
 Real
 Number
 Any

julia> supertypes(Int64)
5-element Array{Type,1}:
 Signed
 Integer
 Real
 Number
 Any
```
"""
function supertypes end

function supertypes(t::Type)
    t0 = t
    t1 = supertype(t)
    result = Type[]
    while t1 ≠ t0
        push!(result, t1)
        t0, t1 = t1, supertype(t1)
    end
    return result
end

supertypes(x::T) where T = [T, supertypes(T)...]

end # module


text:
module EmbeddingsAnalysis

    using LinearAlgebra
    using Statistics
    using Languages
    using Word2Vec
    using ConceptnetNumberbatch
    using StatsBase
    using MultivariateStats
    using Distances
    using QuantizedArrays

    import Base: size
    import Word2Vec: analogy_words

    export conceptnet2wv,
           CompressedWordVectors,
           compressedwordvectors,
           compress,
           analogy_words,
           write2disk,
           similarity_order,
           pca_reduction

    include("defaults.jl")          # defaults
    include("conceptnet2wv.jl")     # convert ConceptNet to WordVectors
    include("cwv.jl")               # CompressedWordVectors
    include("write2disk.jl")        # save WordVectors to disk
    include("similarity_order.jl")  # preprocess WordVectors
    include("pca_reduction.jl")     # preprocess/reduce dimensionality of WordVectors

end # module


text:
# ---
# title: 1195. Fizz Buzz Multithreaded
# id: problem1195
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/fizz-buzz-multithreaded/description/>
# hidden: true
# ---
#
# Write a program that outputs the string representation of numbers from 1 to
# _n_ , however:
#
#   * If the number is divisible by 3, output "fizz".
#   * If the number is divisible by 5, output "buzz".
#   * If the number is divisible by both 3 and 5, output "fizzbuzz".
#
# For example, for `n = 15`, we output: `1, 2, fizz, 4, buzz, fizz, 7, 8, fizz,
# buzz, 11, fizz, 13, 14, fizzbuzz`.
#
# Suppose you are given the following code:
#
#
#
#     class FizzBuzz {
#       public FizzBuzz(int n) { ... }               // constructor
#       public void fizz(printFizz) { ... }          // only output "fizz"
#       public void buzz(printBuzz) { ... }          // only output "buzz"
#       public void fizzbuzz(printFizzBuzz) { ... }  // only output "fizzbuzz"
#       public void number(printNumber) { ... }      // only output the numbers
#     }
#
# Implement a multithreaded version of `FizzBuzz` with **four** threads. The
# same instance of `FizzBuzz` will be passed to four different threads:
#
#   1. Thread A will call `fizz()` to check for divisibility of 3 and outputs `fizz`.
#   2. Thread B will call `buzz()` to check for divisibility of 5 and outputs `buzz`.
#   3. Thread C will call `fizzbuzz()` to check for divisibility of 3 and 5 and outputs `fizzbuzz`.
#   4. Thread D will call `number()` which should only output the numbers.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using BinaryDecisionDiagrams
using Test

a=BDD("a",BDD(true),BDD(false))
b=BDD("a",BDD(false),BDD(true))
c=BDD("b",b,BDD(false))
d=BDD("c",BDD(true),a)
e=BDD("c",BDD(true),a)
f=BDD("a",BDD(true),
          BDD("b",BDD("c",
                      BDD(false),
                      BDD(true)),
                  BDD(false)))

@testset "BDD tests" begin
  @test string(a) == "~a"
  @test string(b) == "a"
  @test string(c) == "(~b & a)"
  @test string(d) == "(~c | (c & ~a))"
  @test BDD(false) == BDD(false)
  @test BDD(true) == BDD(true)
  @test e == d
end

@testset "OBDD tests" begin
  ordering=["c","b","a"]

  oa=OBDD(ordering,c)
  ob=OBDD(ordering,d)
  oc=OBDD("(c,b,a)-> (~(c|true) || !1) | (true&(b|a))")

  otrue=OBDD(oa.ordering,BDD(true))
  ofalse=OBDD(oa.ordering,BDD(false))

  BDDnodes=(()->Set(["$(node)" for node in union(ancestors(BDD(0)),
                                                 ancestors(BDD(1)))]))

  GC.gc()
  before_od=BDDnodes()

  od=OBDD(["d","c","b","a"],BDD("d",c,d))
  with_od=BDDnodes()

  od=nothing
  GC.gc()
  after_od=BDDnodes()

  @test_broken before_od == after_od
  @test before_od != with_od
  @test issubset(before_od,with_od)

  @testset "Parsing" begin
    @test oc == OBDD(ordering,"b|a")
    @test OBDD(string(oa)) == oa
    @test OBDD(string(oc)) == oc
    @test OBDD(string(~oa)) == ~oa
    @test OBDD(string(~ob)) == ~ob
    @test OBDD(string(~oc)) == ~oc
    @test OBDD(string(oa|ob)) == oa|ob

    oe = OBDD(["c","a"], d)
    @test oe =="(~c || (c && ~a))"
    @test oe == "~c + (c*~a)"

    @test OBDD(["a","b","c"], f) == "a->(~b&c)"
    @test OBDD(string(oa|ob)) == OBDD(string(oa))|OBDD(string(ob))

    @test_throws Meta.ParseError OBDD(ordering,"(a,b)-(a&b)")
  end

  @testset "Bitwise operators" begin
    @test (~~oa).root == oa.root
    @test ~~oa == oa
    @test oa&oa == oa
    @test oa|oa == oa
    @test ~oa&oa == OBDD(oa.ordering,BDD(false))
    @test ~oa|oa == OBDD(oa.ordering,BDD(true))
    @test oa&BDD(false) == BDD(false)&oa
    @test oa&false == false&oa
    @test oa&0 == 0&oa
    @test oa&BDD(false) == ofalse
    @test oa&BDD(true) == BDD(true)&oa
    @test BDD(true)&oa == oa&BDD(true)
    @test true&oa == oa&true
    @test 1&oa == oa&1
    @test oa&true == oa
    @test oa|true == otrue
    @test oa|false == oa
    @test oa|ob == ~(~oa & ~ob)
    @test oa|ob == ob|oa
    @test oa&ob == ~(~oa | ~ob)
    @test oa&ob == ob&oa
    @test oa⊻ob == ob⊻oa
    @test otrue⊻otrue == ofalse
    @test ofalse⊻ofalse == ofalse
    @test otrue⊻ofalse == otrue
    @test ofalse⊻otrue == otrue
  end

  ab = SAT_assignments(ob)
  aa = SAT_assignments(oa)
  @testset "SAT" begin
    @test Set(ab) == Set([Dict("c" => 0), Dict("a" => 0)])
    @test Set(aa) == Set([Dict("b" =>0,"a" => 1)])
  end

  @testset "Restrictions" begin
    @test restrict(oc,"b",1) == otrue
    @test restrict(oc,"b",0) == OBDD(oc.ordering,"a")

    for a in ab
      @test restrict(ob, a) == otrue
    end

    @test restrict(ob, aa[1]) == OBDD("(c,b,a)-> ~c")
  end

  binaryexp="a|~c&b"
  O1=['c','a','b']
  O2=['a','b','c']

  oe=OBDD(O1,binaryexp)
  of=OBDD(O2,binaryexp)

  @testset "Variable orderings" begin
    @test_throws ArgumentError oe==of
    @test oe==changeordering(of,O1)

    set_dynamic_ordering()

    @test oe==of
    @test OBDD(["x"],"x")|(~oe&of)==(of&~oe)|OBDD(["x"],"x")

    set_static_ordering()

    @test_throws ArgumentError oe==of
  end
end


text:
# Gaussian
using Distributions
using LinearAlgebra: Diagonal, ⋅, checksquare, cholesky, logdet
using Statistics

"""
    Gaussian(μ, Σ) -> P

Gaussian distribution with mean `μ` and covariance `Σ`. Defines `rand(P)` and `(log-)pdf(P, x)`. The type is designed to work with `StaticArrays`. Internally, a representation of the matrix square root `Σsqrt` of `Σ` is stored. The square root has to satisfy that that `Σ ≈ Σsqrt'*Σsqrt` which means that a Cholesky factor can be used but it would also be possible to store e.g. a `Eigen` factorization. The default is to compute the Choleksy factor when a matrix is passed to the constructor.
"""
struct Gaussian{T,S}
    μ::T
    Σsqrt::S
end

function Gaussian(μ::AbstractVector, Σ::AbstractMatrix)
    Σsqrt = cholesky(Σ).U
    Gaussian{typeof(μ),typeof(Σsqrt)}(μ, Σsqrt)
end
Gaussian(Σ::SMatrix{M,N}) where {M,N} = Gaussian(zero(similar_type(Σ, Size(M))), Σ)
Gaussian(Σ::Diagonal{<:Any,<:SVector}) = Gaussian(zero(similar_type(parent(Σ))), Σ)

Base.length(P::Gaussian) = length(P.μ)

Base.rand(P::Gaussian) = P.μ + P.Σsqrt'*randn(typeof(P.μ))
Base.rand(P::Gaussian, n::Integer) = [rand(P)' for i in 1:n]

function Distributions.logpdf(P::Gaussian, x)
    ε = P.Σsqrt'\(x - P.μ)
    lp  = -logdet(P.Σsqrt) - (ε⋅ε)/2
    lp -= log(2*typeof(lp)(π))*length(P)/2
    return lp
end
Distributions.pdf(P::Gaussian, x) = exp(logpdf(P::Gaussian, x))

Statistics.mean(P::Gaussian) = P.μ
Statistics.cov(P::Gaussian) = P.Σsqrt'*P.Σsqrt

max_stars_repo_name: andreasnoack/MultivariateDistributions

max_issues_repo_name: andreasnoack/MultivariateDistributions

max_forks_repo_name: andreasnoack/MultivariateDistributions


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Statistics: mean

struct Policy{S, A, V}
    shared::S
    action::A
    value::V
end
(m::Policy)(s) = (l = m.shared(s); (softmax(m.action(l)), m.value(l)))
Flux.@functor Policy

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

    for i in 1:epochs
        episode!(episode, env, maxt) do s
            p, v = policy(s)
            return p, v[1]
        end
        nstepr!(episode) #discount!(episode)
        adv!(episode, true)

        Flux.train!(params(policy), [batch(episode)], opt) do s, a, nr, adv
            p, v = policy(s)
            return -mean(log.(gather(p, a)) .* adv) + 0.5f0 * Flux.mse(vec(v), nr)
        end
        i % 5 == 0 && (rewards[div(i, 5)] = test(s -> Flux.onecold(policy(s)[1]), env))
    end
    return rewards
end

env = CartPole{Float32}()

policy = Policy(
    Dense(length(env.observationspace), 32, relu),
    Dense(32, length(env.actionspace)),
    Dense(32, 1)
)

r = run!(policy, ADAM(0.005), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "a2c")


text:
using IndexFunArrays, ROIViews

# a random position within the size limit sz
randpos(sz, dims) = Tuple((d <= length(sz)) ? rand(1:sz[d]) : 1 for d in 1:dims)

@testset "ROIView" begin
    for d in 1:5
        sz = Tuple(9 .+ ones(Int,d));
        data = idx(sz,offset=CtrCorner);
        N = 3
        for s = 1:d # smaller ROI dimensions. s is ROI dimensionality
            ROIpos = Tuple(Tuple(rand(1:10, s)) for n in 1:N)
            ROIsize = Tuple(rand(1:10,s))
            pad_val = 0 .* data[1];
            rois = ROIView(data, ROIpos, ROIsize, pad_val=pad_val);
            for p = 1:N # select a ROI number
                p_roi = randpos(ROIsize,d) # position in the ROI data
                @test size(rois) == (expand_size(ROIsize,size(data))...,N)
                roi_val = rois[p_roi...,p]
                pos = p_roi .+ Tuple((n<=s) ? ROIpos[p][n] .- ROIsize[n].÷2 : 0 for n=1:d)
                if Base.checkbounds(Bool, data, pos...)
                    data_val = data[pos...]
                    @test roi_val .+ 1 == pos
                else
                    data_val = pad_val
                end
                @test data_val == roi_val
            end
        end
    end
end


text:
struct RandomGuessing <: Algorithm end

function train(learner::Learner{RandomGuessing})
    num_samples = sample_count(learner.train_data)
    num_features = feature_count(learner.train_data)
    β_values = _unif_random(learner.param_lb, learner.param_ub)
    offset = rand(Uniform(learner.offset_lb, learner.offset_ub))
    _audit_weights(learner, β_values)
    return TrainingResults(learner, LinearModel(β_values, offset))
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
using Random:
    AbstractRNG, MersenneTwister, randperm, seed!, shuffle!
using Statistics: mean

using LightGraphs:
    getRNG, sample!

"""
    SimpleGraph{T}(nv, ne; seed=-1)

Construct a random `SimpleGraph{T}` with `nv` vertices and `ne` edges.
The graph is sampled uniformly from all such graphs.
If `seed >= 0`, a random generator is seeded with this value.
If not specified, the element type `T` is the type of `nv`.

### See also
[`erdos_renyi`](@ref)

## Examples
```jldoctest
julia> SimpleGraph(5, 7)
{5, 7} undirected simple Int64 graph
```
"""
function SimpleGraph{T}(nv::Integer, ne::Integer; seed::Int=-1) where T <: Integer
    tnv = T(nv)
    maxe = div(Int(nv) * (nv - 1), 2)
    @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe")
    ne > (2 / 3) * maxe && return complement(SimpleGraph(nv, maxe - ne))

    rng = getRNG(seed)
    g = SimpleGraph(tnv)

    while g.ne < ne
        source = rand(rng, one(T):tnv)
        dest = rand(rng, one(T):tnv)
        source != dest && add_edge!(g, source, dest)
    end
    return g
end

SimpleGraph(nv::T, ne::Integer; seed::Int=-1) where T <: Integer =
    SimpleGraph{T}(nv, ne, seed=seed)

"""
    SimpleDiGraph{T}(nv, ne; seed=-1)

Construct a random `SimpleDiGraph{T}` with `nv` vertices and `ne` edges.
The graph is sampled uniformly from all such graphs.
If `seed >= 0`, a random generator is seeded with this value.
If not specified, the element type `T` is the type of `nv`.

### See also
[`erdos_renyi`](@ref)

## Examples
```jldoctest
julia> SimpleDiGraph(5, 7)
{5, 7} directed simple Int64 graph
```
"""
function SimpleDiGraph{T}(nv::Integer, ne::Integer; seed::Int=-1) where T <: Integer
    tnv = T(nv)
    maxe = Int(nv) * (nv - 1)
    @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe")
    ne > (2 / 3) * maxe && return complement(SimpleDiGraph{T}(nv, maxe - ne))

    rng = getRNG(seed)
    g = SimpleDiGraph(tnv)
    while g.ne < ne
        source = rand(rng, one(T):tnv)
        dest = rand(rng, one(T):tnv)
        source != dest && add_edge!(g, source, dest)
    end
    return g
end

SimpleDiGraph(nv::T, ne::Integer; seed::Int=-1) where T <: Integer =
    SimpleDiGraph{Int}(nv, ne, seed=seed)

"""
    randbn(n, p, seed=-1)

Return a binomally-distribted random number with parameters `n` and `p` and optional `seed`.

### References
- "Non-Uniform Random Variate Generation," Luc Devroye, p. 522. Retrieved via http://www.eirene.de/Devroye.pdf.
- http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java
"""
function randbn(n::Integer, p::Real, seed::Integer=-1)
    rng = getRNG(seed)
    log_q = log(1.0 - p)
    x = 0
    sum = 0.0
    while true
        sum += log(rand(rng)) / (n - x)
        sum < log_q && break
        x += 1
    end
    return x
end

"""
    erdos_renyi(n, p)

Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model)
random graph with `n` vertices. Edges are added between pairs of vertices with
probability `p`.

### Optional Arguments
- `is_directed=false`: if true, return a directed graph.
- `seed=-1`: set the RNG seed.
"""
function erdos_renyi(n::Integer, p::Real; is_directed=false, seed::Integer=-1)
    p >= 1 && return is_directed ? CompleteDiGraph(n) : CompleteGraph(n)
    m = is_directed ? n * (n - 1) : div(n * (n - 1), 2)
    ne = randbn(m, p, seed)
    return is_directed ? SimpleDiGraph(n, ne, seed=seed) : SimpleGraph(n, ne, seed=seed)
end

"""
    erdos_renyi(n, ne)

Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model) random
graph with `n` vertices and `ne` edges.

### Optional Arguments
- `is_directed=false`: if true, return a directed graph.
- `seed=-1`: set the RNG seed.
"""
function erdos_renyi(n::Integer, ne::Integer; is_directed=false, seed::Integer=-1)
    return is_directed ? SimpleDiGraph(n, ne, seed=seed) : SimpleGraph(n, ne, seed=seed)
end

"""
    expected_degree_graph(ω)

Given a vector of expected degrees `ω` indexed by vertex, create a random undirected graph in which vertices `i` and `j` are
connected with probability `ω[i]*ω[j]/sum(ω)`.

### Optional Arguments
- `seed=-1`: set the RNG seed.

### Implementation Notes
The algorithm should work well for `maximum(ω) << sum(ω)`. As `maximum(ω)` approaches `sum(ω)`, some deviations
from the expected values are likely.

### References
- Connected Components in Random Graphs with Given Expected Degree Sequences, Linyuan Lu and Fan Chung. [https://link.springer.com/article/10.1007%2FPL00012580](https://link.springer.com/article/10.1007%2FPL00012580)
- Efficient Generation of Networks with Given Expected Degrees, Joel C. Miller and Aric Hagberg. [https://doi.org/10.1007/978-3-642-21286-4_10](https://doi.org/10.1007/978-3-642-21286-4_10)
"""
function expected_degree_graph(ω::Vector{T}; seed::Int=-1) where T <: Real
    g = SimpleGraph(length(ω))
    expected_degree_graph!(g, ω, seed=seed)
end

function expected_degree_graph!(g::SimpleGraph, ω::Vector{T}; seed::Int=-1) where T <: Real
    n = length(ω)
    @assert all(zero(T) .<= ω .<= n - one(T)) "Elements of ω needs to be at least 0 and at most n-1"

    π = sortperm(ω, rev=true)
    rng = getRNG(seed)

    S = sum(ω)

    for u = 1:(n - 1)
        v = u + 1
        p = min(ω[π[u]] * ω[π[v]] / S, one(T))
        while v <= n && p > zero(p)
            if p != one(T)
                v += floor(Int, log(rand(rng)) / log(one(T) - p))
            end
            if v <= n
                q = min(ω[π[u]] * ω[π[v]] / S, one(T))
                if rand(rng) < q / p
                    add_edge!(g, π[u], π[v])
                end
                p = q
                v += 1
            end
        end
    end
    return g
end

"""
    watts_strogatz(n, k, β)

Return a [Watts-Strogatz](https://en.wikipedia.org/wiki/Watts_and_Strogatz_model)
small model random graph with `n` vertices, each with degree `k`. Edges are
randomized per the model based on probability `β`.

### Optional Arguments
- `is_directed=false`: if true, return a directed graph.
- `seed=-1`: set the RNG seed.
"""
function watts_strogatz(n::Integer, k::Integer, β::Real; is_directed=false, seed::Int=-1)
    @assert k < n / 2
    if is_directed
        g = SimpleDiGraph(n)
    else
        g = SimpleGraph(n)
    end
    rng = getRNG(seed)
    for s in 1:n
        for i in 1:(floor(Integer, k / 2))
            target = ((s + i - 1) % n) + 1
            if rand(rng) > β && !has_edge(g, s, target) # TODO: optimize this based on return of add_edge!
                add_edge!(g, s, target)
            else
                while true
                    d = target
                    while d == target
                        d = rand(rng, 1:(n - 1))
                        if s < d
                            d += 1
                        end
                    end
                    if s != d
                        add_edge!(g, s, d) && break
                    end
                end
            end
        end
    end
    return g
end

function _suitable(edges::Set{SimpleEdge{T}}, potential_edges::Dict{T,T}) where T <: Integer
    isempty(potential_edges) && return true
    list = keys(potential_edges)
    for s1 in list, s2 in list
        s1 >= s2 && continue
        (SimpleEdge(s1, s2) ∉ edges) && return true
    end
    return false
end

_try_creation(n::Integer, k::Integer, rng::AbstractRNG) = _try_creation(n, fill(k, n), rng)

function _try_creation(n::T, k::Vector{T}, rng::AbstractRNG) where T <: Integer
    edges = Set{SimpleEdge{T}}()
    m = 0
    stubs = zeros(T, sum(k))
    for i = one(T):n
        for j = one(T):k[i]
            m += 1
            stubs[m] = i
        end
    end
    # stubs = vcat([fill(i, k[i]) for i = 1:n]...) # slower

    while !isempty(stubs)
        potential_edges =  Dict{T,T}()
        shuffle!(rng, stubs)
        for i in 1:2:length(stubs)
            s1, s2 = stubs[i:(i + 1)]
            if (s1 > s2)
                s1, s2 = s2, s1
            end
            e = SimpleEdge(s1, s2)
            if s1 != s2 && ∉(e, edges)
                push!(edges, e)
            else
                potential_edges[s1] = get(potential_edges, s1, 0) + 1
                potential_edges[s2] = get(potential_edges, s2, 0) + 1
            end
        end

        if !_suitable(edges, potential_edges)
            return Set{SimpleEdge{T}}()
        end

        stubs = Vector{Int}()
        for (e, ct) in potential_edges
            append!(stubs, fill(e, ct))
        end
    end
    return edges
end

"""
    barabasi_albert(n, k)

Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model)
random graph with `n` vertices. It is grown by adding new vertices to an initial
graph with `k` vertices. Each new vertex is attached with `k` edges to `k`
different vertices already present in the system by preferential attachment.
Initial graphs are undirected and consist of isolated vertices by default.

### Optional Arguments
- `is_directed=false`: if true, return a directed graph.
- `complete=false`: if true, use a complete graph for the initial graph.
- `seed=-1`: set the RNG seed.
"""
barabasi_albert(n::Integer, k::Integer; keyargs...) =
barabasi_albert(n, k, k; keyargs...)

"""
    barabasi_albert(n::Integer, n0::Integer, k::Integer)

Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model)
random graph with `n` vertices. It is grown by adding new vertices to an initial
graph with `n0` vertices. Each new vertex is attached with `k` edges to `k`
different vertices already present in the system by preferential attachment.
Initial graphs are undirected and consist of isolated vertices by default.

### Optional Arguments
- `is_directed=false`: if true, return a directed graph.
- `complete=false`: if true, use a complete graph for the initial graph.
- `seed=-1`: set the RNG seed.
"""
function barabasi_albert(n::Integer, n0::Integer, k::Integer; is_directed::Bool=false, complete::Bool=false, seed::Int=-1)
    if complete
        g = is_directed ? CompleteDiGraph(n0) : CompleteGraph(n0)
    else
        g = is_directed ? SimpleDiGraph(n0) : SimpleGraph(n0)
    end

    barabasi_albert!(g, n, k; seed=seed)
    return g
end

"""
    barabasi_albert!(g::AbstractGraph, n::Integer, k::Integer)

Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model)
random graph with `n` vertices. It is grown by adding new vertices to an initial
graph `g`. Each new vertex is attached with `k` edges to `k` different vertices
already present in the system by preferential attachment.

### Optional Arguments
- `seed=-1`: set the RNG seed.
"""
function barabasi_albert!(g::AbstractGraph, n::Integer, k::Integer; seed::Int=-1)
    n0 = nv(g)
    1 <= k <= n0 <= n ||
        throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n"))
    n0 == n && return g

    # seed random number generator
    seed > 0 && seed!(seed)

    # add missing vertices
    sizehint!(g.fadjlist, n)
    add_vertices!(g, n - n0)

    # if initial graph doesn't contain any edges
    # expand it by one vertex and add k edges from this additional node
    if ne(g) == 0
        # expand initial graph
        n0 += one(n0)

        # add edges to k existing vertices
        for target in sample!(collect(1:(n0 - 1)), k)
            add_edge!(g, n0, target)
        end
    end

    # vector of weighted vertices (each node is repeated once for each adjacent edge)
    weightedVs = Vector{Int}(undef, 2 * (n - n0) * k + 2 * ne(g))

    # initialize vector of weighted vertices
    offset = 0
    for e in edges(g)
        weightedVs[offset += 1] = src(e)
        weightedVs[offset += 1] = dst(e)
    end

    # array to record if a node is picked
    picked = fill(false, n)

    # vector of targets
    targets = Vector{Int}(undef, k)

    for source in (n0 + 1):n
        # choose k targets from the existing vertices
        # pick uniformly from weightedVs (preferential attachement)
        i = 0
        while i < k
            target = weightedVs[rand(1:offset)]
            if !picked[target]
                targets[i += 1] = target
                picked[target] = true
            end
        end

        # add edges to k targets
        for target in targets
            add_edge!(g, source, target)

            weightedVs[offset += 1] = source
            weightedVs[offset += 1] = target
            picked[target] = false
        end
    end
    return g
end

"""
    static_fitness_model(m, fitness)

Generate a random graph with ``|fitness|`` vertices and `m` edges,
in which the probability of the existence of ``Edge_{ij}`` is proportional
to ``fitness_i × fitness_j``.

### Optional Arguments
- `seed=-1`: set the RNG seed.

### Performance
Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``.

### References
- Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution in scale-free networks. Phys Rev Lett 87(27):278701, 2001.
"""
function static_fitness_model(m::Integer, fitness::Vector{T}; seed::Int=-1) where T <: Real
    m < 0 && throw(ArgumentError("number of edges must be positive"))
    n = length(fitness)
    m == 0 && return SimpleGraph(n)
    nvs = 0
    for f in fitness
        # sanity check for the fitness
        f < zero(T) && throw(ArgumentError("fitness scores must be non-negative"))
        f > zero(T) && (nvs += 1)
    end
    # avoid getting into an infinite loop when too many edges are requested
    max_no_of_edges = div(nvs * (nvs - 1), 2)
    m > max_no_of_edges && throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)"))
    # calculate the cumulative fitness scores
    cum_fitness = cumsum(fitness)
    g = SimpleGraph(n)
    _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, seed)
    return g
end

"""
    static_fitness_model(m, fitness_out, fitness_in)

Generate a random graph with ``|fitness\\_out + fitness\\_in|`` vertices and `m` edges,
in which the probability of the existence of ``Edge_{ij}`` is proportional with
respect to ``i ∝ fitness\\_out`` and ``j ∝ fitness\\_in``.

### Optional Arguments
- `seed=-1`: set the RNG seed.

### Performance
Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``.

### References
- Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution in scale-free networks. Phys Rev Lett 87(27):278701, 2001.
"""
function static_fitness_model(m::Integer, fitness_out::Vector{T}, fitness_in::Vector{S}; seed::Int=-1) where T <: Real where S <: Real
    m < 0 && throw(ArgumentError("number of edges must be positive"))
    n = length(fitness_out)
    length(fitness_in) != n && throw(ArgumentError("fitness_in must have the same size as fitness_out"))
    m == 0 && return SimpleDiGraph(n)
    # avoid getting into an infinite loop when too many edges are requested
    noutvs = ninvs = nvs = 0
    @inbounds for i = 1:n
        # sanity check for the fitness
        (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && error("fitness scores must be non-negative") # TODO 0.7: change to DomainError?
        fitness_out[i] > zero(T) && (noutvs += 1)
        fitness_in[i] > zero(S) && (ninvs += 1)
        (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nvs += 1)
    end
    max_no_of_edges = noutvs * ninvs - nvs
    m > max_no_of_edges && throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)"))
    # calculate the cumulative fitness scores
    cum_fitness_out = cumsum(fitness_out)
    cum_fitness_in = cumsum(fitness_in)
    g = SimpleDiGraph(n)
    _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, seed)
    return g
end

function _create_static_fitness_graph!(g::AbstractGraph, m::Integer, cum_fitness_out::Vector{T}, cum_fitness_in::Vector{S}, seed::Int) where T <: Real where S <: Real
    rng = getRNG(seed)
    max_out = cum_fitness_out[end]
    max_in = cum_fitness_in[end]
    while m > 0
        source = searchsortedfirst(cum_fitness_out, rand(rng) * max_out)
        target = searchsortedfirst(cum_fitness_in, rand(rng) * max_in)
        # skip if loop edge
        (source == target) && continue
        # is there already an edge? If so, try again
        add_edge!(g, source, target) || continue
        m -= one(m)
    end
end

"""
    static_scale_free(n, m, α)

Generate a random graph with `n` vertices, `m` edges and expected power-law
degree distribution with exponent `α`.

### Optional Arguments
- `seed=-1`: set the RNG seed.
- `finite_size_correction=true`: determines whether to use the finite size correction
proposed by Cho et al.

### Performance
Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``.

### References
- Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution in scale-free networks. Phys Rev Lett 87(27):278701, 2001.
- Chung F and Lu L: Connected components in a random graph with given degree sequences. Annals of Combinatorics 6, 125-145, 2002.
- Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009.
"""
function static_scale_free(n::Integer, m::Integer, α::Real; seed::Int=-1, finite_size_correction::Bool=true)
    n < 0 && throw(ArgumentError("number of vertices must be positive"))
    α < 2 && throw(ArgumentError("out-degree exponent must be >= 2"))
    fitness = _construct_fitness(n, α, finite_size_correction)
    static_fitness_model(m, fitness, seed=seed)
end

"""
    static_scale_free(n, m, α_out, α_in)

Generate a random graph with `n` vertices, `m` edges and expected power-law
degree distribution with exponent `α_out` for outbound edges and `α_in` for
inbound edges.

### Optional Arguments
- `seed=-1`: set the RNG seed.
- `finite_size_correction=true`: determines whether to use the finite size correction
proposed by Cho et al.

### Performance
Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``.

### References
- Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution in scale-free networks. Phys Rev Lett 87(27):278701, 2001.
- Chung F and Lu L: Connected components in a random graph with given degree sequences. Annals of Combinatorics 6, 125-145, 2002.
- Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009.
"""
function static_scale_free(n::Integer, m::Integer, α_out::Real, α_in::Float64; seed::Int=-1, finite_size_correction::Bool=true)
    n < 0 && throw(ArgumentError("number of vertices must be positive"))
    α_out < 2 && throw(ArgumentError("out-degree exponent must be >= 2"))
    α_in < 2 && throw(ArgumentError("out-degree exponent must be >= 2"))
    # construct the fitness
    fitness_out = _construct_fitness(n, α_out, finite_size_correction)
    fitness_in = _construct_fitness(n, α_in, finite_size_correction)
    # eliminate correlation
    shuffle!(fitness_in)
    static_fitness_model(m, fitness_out, fitness_in, seed=seed)
end

function _construct_fitness(n::Integer, α::Real, finite_size_correction::Bool)
    α = -1 / (α - 1)
    fitness = zeros(n)
    j = float(n)
    if finite_size_correction && α < -0.5
        # See the Cho et al paper, first page first column + footnote 7
        j += n^(1 + 1 / 2 * α) * (10 * sqrt(2) * (1 + α))^(-1 / α) - 1
    end
    j = max(j, n)
    @inbounds for i = 1:n
        fitness[i] = j^α
        j -= 1
    end
    return fitness
end

"""
    random_regular_graph(n, k)

Create a random undirected
[regular graph](https://en.wikipedia.org/wiki/Regular_graph) with `n` vertices,
each with degree `k`.

### Optional Arguments
- `seed=-1`: set the RNG seed.

### Performance
Time complexity is approximately ``\\mathcal{O}(nk^2)``.

### Implementation Notes
Allocates an array of `nk` `Int`s, and . For ``k > \\frac{n}{2}``, generates a graph of degree
``n-k-1`` and returns its complement.
"""
function random_regular_graph(n::Integer, k::Integer; seed::Int=-1)
    !iseven(n * k) && throw(ArgumentError("n * k must be even"))
    !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied"))
    if k == 0
        return SimpleGraph(n)
    end
    if (k > n / 2) && iseven(n * (n - k - 1))
        return complement(random_regular_graph(n, n - k - 1, seed=seed))
    end

    rng = getRNG(seed)

    edges = _try_creation(n, k, rng)
    while isempty(edges)
        edges = _try_creation(n, k, rng)
    end

    g = SimpleGraph(n)
    for edge in edges
        add_edge!(g, edge)
    end

    return g
end

"""
    random_configuration_model(n, ks)

Create a random undirected graph according to the [configuration model]
(http://tuvalu.santafe.edu/~aaronc/courses/5352/fall2013/csci5352_2013_L11.pdf)
containing `n` vertices, with each node `i` having degree `k[i]`.

### Optional Arguments
- `seed=-1`: set the RNG seed.
- `check_graphical=false`: if true, ensure that `k` is a graphical sequence
(see [`isgraphical`](@ref)).

### Performance
Time complexity is approximately ``\\mathcal{O}(n \\bar{k}^2)``.
### Implementation Notes
Allocates an array of ``n \\bar{k}`` `Int`s.
"""
function random_configuration_model(n::Integer, k::Array{T}; seed::Int=-1, check_graphical::Bool=false) where T <: Integer
    n != length(k) && throw(ArgumentError("a degree sequence of length n must be provided"))
    m = sum(k)
    !iseven(m) && throw(ArgumentError("sum(k) must be even"))
    !all(0 .<= k .< n) && throw(ArgumentError("the 0 <= k[i] < n inequality must be satisfied"))
    if check_graphical
        isgraphical(k) || throw(ArgumentError("degree sequence must be graphical"))
    end
    rng = getRNG(seed)

    edges = _try_creation(n, k, rng)
    while m > 0 && isempty(edges)
        edges = _try_creation(n, k, rng)
    end

    g = SimpleGraphFromIterator(edges)
    if nv(g) < n
        add_vertices!(g, n - nv(g))
    end
    return g
end

"""
    random_regular_digraph(n, k)

Create a random directed [regular graph](https://en.wikipedia.org/wiki/Regular_graph)
with `n` vertices, each with degree `k`.

### Optional Arguments
- `dir=:out`: the direction of the edges for degree parameter.
- `seed=-1`: set the RNG seed.

### Implementation Notes
Allocates an ``n × n`` sparse matrix of boolean as an adjacency matrix and
uses that to generate the directed graph.
"""
function random_regular_digraph(n::Integer, k::Integer; dir::Symbol=:out, seed::Int=-1)
    #TODO remove the function sample from StatsBase for one allowing the use
    # of a local rng
    !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied"))

    if k == 0
        return SimpleDiGraph(n)
    end
    if (k > n / 2) && iseven(n * (n - k - 1))
        return complement(random_regular_digraph(n, n - k - 1, dir=dir, seed=seed))
    end
    rng = getRNG(seed)
    cs = collect(2:n)
    i = 1
    I = Vector{Int}(undef, n * k)
    J = Vector{Int}(undef, n * k)
    V = fill(true, n * k)
    for r in 1:n
        l = ((r - 1) * k + 1):(r * k)
        I[l] .= r
        J[l] = sample!(rng, cs, k, exclude=r)
    end

    if dir == :out
        return SimpleDiGraph(sparse(I, J, V, n, n))
    else
        return SimpleDiGraph(sparse(I, J, V, n, n)')
    end
end

"""
    random_tournament_digraph(n)

Create a random directed [tournament graph]
(https://en.wikipedia.org/wiki/Tournament_%28graph_theory%29)
with `n` vertices.

### Optional Arguments
- `seed=-1`: set the RNG seed.
"""
function random_tournament_digraph(n::Integer; seed::Int=-1)

    rng = getRNG(seed)
    g = SimpleDiGraph(n)

    for i = 1:n, j = i + 1:n
        rand(rng, Bool) ? add_edge!(g, SimpleEdge(i, j)) : add_edge!(g, SimpleEdge(j, i))
    end

    return g
end

"""
    stochastic_block_model(c, n)

Return a Graph generated according to the Stochastic Block Model (SBM).

`c[a,b]` : Mean number of neighbors of a vertex in block `a` belonging to block `b`.
           Only the upper triangular part is considered, since the lower traingular is
           determined by ``c[b,a] = c[a,b] * \\frac{n[a]}{n[b]}``.
`n[a]` : Number of vertices in block `a`

### Optional Arguments
- `seed=-1`: set the RNG seed.

For a dynamic version of the SBM see the [`StochasticBlockModel`](@ref) type and
related functions.
"""
function stochastic_block_model(c::Matrix{T}, n::Vector{U}; seed::Int=-1) where T <: Real where U <: Integer
    size(c, 1) == size(c, 2) == length(n) || throw(ArgumentError("matrix-vector size mismatch"))

    # init dsfmt generator without altering GLOBAL_RNG
    rng = getRNG(seed)
    N = sum(n)
    K = length(n)
    nedg = zeros(Int, K, K)
    g = SimpleGraph(N)
    cum = [sum(n[1:a]) for a = 0:K]
    for a = 1:K
        ra = (cum[a] + 1):cum[a + 1]
        for b = a:K
            ((a == b) && !(c[a, b] <= n[b] - 1)) || ((a != b) && !(c[a, b] <= n[b])) &&
                error("Mean degree cannot be greater than available neighbors in the block.") # TODO 0.7: turn into some other error?

            m = a == b ? div(n[a] * (n[a] - 1), 2) : n[a] * n[b]
            p = a == b ? n[a] * c[a, b] / (2m) : n[a] * c[a, b] / m
            nedg = randbn(m, p, seed)
            rb = (cum[b] + 1):cum[b + 1]
            i = 0
            while i < nedg
                source = rand(rng, ra)
                dest = rand(rng, rb)
                if source != dest
                    if add_edge!(g, source, dest)
                        i += 1
                    end
                end
            end
        end
    end
    return g
end

"""
    stochastic_block_model(cint, cext, n)

Return a Graph generated according to the Stochastic Block Model (SBM), sampling
from an SBM with ``c_{a,a}=cint``, and ``c_{a,b}=cext``.
"""
function stochastic_block_model(cint::T, cext::T, n::Vector{U}; seed::Int=-1) where T <: Real where U <: Integer
    K = length(n)
    c = [ifelse(a == b, cint, cext) for a = 1:K, b = 1:K]
    stochastic_block_model(c, n, seed=seed)
end

"""
    StochasticBlockModel{T,P}

A type capturing the parameters of the SBM.
Each vertex is assigned to a block and the probability of edge `(i,j)`
depends only on the block labels of vertex `i` and vertex `j`.

The assignement is stored in nodemap and the block affinities a `k` by `k`
matrix is stored in affinities.

`affinities[k,l]` is the probability of an edge between any vertex in
block `k` and any vertex in block `l`.

### Implementation Notes
Graphs are generated by taking random ``i,j ∈ V`` and
flipping a coin with probability `affinities[nodemap[i],nodemap[j]]`.
"""
mutable struct StochasticBlockModel{T <: Integer,P <: Real}
    n::T
    nodemap::Array{T}
    affinities::Matrix{P}
    rng::MersenneTwister
end

==(sbm::StochasticBlockModel, other::StochasticBlockModel) =
    (sbm.n == other.n) && (sbm.nodemap == other.nodemap) && (sbm.affinities == other.affinities)

# A constructor for StochasticBlockModel that uses the sizes of the blocks
# and the affinity matrix. This construction implies that consecutive
# vertices will be in the same blocks, except for the block boundaries.
function StochasticBlockModel(sizes::AbstractVector, affinities::AbstractMatrix; seed::Int=-1)
    csum = cumsum(sizes)
    j = 1
    nodemap = zeros(Int, csum[end])
    for i in 1:csum[end]
        if i > csum[j]
            j += 1
        end
        nodemap[i] = j
    end
    return StochasticBlockModel(csum[end], nodemap, affinities, getRNG(seed))
end

### TODO: This documentation needs work. sbromberger 20170326
"""
    sbmaffinity(internalp, externalp, sizes)

Produce the sbm affinity matrix with internal probabilities `internalp`
and external probabilities `externalp`.
"""
function sbmaffinity(internalp::Vector{T}, externalp::Real, sizes::Vector{U}) where T <: Real where U <: Integer
    numblocks = length(sizes)
    numblocks == length(internalp) || throw(ArgumentError("Inconsistent input dimensions: internalp, sizes"))
    B = diagm(0 => internalp) + externalp * (ones(numblocks, numblocks) - I)
    return B
end

function StochasticBlockModel(internalp::Real,
                              externalp::Real,
                              size::Integer,
                              numblocks::Integer;
                              seed::Int=-1)
    sizes = fill(size, numblocks)
    B = sbmaffinity(fill(internalp, numblocks), externalp, sizes)
    StochasticBlockModel(sizes, B, seed=seed)
end

function StochasticBlockModel(internalp::Vector{T}, externalp::Real,
    sizes::Vector{U}; seed::Int=-1) where T <: Real where U <: Integer
    B = sbmaffinity(internalp, externalp, sizes)
    return StochasticBlockModel(sizes, B, seed=seed)
end

const biclique = ones(2, 2) - Matrix{Float64}(I, 2, 2)

#TODO: this documentation needs work. sbromberger 20170326
"""
    nearbipartiteaffinity(sizes, between, intra)

Construct the affinity matrix for a near bipartite SBM.
`between` is the affinity between the two parts of each bipartite community.
`intra` is the probability of an edge within the parts of the partitions.

This is a specific type of SBM with ``\\frac{k}{2} blocks each with two halves.
Each half is connected as a random bipartite graph with probability `intra`
The blocks are connected with probability `between`.
"""
function nearbipartiteaffinity(sizes::Vector{T}, between::Real, intra::Real) where T <: Integer
    numblocks = div(length(sizes), 2)
    return kron(between * Matrix{Float64}(I, numblocks, numblocks), biclique) + Matrix{Float64}(I, 2 * numblocks, 2 * numblocks) * intra
end

#Return a generator for edges from a stochastic block model near-bipartite graph.
nearbipartiteaffinity(sizes::Vector{T}, between::Real, inter::Real, noise::Real) where T <: Integer =
    nearbipartiteaffinity(sizes, between, inter) .+ noise

nearbipartiteSBM(sizes, between, inter, noise; seed::Int=-1) =
    StochasticBlockModel(sizes, nearbipartiteaffinity(sizes, between, inter, noise), seed=seed)

"""
    random_pair(rng, n)

Generate a stream of random pairs in `1:n` using random number generator `RNG`.
"""
function random_pair(rng::AbstractRNG, n::Integer)
    f(ch) = begin
        while true
            put!(ch, SimpleEdge(rand(rng, 1:n), rand(rng, 1:n)))
        end
    end
    return f
end

"""
    make_edgestream(sbm)

Take an infinite sample from the Stochastic Block Model `sbm`.
Pass to `Graph(nvg, neg, edgestream)` to get a Graph object based on `sbm`.
"""
function make_edgestream(sbm::StochasticBlockModel)
    pairs = Channel(random_pair(sbm.rng, sbm.n), ctype=SimpleEdge, csize=32)
    edges(ch) = begin
        for e in pairs
            i, j = Tuple(e)
            i == j && continue
            p = sbm.affinities[sbm.nodemap[i], sbm.nodemap[j]]
            if rand(sbm.rng) < p
                put!(ch, e)
            end
        end
    end
    return Channel(edges, ctype=SimpleEdge, csize=32)
end

"""
    SimpleGraph{T}(nv, ne, edgestream::Channel)

Construct a `SimpleGraph{T}` with `nv` vertices and `ne` edges from `edgestream`.
Can result in less than `ne` edges if the channel `edgestream` is closed prematurely.
Duplicate edges are only counted once.
The element type is the type of `nv`.
"""
function SimpleGraph(nvg::Integer, neg::Integer, edgestream::Channel)
    g = SimpleGraph(nvg)
    # println(g)
    for e in edgestream
        add_edge!(g, e)
        ne(g) >= neg && break
    end
    # println(g)
    return g
end

"""
    SimpleGraph{T}(nv, ne, smb::StochasticBlockModel)

Construct a random `SimpleGraph{T}` with `nv` vertices and `ne` edges.
The graph is sampled according to the stochastic block model `smb`.
The element type is the type of `nv`.
"""
SimpleGraph(nvg::Integer, neg::Integer, sbm::StochasticBlockModel) =
    SimpleGraph(nvg, neg, make_edgestream(sbm))

#TODO: this documentation needs work. sbromberger 20170326
"""
    blockcounts(sbm, A)

Count the number of edges that go between each block.
"""
function blockcounts(sbm::StochasticBlockModel, A::AbstractMatrix)
    I = collect(1:sbm.n)
    J =  [sbm.nodemap[i] for i in 1:sbm.n]
    V =  ones(sbm.n)
    Q = sparse(I, J, V)
    return (Q'A) * Q
end

function blockcounts(sbm::StochasticBlockModel, g::AbstractGraph)
    return blockcounts(sbm, adjacency_matrix(g))
end

function blockfractions(sbm::StochasticBlockModel, g::Union{AbstractGraph,AbstractMatrix})
    bc = blockcounts(sbm, g)
    bp = bc ./ sum(bc)
    return bp
end

"""
    kronecker(SCALE, edgefactor, A=0.57, B=0.19, C=0.19)

Generate a directed [Kronecker graph](https://en.wikipedia.org/wiki/Kronecker_graph)
with the default Graph500 parameters.

###
References
- http://www.graph500.org/specifications#alg:generator
"""
function kronecker(SCALE, edgefactor, A=0.57, B=0.19, C=0.19)
    N = 2^SCALE
    M = edgefactor * N
    ij = ones(Int, M, 2)
    ab = A + B
    c_norm = C / (1 - (A + B))
    a_norm = A / (A + B)

    for ib = 1:SCALE
        ii_bit = rand(M) .> (ab)  # bitarray
        jj_bit = rand(M) .> (c_norm .* (ii_bit) + a_norm .* .!(ii_bit))
        ij .+= 2^(ib - 1) .* (hcat(ii_bit, jj_bit))
    end

    p = randperm(N)
    ij = p[ij]

    p = randperm(M)
    ij = ij[p, :]

    g = SimpleDiGraph(N)
    for (s, d) in zip(@view(ij[:, 1]), @view(ij[:, 2]))
        add_edge!(g, s, d)
    end
    return g
end

"""
    dorogovtsev_mendes(n)

Generate a random `n` vertex graph by the Dorogovtsev-Mendes method (with `n \\ge 3`).

The Dorogovtsev-Mendes process begins with a triangle graph and inserts `n-3` additional vertices.
Each time a vertex is added, a random edge is selected and the new vertex is connected to the two
endpoints of the chosen edge. This creates graphs with a many triangles and a high local clustering coefficient.

It is often useful to track the evolution of the graph as vertices are added, you can access the graph from
the `t`th stage of this algorithm by accessing the first `t` vertices with `g[1:t]`.

### References
- http://graphstream-project.org/doc/Generators/Dorogovtsev-Mendes-generator/
- https://arxiv.org/pdf/cond-mat/0106144.pdf#page=24
"""
function dorogovtsev_mendes(n::Integer; seed::Int=-1)
    n < 3 && throw(DomainError("n=$n must be at least 3"))
    rng = getRNG(seed)
    g = CycleGraph(3)

    for iteration in 1:(n-3)
        chosenedge = rand(rng, 1:(2*ne(g))) # undirected so each edge is listed twice in adjlist
        u, v = -1, -1
        for i in 1:nv(g)
            edgelist = outneighbors(g, i)
            if chosenedge > length(edgelist)
                chosenedge -= length(edgelist)
            else
                u = i
                v = edgelist[chosenedge]
                break
            end
        end

        add_vertex!(g)
        add_edge!(g, nv(g), u)
        add_edge!(g, nv(g), v)
    end
    return g
end

"""
    random_orientation_dag(g)

Generate a random oriented acyclical digraph. The function takes in a simple
graph and a random number generator as an argument. The probability of each
directional acyclic graph randomly being generated depends on the architecture
of the original directed graph.

DAG's have a finite topological order; this order is randomly generated via "order = randperm()".
"""
function random_orientation_dag(g::SimpleGraph{T}, seed::Int=-1) where T <: Integer
    nvg = length(g.fadjlist)
    rng = getRNG(seed)
    order = randperm(rng, nvg)
    g2 = SimpleDiGraph(nv(g))
    @inbounds for i in vertices(g)
        for j in outneighbors(g, i)
            if order[i] < order[j]
                add_edge!(g2, i, j)
            end
        end
    end
    return g2
end