text:
using Documenter
using Plots
using CapacityExpansion

makedocs(
    sitename = "CapacityExpansion.jl",
    authors = "Elias Kuepper, and Holger Teichgraeber",
    pages = [
        "Introduction" => "index.md",
        "Quickstart" => "quickstart.md",
        "Workflow" => "workflow.md",
        "Data" => [
            "preparing_clust_data.md",
            "preparing_opt_cep_data.md",
            "csv_structure.md",
            "README_GER_18.md",
            "README_GER_1.md",
            "README_CA_14.md",
            "README_CA_1.md",
            "README_TX_1.md",
        ],
        "Optimization" => ["opt_cep.md", "results_opt.md", "opt_cep_examples.md"],
        "Teaching" => ["teaching.md", "homework.md"],
        "Testing" => "testing.md",
    ],
    format = Documenter.HTML(
        assets = [
            asset("assets/cep_text.svg", class = :ico, islocal = true),
            asset("assets/opt_cep.svg", class = :ico, islocal = true),
            asset("assets/workflow.svg", class = :ico, islocal = true),
            asset("assets/GER_1.svg", class = :ico, islocal = true),
            asset("assets/GER_18.svg", class = :ico, islocal = true),
            asset("assets/CA_1.svg", class = :ico, islocal = true),
            asset("assets/CA_14.svg", class = :ico, islocal = true),
            asset("assets/opt_cep_cap_plot.svg", class = :ico, islocal = true),
            asset("assets/preparing_clust_data_load.svg", class = :ico, islocal = true),
            asset("assets/preparing_clust_data_agg.svg", class = :ico, islocal = true),
        ],
    ),
)

deploydocs(repo = "github.com/YoungFaithful/CapacityExpansion.jl.git", devbranch = "dev")


text:
#push!(LOAD_PATH, "../src/")

using Documenter, ImagePhaseCongruency

makedocs(
    sitename = "ImagePhaseCongruency",
    pages = [
        "index.md",
        "examples.md",
        "functions.md"
    ]
)

deploydocs(
    repo = "github.com/peterkovesi/ImagePhaseCongruency.jl.git",
)

ext: jl


text:
using Sockets

mutable struct Node
    host::Union{String,IPv4}
    port::Integer
    testnet::Bool
    logging::Bool
    acknowledged::Bool
    sock::TCPSocket
    Node(host::Union{String,IPv4}, port::Integer, testnet::Bool=false, logging::Bool=false) = new(host, port, testnet, logging, false)
end

Node(host::Union{String,IPv4}, testnet::Bool=false) = Node(host, DEFAULT["port"][testnet], testnet)

"""
    connect!(node::Node) -> TCPSocket

Connect to the host `node.host` on port `node.port`.
"""
function connect!(node::Node)
    try
        if !isopen(node.sock)
            node.sock = connect(node.host, node.port)
        else
            node.sock
        end
    catch
        node.sock = connect(node.host, node.port)
    end
end

"""
    close!(node::Node) -> TCPSocket

Close Node's TCPSocket
"""
function close!(node::Node)
    close(node.sock)
    node.sock
end

"""
    send2node(node::Node, message::T) -> Integers

Send a message to the connected node, returns the numbers of bytes sent
"""
function send2node(node::Node, message::T) where {T<:AbstractMessage}
    envelope = NetworkEnvelope(message.command,
                               serialize(message),
                               node.testnet)
    if node.logging
        println("Sending : ", envelope)
    end
    write(node.sock, serialize(envelope))
end

"""
handshake(node::Node) -> Bool

Do a handshake with the other node, returns true if successful
Handshake is sending a version message and getting a verack back.
"""
function handshake(node::Node, messages::Channel)
    println("Connecting to node...")
    connect!(node)
    @async read(node.sock)
    @async read_messages(node, messages)
    println("Sending version message...")
    send2node(node, VersionMessage())
    i = 0
    while !node.acknowledged
        if eof(node.sock) || i == 2
            println("Failed handshake!")
            return false
        end
        msg = take!(messages)
        answer(node, msg)
        i += 1
    end
    println("Handshake successfull")
    return true
end

"""
    getheaders(node::Node, stop::Integer, start::Integer=1) -> Vector{BlockHeader}

Returns a list of blockheaders, from `start` to `stop` height
!!!Experimental function, not tested as it should
"""
function getheaders(node::Node, stop::Integer, start::Integer=1)
    handshake(node)
    last_block_hash = GENESIS_BLOCK_HASH[node.testnet]
    current_height = start
    result = BlockHeader[]
    while current_height <= stop
        try
            msg = GetHeadersMessage(last_block_hash)
            send2node(node, msg)
            if bytesavailable(node.sock) > 0
                raw = read(node.sock.buffer)
                envelopes = io2envelopes(raw, node.testnet)
                for envelope in envelopes
                    if envelope.command == "headers"
                        headers = PARSE_PAYLOAD[envelope.command](envelope.payload)
                        for header in headers.headers
                            if !check_pow(header)
                                error("bad proof of work at block ", id(header))
                            end
                            if (last_block_hash != GENESIS_BLOCK_HASH) && (header.prev_block != last_block_hash)
                                error("discontinuous block at ", id(header))
                            end
                            if current_height % 2016 == 0
                                println(id(header))
                            end
                            last_block_hash = hash(header)
                            current_height += 1
                            push!(result, header)
                            if node.logging
                                println(header)
                            end
                        end
                    end
                end
            else
                sleep(0.01)
            end
        catch e
            if typeof(e) == InterruptException
                return error("Interrupted")
            end
            println("Error ", e, " was raised, retrying...")
            sleep(1)
        end
    end
    result
end

function read_messages(node::Node, net_msg::Channel)
    while !eof(node.sock)
        envelope = io2envelope(node.sock.buffer, node.testnet)
        command = envelope.command
        msg = haskey(PARSE_PAYLOAD, command) ?
              PARSE_PAYLOAD[command](envelope.payload) :
              bytes2hex(envelope.payload)
        put!(net_msg, msg)
    end
end

function answer(node::Node, msg::HeadersMessage)
    println("\n\n --- Headers received --- \n\n")
    getdata = GetDataMessage()
    for b in msg.headers
        if check_pow(b)
            append!(getdata, DATA_MESSAGE_TYPE["MSG_FILTERED_BLOCK"], hash(b))
        else
            error("Invalid proof of work for block ", b)
        end
    end
    send2node(node, getdata)
end

function answer(node::Node, msg::MerkleBlockMessage)
    if !is_valid(msg)
        error("Merkle proof is *not* valid!")
    else
        println("Merkle proof is valid.")
    end
end

function answer(node::Node, msg::VersionMessage)
    send2node(node, VerAckMessage())
end

function answer(node::Node, msg::VerAckMessage)
    node.acknowledged = true
    println("\n--------\nVersion Acknowledgement received from\n", node)
    return
end

function answer(node::Node, msg::PingMessage)
    send2node(node, PongMessage(serialize(msg)))
end

function answer(node::Node, msg::AbstractMessage)
    println("Ignoring...")
end

function get_tx_of_interest(node::Node, net_msg::Channel, adr::String)
    found = false
    while !found
        msg = take!(net_msg)
        if typeof(msg) <: Tx
            println("\n\n --- Tx received --- \n\n")
            i = 0
            for tx_out in msg.tx_outs
                if script2address(tx_out.script_pubkey, true) == adr
                    println("found: ", txid(msg), " ", i)
                    found = true
                    break
                end
                i += 1
            end
        else
            println("Processing ", msg)
            answer(node, msg)
        end
    end
end


text:
module Patches

using AWS
using Dates
using Downloads: Downloads
using HTTP
using JSON
using GitHub
using Mocking
using OrderedCollections: LittleDict

version = v"1.1.0"
status = 200
headers = Pair[
    "x-amz-id-2" => "x-amz-id-2",
    "x-amz-request-id" => "x-amz-request-id",
    "Date" => "Tue, 16 Jun 2020 21:29:18 GMT",
    "x-amz-bucket-region" => "us-east-1",
    "Content-Type" => "application/xml",
    "Transfer-Encoding" => "chunked",
    "Server" => "AmazonS3",
]

body = """
    <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n
    <ListBucketResult xmlns=\"http://s3.amazonaws.com/doc/2006-03-01/\">
        <Name>sample-bucket</Name>
        <Prefix></Prefix>
        <Marker></Marker>
        <MaxKeys>1000</MaxKeys>
        <IsTruncated>false</IsTruncated>
        <Contents>
            <Key>test.txt</Key>
            <LastModified>2020-06-16T21:37:34.000Z</LastModified>
            <ETag>"d41d8cd98f00b204e9800998ecf8427e"</ETag>
            <Size>0</Size>
            <Owner>
                <ID>id</ID>
                <DisplayName>matt.brzezinski</DisplayName>
            </Owner>
            <StorageClass>STANDARD</StorageClass>
        </Contents>
    </ListBucketResult>
    """

function _response(;
    version::VersionNumber=version,
    status::Int64=status,
    headers::Array=headers,
    body::String=body,
)
    response = HTTP.Messages.Response()

    response.version = version
    response.status = status
    response.headers = headers
    response.body = b"[Message Body was streamed]"

    b = IOBuffer(body)

    return AWS.Response(response, b)
end

function _aws_http_request_patch(response::AWS.Response=_response())
    p = @patch AWS._http_request(::AWS.AbstractBackend, request::Request, ::IO) = response
    return p
end

_cred_file_patch = @patch function dot_aws_credentials_file()
    return ""
end

_config_file_patch = @patch function dot_aws_config_file()
    return ""
end

_assume_role_patch = function (
    op;
    access_key="access_key",
    secret_key="secret_key",
    session_token="token",
    role_arn="arn:aws:sts:::assumed-role/role-name",
    expiry=duration -> now(UTC) + duration,
    token_code_ref=nothing,
)
    @patch function AWSServices.sts(op, params; aws_config, feature_set)
        duration = Second(parse(Int, get(params, "DurationSeconds", "3600")))
        expiration = expiry(duration)
        if token_code_ref !== nothing
            token_code_ref[] = params["TokenCode"]
        end
        xml = """
            <$(op)Response xmlns="https://sts.amazonaws.com/doc/2011-06-15/">
              <$(op)Result>
                <AssumedRoleUser>
                  <Arn>$(role_arn)/$(params["RoleSessionName"])</Arn>
                </AssumedRoleUser>
                <Credentials>
                  <AccessKeyId>$access_key</AccessKeyId>
                  <SecretAccessKey>$secret_key</SecretAccessKey>
                  <SessionToken>$session_token</SessionToken>
                  <Expiration>$expiration</Expiration>
                </Credentials>
              </$(op)Result>
            </$(op)Response>
            """

        r = _response(; body=xml)
        return feature_set.use_response_type ? r : parse(r)::AbstractDict
    end
end

_getpass_patch = function (; secret="the_secret")
    @patch function Base.getpass(prompt)
        return Base.SecretBuffer(secret)
    end
end

_github_tree_patch = @patch function tree(repo, tree_obj; kwargs...)
    if tree_obj == "master"
        tree = [Dict("path" => "apis", "sha" => "apis-sha")]
        return Tree("test-sha", HTTP.URI(), tree, false)
    else
        tree = [Dict("path" => "test-2020-01-01.normal.json", "sha" => "test-sha")]
        return Tree("test-sha", HTTP.URI(), tree, false)
    end
end

_instance_metadata_timeout_patch = @patch function HTTP.request(
    method::String, url; kwargs...
)
    return throw(HTTP.ConnectionPool.ConnectTimeout("169.254.169.254", "80"))
end

# This patch causes `HTTP.request` to return all of its keyword arguments
# except `require_ssl_verification` and `response_stream`. This is used to
# test which other options are being passed to `HTTP.Request` inside of
# `_http_request`.
_http_options_patches = [
    @patch function HTTP.request(args...; kwargs...)
        options = Dict(kwargs)
        delete!(options, :require_ssl_verification)
        delete!(options, :response_stream)
        return options
    end
    @patch AWS.Response(options, args...) = options
]

get_profile_settings_empty_patch = @patch function aws_get_profile_settings(profile, ini)
    return nothing
end

get_profile_settings_patch = @patch function aws_get_profile_settings(profile, ini)
    return Dict("foo" => "bar")
end

# Simulate the HTTP.request behaviour with a HTTP 400 response
function gen_http_options_400_patches(message)
    body = "{\"__type\":\"AccessDeniedException\",\"Message\":\"$message\"}"
    headers = [
        "Content-Type" => "application/x-amz-json-1.1",
        "Content-Length" => string(sizeof(body)),
    ]

    return [
        @patch function HTTP.request(
            args...; status_exception=true, response_stream=nothing, kwargs...
        )
            request = HTTP.Request("GET", "/")

            if response_stream !== nothing
                write(response_stream, body)
                close(response_stream)  # Simulating current HTTP.jl 0.9.14 behaviour
                body = HTTP.MessageRequest.body_was_streamed
            end

            response = HTTP.Response(400, headers; body=body, request=request)
            exception = HTTP.StatusError(400, response)
            return !status_exception ? response : throw(exception)
        end
        @patch function Downloads.request(args...; output=nothing, kwargs...)
            if output !== nothing
                write(output, body)
            end

            return Downloads.Response(
                "https",
                "https://region.amazonaws.com/",
                400,
                "HTTP/1.1 400 Bad Request",
                headers,
            )
        end
    ]
end

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -2.3024137037198598 <= q <= 2.756640266423038)

                     add_NL_constraint(m, :(log(1 + exp(0.8745011162287222 + -0.4747792679508138*log(1 + exp(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2]))) + 0.06822380173369114*log(1 + exp(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))))) + log(1 + exp(-0.5284568968024108 + 0.8782860084347535*log(1 + exp(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2]))) + 0.990986450920027*log(1 + exp(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#==============================================================================#
# RJNN: Random network and data generation
#==============================================================================#

using Random

"Returns a neural network with random edges and weights."
function randnet(ins::Integer, outs::Integer, hidden::Vector{<:Integer};
        batch_size=1024, density=0.5)
    sizes = [ins, hidden..., outs]
    dims = zip(sizes[1:end-1], sizes[2:end])
    weights = [zeros(m, n) for (m, n) in dims]
    biases = [zeros(1, n) for n in sizes[2:end]]

    # Initialise a random subset of the weights and biases, and ensure that
    # every unit is connected to its neighbouring layers by at least one edge --
    # i.e., that the weight matrices have no empty rows or columns.
    for (W, b) in zip(weights, biases)
        randnfill!(W, density)
        randnfill!(b, density)
        randnconnect!(W)
    end

    Net(weights, biases, batch_size)
end

"Initialises a random subset of an array's elements using a standard normal
 distribution."
function randnfill!(A::AbstractArray{<:AbstractFloat}, p)
    js = randsubseq(eachindex(A), p) # each position is chosen with independent
    randn!(@view A[js])              # probability p (Bernoulli sampling).
end

"Ensures that every row and column of a matrix contains at least one non-zero
 value."
function randnconnect!(A::AbstractMatrix{<:AbstractFloat})
    for i in axes(A, 1)
        if iszero(@view A[i,:])
            A[i, rand(axes(A, 2))] = randn()
        end
    end
    for j in axes(A, 2)
        if iszero(@view A[:,j])
            A[rand(axes(A, 1)), j] = randn()
        end
    end
end

max_stars_repo_name: kevdur/rjnn

max_issues_repo_name: kevdur/rjnn

max_forks_repo_name: kevdur/rjnn


text:
const data_dir    = joinpath("..", "data")
const result_dir  = joinpath("..", "result")
const model_dir   = joinpath("..", "model")
const raw_samples = joinpath(data_dir, "raw_samples")
const sample_meta = joinpath(data_dir, "sample_meta")
const raw_field_datasets = joinpath(data_dir,"raw_field_datasets")

const miss = "MISSING"
const test = false

immutable Atom
    idx_row::Int64
    idx_col::Int64
    field_val::ASCIIString
end
function row(atom::Atom)
    atom.idx_row
end
function col(atom::Atom)
    atom.idx_col
end
function field_val(atom::Atom)
    atom.field_val
end
function convert(Sample, atom::Atom)
    Sample(atom.idx_row, Int64[atom.idx_col], ASCIIString[field_val(atom)])
end
function col(atom_1::Atom,atom_2::Atom)
    Int64[atom_1.idx_col, atom_2.idx_col]
end
function field_val(atom_1::Atom, atom_2::Atom)
    ASCIIString[atom_1.field_val, atom_2.field_val]
end
function Atom(idx_row,idx_col,field_val::ASCIIString)
    @assert idx_row > 0
    @assert idx_col > 0
    if length(field_val) == 0
        field_val = miss
    end
    Atom(idx_row, idx_col, field_val)
end

immutable Sample
    idx_row::Int64
    idx_col::Array{Int64,1}
    field_val::Array{ASCIIString,1}
end
function field_val(sample::Sample)
    sample.field_val
end
function row(sample::Sample)
    sample.idx_row
end
function col(sample::Sample)
    sample.idx_col
end
function Sample(idx_row,idx_col,field_val)
    @assert idx_row > 0
    @assert idx_col > 0
    if length(field_val) == 0
        field_val = miss
    end
    Sample(idx_row, idx_col, field_val)
end
function Sample(atom_1::Atom, atom_2::Atom)
    @assert row(atom_1) == row(atom_2)
    Sample(row(atom_1), col(atom_1,atom_2),field_val(atom_1,atom_2))
end
function Sample(sample::Sample, atom::Atom)
    @assert row(sample) == row(atom)
    idx_cols = vcat(col(sample),col(atom))
    field_vals = vcat(field_val(sample),field_val(atom))
    Sample(row(sample),idx_cols,field_vals)
end
function Sample(sample_1::Sample, sample_2::Sample)
    @assert row(sample_1) == row(sample_2)
    idx_cols = vcat(col(sample_1),col(sample_2))
    field_vals = vcat(field_val(sample_1),field_val(sample_2))
    Sample(row(sample_1),idx_cols,field_vals)
end

function show(sample::Sample)
    println(sample.idx_row,sample.idx_col)
    println(sample.field_val)
end
function tet()
end
function row(samples::Array{Sample})
    map(row, samples)
end
type OrderedSampleMatrix
    num_samples
    samples::Array{Sample} # sorted idx_IDsample
end
function OrderedSampleMatrix(samples::Array{Sample}) # for null ordered sample matrix
    idxs = row(samples)
    num_samples = length(samples)
    perm = sortperm(idxs)
    OrderedSampleMatrix(num_samples, samples[perm])
end
function start(orspmat::OrderedSampleMatrix)
    1
end
function start(sp::Sample)
    1
end
function next(sp::Sample, state)
    sp.samples[state], state+1
end
function done(orspmat::OrderedSampleMatrix, state)
    state > orspmat.num_samples
end
function convert(::Type{Array{Sample,1}}, samples::Array{Any,1})
    ret = Sample[]
    for sample in samples
        push!(ret, sample)
    end
    ret
end

#@doc """  build selected fields into OrderedSampleMatrices at field_dir
#""" ->
function build_field_datasets()
    header = readdlm(joinpath(sample_meta, "header.tsv"), '\t', ASCIIString)
    header = header[:,1]
    @show header
    deleteat!(header, 1)
    field_names = [(name,idx) for (idx,name) in enumerate(header)]
    map(build_field_dataset, field_names)
end
function build_field_dataset(field_meta::Tuple{ASCIIString,Int64})
    build_field_dataset(field_meta[1],field_meta[2])
end
#@doc """ the simplist way of build machine learning dataset only using one feature: Mutation_ID
#         ID_sample: [(Gene_name,feat),(),...]
#         hash(Gene_name) --> col_idx
#         hash(Sample_ID) --> row_idx
#""" ->
function build_field_dataset(field_name::ASCIIString, idx_field::Int64)
    field_dir = joinpath(data_dir, "field")
    !isdir(field_dir) && mkdir(field_dir)

    # create hashes for sample_ID and single_field

    ### genes:   all the genes then filter
    ### samples: all the samples then filter
    ### field: single_field
    sample_meta   = joinpath(data_dir, "sample_meta")
    !isdir(sample_meta) && @error("can't find sample_meta")

    field_vals = readcsv(joinpath(sample_meta, string(field_name, ".csv")))
    gene_names = readcsv(joinpath(sample_meta, "Gene name.csv"), ASCIIString)
    datasets = ["train.csv","val.csv","test.csv"]
    num = 0
    for dataset in datasets
        println("dataset is ", dataset)
        sample_fls = readcsv(joinpath(sample_meta,dataset))
        ID_samples = map(x->x[1:end-4],sample_fls)
        hash_sample = Dict()
        for idx = 1:length(sample_fls)
            hash_sample[ID_samples[idx]] = idx
        end

        #the name of columns is Gene_name

        hash_gene = Dict()
        for idx = 1:length(gene_names)
            hash_gene[gene_names[idx]] = idx
        end
        @show hash_gene["39340"]
        function merge_atom(atom_1::Atom,atom_2::Atom)
            Sample(atom_1, atom_2)
        end
        function merge_atom(atom_1::Sample,atom_2::Atom)
            Sample(atom_1, atom_2)
        end
        function merge_atom(atom_1::Sample,atom_2::Sample)
            Sample(atom_1, atom_2)
        end

        samples = @parallel vcat for sample_fl in sample_fls
            f = open(joinpath(raw_samples, sample_fl))
            lines = readlines(f)
            function atom_mapper(line)
                line = strip(line,'\n')
                row = split(line, '\t')
                @assert length(row) == 34
                # row[6] is ID_sample
                row_idx = try
                            hash_sample[row[6]]
                          catch e
                            println("==========Error Info===========")
                            println(e)
                            println(row)
                          end

                col_idx = try
                            hash_gene[strip(row[1])]
                          catch e
                            println("==========Error Info===========")
                            println(e)
                            println(row)
                          end

                # Mutation_ID:#17 row_idx col_idx
                atom = Atom(row_idx, col_idx, row[idx_field])
            end
            sample = reduce(merge_atom, pmap(atom_mapper, lines)) # short filea
            if typeof(sample) <: Atom
                sample = convert(Sample, sample)
            end
            close(f)
            sample
        end
        @debug("here ========")
        sample_matrix = OrderedSampleMatrix(samples)
        dataset_dir = joinpath(field_dir, filename(dataset))
        isdir(dataset_dir) || mkdir(dataset_dir)
        dataset_name = joinpath(dataset_dir, string(field_name,".tsv"))
        @debug("bugs above ====")
        jldopen(dataset_name, "w") do file
            println("save $dataset into file $dataset_name")
            write(file, "sample_matrix", sample_matrix)
        end
    end
end

@doc """ Since cosmic raw data file is large, we tend to build all the fields with only one exaust search
""" ->
function build_fields_meta(cosmic_path::ASCIIString)
    f = open(cosmic_path)
    header = readline(f)
    header = split(header, '\t')
    data   = readdlm(f, ASCIIString)
    map!(x -> length(x) == 0 ? miss : x, data)

    if isdir(sample_meta)

    else
        mkdir(sample_meta)
    end
    for (idx,field) in enumerate(header) #@simd
        writecsv(joinpath(sample_meta, string(strip(header[idx]),".csv")), sort(unique(data[:,idx])))
    end

end

@doc """ random build raw train, validation, evaluation dataset based on the raw_samples
""" ->
function split_samples(;ratios=(0.6,0.2,0.2))
    raw_samples = joinpath(data_dir, "raw_samples")
    sample_fls = try
        readdir(raw_samples)
    catch
        @warn("$raw_samples is not exits, program terminated") #TODO
        exit(-1)
    end
    num_samples = length(sample_fls) - 1
    @assert length(readdir(raw_samples)) > 0

    tr_val_point  = round(Int64, num_samples * ratios[1])
    val_te_point  = round(Int64, num_samples * (ratios[1] + ratios[2]))

    shuffle!(sample_fls)
    train_samples = sample_fls[1:tr_val_point]
    val_samples   = sample_fls[tr_val_point+1:val_te_point]
    te_samples    = sample_fls[val_te_point+1:end]

    sample_meta   = joinpath(data_dir, "sample_meta")
    if !isdir(sample_meta)
        mkdir(sample_meta)
    end
    train_fl,val_fl,test_fl = "train.csv","val.csv","test.csv"
    for (fl,samples) in [(train_fl,train_samples), (val_fl,val_samples),(test_fl, te_samples)]
        fl = joinpath(sample_meta, fl)
        writecsv(fl, samples)
    end
end

@doc """ build raw samples to disk: ../data/raw_samples from raw cosmic file
""" ->
function build_raw_samples(cosmic_path::ASCIIString)
    raw_samples = joinpath(data_dir,"raw_samples")
    if !isdir(raw_samples)
        mkdir(raw_samples)
    else
        println("removing files under $raw_samples")
        rm(raw_samples,recursive=true)
        mkdir(raw_samples)
    end
    f = open(cosmic_path)
    header = readline(f)
    write(joinpath(sample_meta, "header.tsv"), split(strip(header, '\t')))
    for line in eachline(f)
        row = split(line,'\t')
        sample_name = strip(row[6]) # #row[6] is ID_sample
        f_sp = open(joinpath(raw_samples,string(sample_name,".tsv")), "a+")
        write(f_sp, line)
        close(f_sp)
    end
    close(f)
end

function build_dir()
    for dir in (data_dir,result_dir,model_dir)
        if !isdir(dir)
            println("buiding directory $data_dir")
            mkdir(dir)
        end
    end
end

function filename(file::ASCIIString)
    split(file,".")[1]
end


text:
import StateSpaceReconstruction:
	DelaunayTriangulation

@recipe function plot_triang(original_pts, DT::DelaunayTriangulation;
            plot_states = true,
            plot_simplices = true,
            evenly_subsample = true,
            k = 3, # The splitting factor, k^dim new points are introduced per simplex
            n_interior_pts = 100,
            insert_pts = false,
            ms = 0.5,
            ms_originalpts = 1.5,
            mc = :black,
            mc_projection = :red,
            mc_originalpts = :black,
            mα = 0.5,
            edges = true,
            vertices = false,
            si = 0, # index of the i'th highlighted simplex
            sj = 0, # index of the j'th highlighted simplex
            ind_projection = 0, #show the simplex corresponding to the forward projection of the `ind_projection`-th simplex
            intersecting_vertices = nothing,
            plot_intersecting_pts = false,
            highlight_ms = 2,
            extrapts = [],
            lc = :black,
            plot_originalpts = true,
            plot_triang = true,
            plot_imagetriang = true,
            plot_fulltriang = false,
            lc_fulltriang = :black,
            lw_fulltriang = 0.5,
            lα_fulltriang = 0.5,
            lc_triang = :black,
            lc_imagetriang = :black,
            lw_imagetriang = 0.5,
            lw_triang = 0.5,
            ind_excluded_pt = 0,
            lw_highlight = 1.2,
            markershape_excluded_pt = :star5,
            mc_excluded_pt = :black,
            ms_excluded_pt = 3,
            lα = 0.2,
            lw = 0.5,
            lc_i = :green,
            lc_j = :red,
            lw_i = 0.5,
            lw_j = 0.5,
            lα_i = 1.0,
            lα_j = 1.0,
            mc_j = :red,
            mc_i = :green,
            mw_i = 0.5,
            mw_j = 0.5,
            mα_i = 1.0,
            mα_j = 1.0,
            lc_projection = :blue,
            lw_projection = 2,
            lα_projection = 1,
            idx_intersect_simplex = si, # the index of the simplex for which to compute the intersection between the simplex and its forward image
            mc_intersecting_pts = :black,
            markershape_intersecting_pts = :pentagon,
            ms_intersecting_pts = 1,
            fill_intersecting_volume = true,
            n_fillpoints = 500,
            ms_fillpoints = 1.0,
            mα_fillpoints = 1.0,
            mc_fillpoints = :black,
            sample_randomly = true
    )
    fulltriang = DelaunayTriangulation(customembed(original_pts).reconstructed_pts)
    n_simplices_full = size(fulltriang.indices, 2)
    n_simplices = size(DT.indices, 2)
    n_points = maximum(size(pts))

    # Embedding center
    #ce = sum(r.points, 2)/size(r.points, 2)
    #lp = r.points[:, end]

    #scatter3d!(p, splitaxes(ce), ms = 3)
    #scatter3d!(p, splitaxes(lp), ms = 3, mc = :black)

    legend --> false
    xaxis --> false
    if length(extrapts) > 0
        @series begin
            seriestype := :scatter
            legend --> false
            mc --> mc
            ms --> ms
            mα --> mα
            splitaxes(extrapts)
        end
    end
    if plot_originalpts
        @series begin
            seriestype := :scatter
            mc --> mc_originalpts
            ms --> ms_originalpts
            splitaxes(original_pts)
        end
        if ind_excluded_pt > 0
            @series begin
                seriestype := :scatter
                mc --> mc_excluded_pt
                markershape --> markershape_excluded_pt
                ms --> ms_excluded_pt
                splitaxes(original_pts[:, ind_excluded_pt])
            end
        end
    end

    if plot_intersecting_pts && (ind_projection > 0)
        sᵢ = getsimplex(original_pts, DT, si)
        sⱼ = forwardmap(original_pts, DT, ind_projection)
        intvol = simplexintersection(sᵢ, sⱼ)

        if intvol > 0 + 1e-8
            intersecting_vertices, fillpoints = generate_fillpoints(sᵢ, sⱼ,
                            n_fillpoints, sample_randomly)

            @series begin
                seriestype := :scatter
                ms --> ms_intersecting_pts
                markershape --> markershape_intersecting_pts
                mc --> mc_intersecting_pts
                splitaxes(transpose(intersecting_vertices))
            end
            @series begin
                seriestype := :scatter
                ms --> ms_fillpoints
                mα --> mα_fillpoints
                mc --> mc_fillpoints
                lc --> mc_fillpoints
                fc --> mc_fillpoints
                splitaxes(fillpoints)
            end
        end
    end

    if plot_fulltriang
        for i = 1:n_simplices
            sᵢ = getsimplex(original_pts, fulltriang, i)

            @series begin
                seriestype := :path
                lc --> lc_fulltriang
                lw --> lw_fulltriang
                lα --> lα_fulltriang
                legend --> false
                preparesimplex(sᵢ)
            end
        end
    end

    for i = 1:n_simplices
        sᵢ = getsimplex(original_pts, DT, i)

        if edges
            if plot_triang
                @series begin
                    seriestype := :path
                    lc --> lc_triang
                    lw --> lw_triang
                    lα --> lα
                    legend --> false
                    preparesimplex(sᵢ)
                end
            end

            if plot_imagetriang
                @series begin
                    seriestype := :path
                    lc --> lc_imagetriang
                    lw --> lw_imagetriang
                    lα --> lα
                    ls --> :dash
                    legend --> false
                    preparesimplex(forwardmap(original_pts, DT, i))
                end
            end
        end

        if vertices
            if plot_triang
                @series begin
                    seriestype := :scatter
                    mc --> mc
                    ms --> ms
                    mα --> mα
                    legend --> false
                end
            end
        end

        if insert_pts
            if !evenly_subsample
                @series begin
                    mc --> mc
                    ms --> ms
                    lα --> lα
                    legend --> false
                    random_interior_pts = interiorpoints(sᵢ, n_interior_pts)
                    splitaxes(random_interior_pts)
                end
            elseif evenly_subsample
                @series begin
                    mc --> mc
                    ms --> ms
                    lα --> lα
                    legend --> false
                    splitaxes(sᵢ * even_sampling_rules(3, k))
                end
            end
        end
    end

    if si > 0
        sᵢ = getsimplex(original_pts, DT, si)

        if edges
            @series begin
                seriestype := :path
                label --> "simplex #$sᵢ"
                lw --> lw_highlight
                lc --> lc_i
                lα --> lα_i
                preparesimplex(sᵢ)
            end
        end

        if vertices
            @series begin
                label --> "simplex #$sᵢ"
                ms --> ms_i
                mc --> mc_i
                preparesimplex(sᵢ)
            end
        end
    end

    if sj > 0
        sⱼ = getsimplex(original_pts, DT, sj)

        if edges
            @series begin
                seriestype := :path
                label --> "simplex #$sⱼ"
                lw --> lw_highlight
                lc --> lc_j
                lα --> lα_j
                ls --> :dot
                preparesimplex(sⱼ)
            end
        end

        if vertices
            @series begin
                label --> "simplex #$sⱼ"
                ms --> ms_j
                mc --> mc_j
                preparesimplex(sⱼ)
            end
        end
    end

    if ind_projection > 0
        @series begin
            seriestype := :path
            label --> "projection of simplex #$ind_projection"
            lw --> lw_highlight
            lc --> lc_projection
            lα --> lα_projection
            preparesimplex(forwardmap(original_pts, DT, ind_projection))
        end
    end
end

"""
Visualize a set of N points (`original_pts`), a delaunay triangulation
(`DT`) constructed from and embedding (`E`) of the first ``N-1`` points
of the points. By assuming a triangulation of the ``N-1`` first points,
the recipe allows for plotting both the original triangulation and the
image of that triangulation under the forward linear map of its vertices,
one step forward in time.

## Arguments
- **``original_pts``**: A set of ``N`` points where each column is a point.
- **``E``**: An embedding constructed from the first ``N-1`` points of
    `original_pts`.
- **``DT``**: A DelaunayTriangulation instance constructed from `E`.

## Keyword arguments
- **``plot_triang``**: Plot the original triangulation?
- **``plot_imagetriang``**: Plot the image of the triangulation under the forward linear
    map of its vertices?
- **``plot_fulltriang``**: Plot a triangulation of all points in `original_pts`?

"""
plot_triang

export plot_triang


text:
module CausalDiscovery
using Reexport

include("Autumn/Autumn.jl")
@reexport using .Autumn

# include("CISC/CISC.jl")
# include("CISC.jl")
# include("MCMC.jl/model.jl")
# include("MCMC.jl/grammar.jl")

end # module


text:
# To execute all statistical tests execute this file from within this folder
# "julia execute_tests_from_paper.jl"

include("stat_tests_mutation.jl") # Mutation scenario
include("stat_tests_between_scenarios.jl") # Mutation vs. no mutation
include("stat_tests_no_mutation.jl") # No mutation scenario
include("stat_tests_no_mutation_pinf125.jl") # Higer infection probability scenario


text:
using CausalityTools, StatsBase

@testset "Cross mapping" begin

	x, y = rand(100), rand(100)
	@test crossmap(x, y, 3, 1) isa Float64
	@test crossmap(x, y, 3, 1, :random) isa Vector{Float64}
	@test crossmap(x, y, 3, 1, :segment) isa Vector{Float64}
end

@testset "Pairwise asymmetric inference" begin
	@test pai(x, y, 3, 1) isa Float64
	@test pai(x, y, 3, 1, :random) isa Vector{Float64}
	@test pai(x, y, 3, 1, :segment) isa Vector{Float64}
end


text:
using JuLIP, Test
using JuLIP.Potentials, JuLIP.Testing
using LinearAlgebra

calc = lennardjones(r0=rnn(:Al))
at = set_pbc!(bulk(:Al) * 2, true)
set_calculator!(at, calc)
println(@test maxnorm(forces(at)) < 1e-12)

h3("Check atoms deform correctly with the cell")
# perturb the cell shape
apply_defm!(at, I + 0.01*rand(JMatF))
# check that under this deformation the atom positions are still
# in a homogeneous lattice (i.e. they are deformed correctly with the cell)
println(@test maxnorm(forces(at)) < 1e-12)
# now perturb the atom positions as well
rattle!(at, 0.1)

# set the constraint >>> this means the deformed cell defines F0 and X0
variablecell!(at)

h3("check that energy, forces, virial, stress, dofs, gradient evaluate ... ")
energy(at)
forces(at)
gradient(at)
virial(at)
println(@test stress(at) ≈ - virial(at) / volume(at))

h3("Check that setting and getting dofs is consistent ... ")
x = dofs(at)
@assert length(x) == 3 * length(at) + 9
y = x + 0.01 * rand(Float64, size(x))
set_dofs!(at, y)
z = dofs(at)
println(@test norm(y - z, Inf) < 1e-14)

# now perform the finite-difference test
set_dofs!(at, x)
println(@test JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1))

h3("Check virial for SW potential")
si = bulk(:Si, cubic=true) * 2
rattle!(si, 0.1)
apply_defm!(at, I + 0.01*rand(JMatF))
sw = StillingerWeber()
variablecell!(si)
println(@test JuLIP.Testing.fdtest(calc, si, verbose=true, rattle=0.1))


text:
using DrWatson
@quickactivate "CDC WONDER"
DrWatson.greet()

max_stars_repo_path: CDC WONDER/scripts/intro.jl

max_stars_repo_name: Ved-Mahajan/DataExplorers

max_issues_repo_path: CDC WONDER/scripts/intro.jl

max_issues_repo_name: Ved-Mahajan/DataExplorers

max_forks_repo_path: CDC WONDER/scripts/intro.jl

max_forks_repo_name: Ved-Mahajan/DataExplorers


text:
@cdef :PyMapping_HasKeyString Cint (PyPtr, Cstring)
@cdef :PyMapping_SetItemString Cint (PyPtr, Cstring, PyPtr)
@cdef :PyMapping_GetItemString PyPtr (PyPtr, Cstring)


text:
struct BatchScalMvNormal{Tm<:AbstractMatrix,Tσ<:AbstractVector} <: ContinuousMatrixDistribution
    μ::Tm
    σ::Tσ
end

struct BatchDiagMvNormal{Tm<:AbstractMatrix,Tσ<:AbstractMatrix} <: ContinuousMatrixDistribution
    μ::Tm
    σ::Tσ
end

const BMN = Union{BatchScalMvNormal,BatchDiagMvNormal}

BatchMvNormal(μ::AbstractMatrix{T}, σ::AbstractVector{T}) where T<:Real = BatchScalMvNormal(μ,σ)
BatchMvNormal(μ::AbstractMatrix{T}, σ::AbstractMatrix{T}) where T<:Real = BatchDiagMvNormal(μ,σ)

#Base.length(d::BatchMvNormal) = size(d.μ,1)
Base.eltype(d::BMN) = eltype(d.μ)
Distributions.params(d::BMN) = (d.μ, d.σ)
Distributions.mean(d::BMN) = d.μ
Distributions.var(d::BatchDiagMvNormal) = d.σ .^2
Distributions.var(d::BatchScalMvNormal) = fillsimilar(d.σ,size(d.μ,1),1) .* reshape(d.σ .^2,1,:)

function Distributions.rand(d::BatchDiagMvNormal)
    μ, σ = d.μ, d.σ
    r = DistributionsAD.adapt_randn(Random.GLOBAL_RNG, μ, size(μ)...)
    μ .+ σ .* r
end

function Distributions.rand(d::BatchScalMvNormal)
    μ, σ = d.μ, reshape(d.σ, 1, :)
    r = DistributionsAD.adapt_randn(Random.GLOBAL_RNG, μ, size(μ)...)
    μ .+ σ .* r
end

function Distributions.logpdf(d::BMN, x::AbstractMatrix{T}) where T<:Real
    n = size(d.μ,1)
    μ = mean(d)
    σ2 = var(d)
    -(vec(sum(((x - μ).^2) ./ σ2 .+ log.(σ2), dims=1)) .+ n*log(T(2π))) ./ 2
end


text:
include("../../utils.jl")

# simulates allocation profile of some dataframes benchmarks
# by repeatedly append to a vector
function append_lots(iters=100*1024, size=1596)
    v = Float64[]
    for i = 1:iters
        append!(v,rand(size))
    end
    return v
end

@gctime append_lots()[end]

max_stars_repo_name: d-netto/GCBenchmarks

max_issues_repo_name: d-netto/GCBenchmarks

max_forks_repo_name: d-netto/GCBenchmarks


text:
using Advent
using Test

include("parse.jl")
include("site.jl")


text:
abstract type ValueFunctionLearner <: Learner end

update(l::ValueFunctionLearner) = l.update

"""
   QLearner(model, num_actions, num_demons)
"""
mutable struct QLearner{F, LU<:LearningUpdate} <: ValueFunctionLearner
    model::F
    update::LU
    num_actions::Int
    num_demons::Int
end

function LinearQLearner(update, num_features, num_actions, num_demons, w_init)
    model =  ones(num_actions*num_demons, num_features) * w_init
    QLearner(model, update, num_actions, num_demons)
end

# LinearQLearner(update, num_features, num_actions, num_demons, w_init; init=(s...)->zeros(s...)) =
#     QLearner(init(num_actions*num_demons, num_features), update, num_actions, num_demons)

predict(l::QLearner{F}, ϕ) where {F <: Matrix} = l.model*ϕ
predict(l::QLearner, ϕ) = l.model(ϕ)

(l::QLearner)(ϕ) = predict(l, ϕ)

predict(l::QLearner, ϕ::AbstractVector, a) = predict(l, ϕ)[a .+ (0:(l.num_demons-1)).*l.num_actions]
predict(l::QLearner, ϕ::AbstractMatrix, a) = predict(l, ϕ)[a .+ (0:(l.num_demons-1)).*l.num_actions, :]

(l::QLearner)(ϕ, a) = predict(l, ϕ, a)

is_linear(l::QLearner) = false
is_linear(l::QLearner{Matrix{<:Number}}) = true

function get_weights(l::QLearner)
    return l.model
end

# mutable struct VLearner <: ValueFunctionLearner
#     model::F
#     update::LU
#     num_demons::Int
# end

# predict(l::VLearner{F}, ϕ) where {F <: Matrix} = l.model*ϕ
# predict(l::VLearner, ϕ) = l.model(ϕ)

# (l::VLearner)(ϕ) = predict(l.model, ϕ)

# is_linear(l::VLearner) = false
# is_linear(l::VLearner{Matrix{<:Number}}) = true

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/525

using OrdinaryDiffEq, StaticArrays, Test

mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

function f(u,p,t) # new out-of-place definition
    SimType([-0.5*u[1] + u.f1,
    -0.5*u[2]],u.f1)
end

function f!(du,u,p,t) # old in-place definition
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

const tstop1 = [5.]
const tstop2 = [8.]

function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

function affect!_oop(integrator)
  integrator.u.f1 = 1.5
end

function affect2!_oop(integrator)
  integrator.u.f1 = 1.5
end

save_positions = (true,true)
cb = DiscreteCallback(condition, affect!, save_positions=save_positions)
save_positions = (false,true)
cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)
cbs = CallbackSet(cb,cb2)

cb_oop = DiscreteCallback(condition, affect!_oop, save_positions=save_positions)
save_positions = (false,true)
cb2_oop = DiscreteCallback(condition2, affect2!_oop, save_positions=save_positions)
cbs_oop = CallbackSet(cb_oop,cb2_oop)

u0 = SimType([10.0;10.0], 0.0)

prob_inplace = ODEProblem(f!,u0,(0.0,10.0))
prob = ODEProblem(f,u0,(0.0,10.0))

const tstop = [5.;8.]

sol = solve(prob_inplace,Tsit5(),callback = cbs, tstops=tstop)
sol = solve(prob,Tsit5(),callback = cbs_oop, tstops=tstop)

# https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/336

const A = SMatrix{2,2}([0 1;
                        0 0])

mutable struct MyStruct{T} <: DEDataVector{T}
    x::MVector{2,T}
    a::SVector{2,T}
end

function dyn(du,u,t,p)
    u.a = SVector{2}(0.0, 0.1)
    du .= A*u.x + u.a

    return nothing
end

u0   = MyStruct(MVector{2}(0.,0.), SVector{2}(0.,0.))
prob = ODEProblem(dyn, u0, (0.,10.))

@test copy(u0) isa MyStruct
@test zero(u0) isa MyStruct
@test similar(u0) isa MyStruct
@test similar(u0,Float64) isa MyStruct
@test similar(u0,Float64,size(u0)) isa MyStruct

sol  = solve(prob, Tsit5())

# https://github.com/JuliaDiffEq/StochasticDiffEq.jl/issues/247

using OrdinaryDiffEq
using StochasticDiffEq

mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

function f(du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

const tstop1 = [5.]
const tstop2 = [8.]

function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

save_positions = (true,true)

cb = DiscreteCallback(condition, affect!, save_positions=save_positions)

save_positions = (false,true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)

u0 = SimType([10.0;10.0], 0.0)
prob = ODEProblem(f,u0,(0.0,10.0))

const tstop = [5.;8.]

# here the new part

function g(du,u,p,t)
  du[1] = 1.0
  du[2] = 1.2
end

dt = 1/2^4

prob2 = SDEProblem(f,g,u0,(0.0,10.0))

# this creates an error

sol = solve(prob2,callback = cbs,tstops=tstop,EM(),dt=dt,saveat=collect(8:0.1:10))

# https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/327

struct SimulationState{T,U} <: DiffEqBase.DEDataArray{T, 1}
    x::Vector{T}
    u::U
end

#function Base.convert(::Type{SimulationState{T, U}}, a::AbstractArray{T}) where {T,U}
#    SimulationState(a, zero(eltype(a)))
#end

function open_loop(state, parameters, time)
    v = state[1] * -0.1 + state.u
    return SimulationState([v],state.u)
end

initial_conditions = SimulationState([10.0], 0.0)
time_span = (0.0, 20.0)
ode_prob = ODEProblem(open_loop, initial_conditions, time_span, nothing)

sol = solve(ode_prob, Tsit5(), reltol=1e-8, abstol=1e-8)


text:
using KSVD
using Test, LinearAlgebra

include("ksvd.jl")
include("matching_pursuit.jl")


text:
module TestUnits

using Test
using Mimi

import Mimi: verify_units, connect_param!, ComponentReference

# Try directly using verify_units
@test verify_units("kg", "kg")
@test !verify_units("kg", "MT")
@test !verify_units("kg", "")

# Create a model with some matching and some mis-matching units
@defcomp Foo begin
    output = Variable(unit="kg")
end

@defcomp Bar begin
    input = Parameter(unit="MT")
end

@defcomp Baz begin
    input = Parameter(unit="kg")
end

m = Model()
set_dimension!(m, :time, [1])
foo = add_comp!(m, Foo)
bar = add_comp!(m, Bar)
baz = add_comp!(m, Baz)

# Check that we cannot connect foo and bar...
@test_throws ErrorException bar[:input] = foo[:output]
#@test_throws ErrorException connect_param!(m, :Bar, :input,  :Foo, :output, ignoreunits=false)

# ...unless we pass ignoreunits=true
connect_param!(m, :Bar, :input,  :Foo, :output, ignoreunits=true)

# But we can connect foo and baz
baz[:input] = foo[:output]

end # module


text:
import HarmonicBalance.load

current_path = @__DIR__

@test load(current_path * "/parametron_result.jld2") isa HarmonicBalance.Result


text:
global const cmsat_lib = Sys.isunix() ? "libcryptominisat5" : "cryptominisat5"

using Libdl
try
    dlopen(cmsat_lib)
catch
    @warn "CryptoMiniSat5 library not loaded:
    please download and build from https://github.com/msoos/cryptominisat"
    rethrow()
end


text:
using Catwalk
using Test
using Random
using BenchmarkTools

Random.seed!(42)

# To run only selected tests, use e.g.:
#
#   using Pkg; Pkg.test("Catwalk", test_args=["scheduling"])
#
enabled_tests = lowercase.(ARGS)
function addtests(fname)
    key = lowercase(splitext(fname)[1])
    if isempty(enabled_tests) || key in enabled_tests
        include(fname)
    end
end

addtests("typelist.jl")
addtests("costmodel.jl")
addtests("scheduling.jl")
addtests("typesweep.jl")


text:
__precompile__(true)

module FEASTSolver

using LinearAlgebra: ldiv!, lu!, tr, dot, LU, Factorization, mul!, qr!, rmul!, lmul!, ldiv!, eigen!, svd!, norm, Diagonal, I, diagm
using IterativeSolvers: bicgstabl
using SparseArrays: similar, AbstractSparseMatrix, sprandn, sprand
using IterativeSolvers: gmres!, bicgstabl!, gmres, bicgstabl
using FastGaussQuadrature: gausslegendre
using Distributed: @distributed
using SharedArrays: SharedArray

import Base: length

export feast!, ifeast!, nlfeast!, nlfeast_opt!, nlfeast_it!, nlfeast_moments!, nlfeast_moments_SS!
export gen_feast!, dual_gen_feast!
export beyn, companion, block_SS!, nlfeast_moments_all!
export in_contour, circular_contour_trapezoidal, circular_contour_gauss, rectangular_contour_gauss, rectangular_contour_trapezoidal
export convergence_info, rational_func
export contour_estimate_eig

include("contour.jl")
include("lapack.jl")
include("utils.jl")
include("beyn.jl")
include("companion.jl")
include("feast.jl")
include("feast_experimental.jl")
include("nlfeast.jl")
include("nlfeast_experimental.jl")
include("stochastic.jl")

end # module


text:
# Some examples of plots to look at output of Gibbs sampling
using MCMCDiagnostics
using Plots
gr()

# median permutation
med_perm = median(spl_out.sigma,2)

scatter(perm_data,med_perm,legend=false,color=log.(PP))

scatter(med_perm,log.(PP),legend=false,color=perm_data)

scatter(perm_data,mean(spl_out.sigma,2),legend=false,color=log.(PP))

# calculate ess per row
perm_ess_factor = [ess_factor_estimate(log.(spl_out.sigma[i,:])) for i in 1:K  ]
perm_ess = [perm_ess_factor[i][1] for i in 1:size(perm_ess_factor,1)]
perm_lag = [perm_ess_factor[i][2] for i in 1:size(perm_ess_factor,1)]
scatter(med_perm,perm_ess,legend=false,color=log.(PP))
scatter(perm_data,perm_ess,legend=false,color=log.(PP))
scatter(perm_lag,legend=false,color=log.(PP))

scatter(PP,perm_ess,legend=false,color=log.(PP))

# calculate ess per arrival time
med_T = median(spl_out.T,2)
plot(spl_out.T,legend=false)
plot!(T_data,lw=2,color="black",legend=false)

plot(cumsum(PP[spl_out.sigma],1)[1:(end-1),:] .- spl_out.T[2:end,:], legend=false)
plot!(cumsum(PP_data)[1:(end-1)] .- T_data[2:end],color="black",lw=2)

T_ess_factor = [ess_factor_estimate(log.(spl_out.T[i,:])) for i in 1:K  ]
T_ess = [T_ess_factor[i][1] for i in 1:size(T_ess_factor,1)]
T_lag = [T_ess_factor[i][2] for i in 1:size(T_ess_factor,1)]
scatter(T_ess,legend=false,color=log.(mean(PP[spl_out.sigma],2)))
scatter(T_lag,legend=false)

plot(PP_data,lw=1.5,legend=false); plot!(median(PP[spl_out.sigma],2),color="black",lw=1.5)

plot(PP_data,lw=1.5,legend=false); plot!(mean(PP[spl_out.sigma],2),color="black",lw=1.5)


text:
module Operators

export
    Δx, Δy, ΔzF, ΔzC,
    Axᵃᵃᶜ, Axᵃᵃᶠ, Ayᵃᵃᶜ, Ayᵃᵃᶠ, Azᵃᵃᵃ,
    Vᵃᵃᶜ, Vᵃᵃᶠ,
    δxᶜᵃᵃ, δxᶠᵃᵃ, δyᵃᶜᵃ, δyᵃᶠᵃ, δzᵃᵃᶜ, δzᵃᵃᶠ,
    Ax_ψᵃᵃᶠ, Ax_ψᵃᵃᶜ, Ay_ψᵃᵃᶠ, Ay_ψᵃᵃᶜ, Az_ψᵃᵃᵃ,
    δᴶxᶜᵃᶜ, δᴶxᶜᵃᶠ, δᴶxᶠᵃᶜ, δᴶxᶠᵃᶠ, δᴶyᵃᶜᶜ, δᴶyᵃᶜᶠ, δᴶyᵃᶠᶜ, δᴶyᵃᶠᶠ, δᴶzᵃᵃᶜ, δᴶzᵃᵃᶠ,
    ∂xᶜᵃᵃ, ∂xᶠᵃᵃ, ∂yᵃᶜᵃ, ∂yᵃᶠᵃ, ∂zᵃᵃᶜ, ∂zᵃᵃᶠ,
    ∂²xᶜᵃᵃ, ∂²xᶠᵃᵃ, ∂²yᵃᶜᵃ, ∂²yᵃᶠᵃ, ∂²zᵃᵃᶜ, ∂²zᵃᵃᶠ,
    ∂⁴xᶜᵃᵃ, ∂⁴xᶠᵃᵃ, ∂⁴yᵃᶜᵃ, ∂⁴yᵃᶠᵃ, ∂⁴zᵃᵃᶜ, ∂⁴zᵃᵃᶠ,
    ℑxᶜᵃᵃ, ℑxᶠᵃᵃ, ℑyᵃᶜᵃ, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, ℑzᵃᵃᶠ,
    ℑᴶxᶜᵃᵃ, ℑᴶxᶠᵃᵃ, ℑᴶyᵃᶜᵃ, ℑᴶyᵃᶠᵃ, ℑᴶzᵃᵃᶜ, ℑᴶzᵃᵃᶠ,
    ℑxyᶜᶜᵃ, ℑxyᶠᶜᵃ, ℑxyᶠᶠᵃ, ℑxyᶜᶠᵃ, ℑxzᶜᵃᶜ, ℑxzᶠᵃᶜ, ℑxzᶠᵃᶠ, ℑxzᶜᵃᶠ, ℑyzᵃᶜᶜ, ℑyzᵃᶠᶜ, ℑyzᵃᶠᶠ, ℑyzᵃᶜᶠ,
    ℑxyzᶜᶜᶠ, ℑxyzᶜᶠᶜ, ℑxyzᶠᶜᶜ, ℑxyzᶜᶠᶠ, ℑxyzᶠᶜᶠ, ℑxyzᶠᶠᶜ,
    hdivᶜᶜᵃ, divᶜᶜᶜ, div_xyᶜᶜᵃ, div_xzᶜᵃᶜ, div_yzᵃᶜᶜ,
    ∇², ∇²hᶜᶜᵃ, ∇²hᶠᶜᵃ, ∇²hᶜᶠᵃ, ∇⁴hᶜᶜᵃ, ∇⁴hᶠᶜᵃ, ∇⁴hᶜᶠᵃ

using Oceananigans.Grids

include("areas_and_volumes.jl")
include("difference_operators.jl")
include("derivative_operators.jl")
include("interpolation_operators.jl")
include("divergence_operators.jl")
include("laplacian_operators.jl")

end


text:
# function test_incidence_to_triangles()
D = [1 -1 0 0 0;
   1 0 -1 0 0;
   0 1 -1 0 0;
   0 1 0 -1 0;
   0 1 0 0 -1;
   0 0 0 1 -1];
triangles = [1 2 3;
      2 4 5];
tri = incidence_to_triangles(D)
@test triangles == tri

D = [1 -1 0 0 0;
   1 0 -1 0 0;
   1 0 0 -1 0;
   1 0 0 0 -1;
   0 1 -1 0 0;
   0 1 0 -1 0;
   0 1 0 0 -1;
   0 0 1 -1 0;
   0 0 1 0 -1;
   0 0 0 1 -1];
triangles = [1 2 3;
      1 2 4;
      1 2 5;
      1 3 4;
      1 3 5;
      1 4 5;
      2 3 4;
      2 3 5;
      2 4 5;
      3 4 5];
tri = incidence_to_triangles(D)
@test triangles == tri

D = [1 -1 0 0 0;
   0 0 1 -1 0;
   0 0 1 0 -1;
   1 0 -1 0 0;
   1 0 0 -1 0;
   1 0 0 0 -1;
   0 1 -1 0 0;
   0 1 0 -1 0;
   0 1 0 0 -1;
   0 0 0 1 -1];
triangles = [1 2 3;
      1 2 4;
      1 2 5;
      1 3 4;
      1 3 5;
      1 4 5;
      2 3 4;
      2 3 5;
      2 4 5;
      3 4 5];
tri = incidence_to_triangles(D)
@test triangles == tri

# function test_incidence_to_dict()
D = [1 -1 0 1 0;
   0 0 1 -1 0];
node_dict = Dict(1 => Set(2), 3 => Set(4))
nd = incidence_to_dict(D)
@test node_dict == nd

D = [1 -1 0 1 0;
   1 -1 1 1 0;
   1 0 1 1 0];
node_dict = Dict(1 => Set([2, 3]))
nd = incidence_to_dict(D)
@test node_dict == nd


text:
using StaticArrays
struct Particles{T,V <: AbstractVector{SVector{3,T}}}
    npar::Int
    charge::T
    mass::T
    positions::V
    momentums::V
    self_efields::V
    self_bfields::V
end

# general interface
function Particles(;pos, mom, charge=-1.0, mass=1.0)
    return Particles(pos, mom; charge=charge, mass=mass)
end

function Particles(pos::Vector{SVector{3,T}}, mom::Vector{SVector{3,T}}; charge=-1.0, mass=1.0) where {T}
    @assert length(pos) == length(mom)
    npar = length(pos)
    return Particles(npar, charge, mass, pos, mom, fill(SVector(0.0, 0.0, 0.0), npar), fill(SVector(0.0, 0.0, 0.0), npar))
end

function Particles(pos::Matrix{T}, mom::Matrix{T}; charge=-1.0, mass=1.0) where {T}
    @assert size(pos, 1) == 3
    @assert size(mom, 1) == 3
    @assert size(pos, 2) == size(mom, 2)
    npar = size(pos, 2)
    pos_new = reinterpret(reshape, SVector{3,T}, pos)
    mom_new = reinterpret(reshape, SVector{3,T}, mom)
    return Particles(npar, charge, mass, pos_new, mom_new, reinterpret(reshape, SVector{3,T}, zeros(3, npar)), reinterpret(reshape, SVector{3,T}, zeros(3, npar)))
end

function update_self_field!(beam::Particles; lambda)
    q = beam.charge
    @inbounds Threads.@threads for i in 1:beam.npar
        beam.self_efields[i] = SVector(0.0, 0.0, 0.0)
        beam.self_bfields[i] = SVector(0.0, 0.0, 0.0)
        xi = beam.positions[i]
        for j in 1:beam.npar
            xj = beam.positions[j]
            pj = beam.momentums[j]
            xij = xi - xj
            kernel = xij / sqrt(dot(xij, xij) + dot(pj, xij)^2 + eps())^3
            beam.self_efields[i] += 2.8179403699772166e-15 * q / lambda * sqrt(1.0 + dot(pj, pj)) * kernel
            beam.self_bfields[i] += 2.8179403699772166e-15 * q / lambda * cross(pj, kernel)
        end
    end
end

function output(beam::Particles, fname)
    pos = beam.positions
    mom = beam.momentums
    open(fname, "w") do io
        writedlm(io, pos)
    end
end


text:
#push!(LOAD_PATH,"../src/")

using Documenter, ANISPROU

makedocs(sitename="ANISPROU",
         modules = [ANISPROU],
         authors = "Andrea Zunino",
         format = Documenter.HTML(prettyurls=get(ENV,"CI",nothing)=="true")  )

deploydocs(
    repo="github.com/inverseproblem/ANISPROU.jl.git",
    devbranch = "main",
    branch = "gh-pages"
)


text:
using Pkg
using Distributed

Pkg.add("Primes")
Pkg.add("Combinatorics")
Pkg.update()

# Read all files with "p{number}.jl"
files = readdir(dirname(@__FILE__()))
files = sort(filter((x)->occursin(r"^p\d*\.jl", x), files))

# Run!
totaltime = @distributed (+) for i=1:length(files)
  "Running file $i of $(length(files)) : $(files[i])" |> println
  @elapsed evalfile(files[i])
end
"Took $totaltime seconds to run $(length(files)) files" |> println

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
## rotation + dilation
N = 256; iter = 100000
θx = 2π*rand(N)
rx = sqrt.(rand(N))
x = polar2cart(rx, θx)
disc_plot(x)

# 1. Learning steps
α = 0.01
anim = @animate for n in 1:iter
    global α
    if n > .99*iter
        α -= 0.02/(.01*iter)
    end
    # update y
    global θx, rx, x
    θ = 2π*rand()
    d = (rand()-0.5)*2*log(16)
    θy = (θx .+ θ .+ (d < 0)*π) .% 2π
    ry = abs.(d) .* rx
    y = polar2cart(ry, θy)
    # find Λⱼ(n)
    nnx = Dict()
    for j in 1:N
        nnx[j] = []
    end
    for i in 1:N
        if ry[i] > 1
            continue
        end
        arr = [norm(x[j,:] - y[i,:]) for j in 1:N]
        min_val = minimum(arr)
        min_ind = findall(arr .== min_val)
        for j in min_ind
            nnx[j] = push!(nnx[j], i)
        end
    end
    for j in 1:N
        if nnx[j] != []
            x[j,:] += α * sum([y[i,:]-x[j,:] for i in nnx[j]])
        end
    end
    rx, θx = cart2polar(x)
    rx = min.(rx, 1)
    x = polar2cart(rx, θx)
end

plt = disc_plot(x)
savefig(plt, "figs/selfSimNN_2D_rotdil_100000final.png")


text:
quick = 1
# assume 2D
function quick_interpolation(field)
    afield = Array(field)
    Nx, Ny = size(field)
    @assert Nx == Ny
    N = Nx
    interpolated_field = zeros(ComplexF64, 512, 512)
    fouriedfield = fft(afield)
    if N > 512
        indices1 = 1:div(512, 2)+1
        indices2 = N-div(512, 2)+2:N
        indices = [indices1..., indices2...]
        vfield = view(fouriedfield, indices, indices)
        interpolated_field .= vfield
        return real.(ifft(interpolated_field))
    elseif N < 512
        indices1 = 1:div(N, 2)+1
        indices2 = 512-div(N, 2)+2:512
        indices = [indices1..., indices2...]
        vfield = view(interpolated_field, indices, indices)
        vfield .= fouriedfield
        return real.(ifft(interpolated_field))
    else
        return real.(afield)
    end

end

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
# # Maximum a Posteriori

using Distributions
using HMMBase
using PyPlot
using Seaborn

rc("axes", xmargin = 0) # hide
set_style("whitegrid")  # hide

# Let's consider a simple time series with one outlier:
y = rand(1000)
y[100] = 10000
figure(figsize = (9,2)) # hide
plot(y)
gcf() # hide
#-

# An MLE approach for the observations distributions parameters
# may fail with a singularity (variance = 0) if an outlier
# becomes the only observation associated to some state:
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(5,1)])
try
    fit_mle(hmm, y, display=:iter);
catch e
    println(e)
end

# We can avoid this by putting a prior on the variance:
import ConjugatePriors: InverseGamma, NormalKnownMu, posterior_canon
import StatsBase: Weights

function fit_map(::Type{<:Normal}, observations, responsibilities)
    μ = mean(observations, Weights(responsibilities))

    ss = suffstats(NormalKnownMu(μ), observations, responsibilities)
    prior = InverseGamma(2,1)
    posterior = posterior_canon(prior, ss)
    σ2 = mode(posterior)

    Normal(μ, sqrt(σ2))
end

hmm, hist = fit_mle(hmm, y, estimator=fit_map, display=:iter)
figure(figsize = (4,3)) # hide
plot(hist.logtots)
xlabel("Iteration") # hide
ylabel("Log-likelihood") # hide
gcf() # hide
#-

hmm.B

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
"""
    closure(f)
A macro that execute the input function on its output until it stays the same.
"""
macro closure(f)
    return :(function _closure(s, args...)
                while true
                    n = $f([deepcopy(s); collect(args)]...)
                    if n == s
                        return s
                    end
                    s = n
                end
            end)
end


text:
using PlanarMaps
using Test

# write your own tests here
N = NeighborCycle([3,7,9,4])

@assert length(N) == 4
@assert over(N,9,3) == 7
@assert rotate(N,9) == [9,4,3,7]
@assert NeighborCycle([1,4,2,7]) == NeighborCycle([2,7,1,4])
@assert NeighborCycle([1,4,2,7]) ≠ NeighborCycle([1,4,7,2])

@assert Face([2,1,5,4,2,3]) == Face([4, 2, 3, 2, 1, 5])

P = PlanarMap([[2,5],[3,4,1],[2],[5,2],[4,1]])
F = face(P,1,2)
@assert F == Face([1,2,4,5])

@assert cw(P,1,2) == 5
@assert cw(P,2,4) == 3
@assert ccw(P,4,5) == 2

W = UWT(100)

draw(W)
draw(W,fillfaces=true,linecolor="black")
draw(P)

# Below is the wooded triangulation from the article
# Schnyder woods, SLE(16), and Liouville Quantum Gravity
# by Li, Sun, and Watson, with the vertices numbered
# in order counterclockwise around the blue tree,
# starting with the first interior vertex and concluding
# with the three exterior vertices (10 = blue, 11 = red,
# 12 = green)
P = PlanarMap([
    [2,12,10,3],
    [12,1,3,4,5],
    [10, 11, 8, 5, 4, 2, 1],
    [5,2,3],
    [6,12,2,4,3,8,9,7],
    [12,5,7],
    [12,6,5,9,11],
    [9,5,3,11],
    [7,5,8,11],
    [1,12,11,3],
    [12,7,9,8,3,10],
    [10,1,2,5,6,7,11]])

bluetree = RootedTree(
    Dict(1=>10,2=>1,3=>10,4=>3,5=>3,6=>5,7=>5,8=>3,9=>8),
    10)

redtree = RootedTree(
    Dict(1=>3,2=>5,3=>11,4=>5,5=>9,6=>7,7=>11,8=>11,9=>11),
    11)

greentree = RootedTree(
    Dict(1=>12,2=>12,3=>2,4=>2,5=>12,6=>12,7=>12,8=>5,9=>7),
    12)

WT = WoodedTriangulation(P,bluetree,redtree,greentree)


text:
#
# A simple view datastructure for atoms of struct of array systems
#
export AtomView
struct AtomView{FS <: AbstractSystem}
    system::FS
    index::Int
end
velocity(v::AtomView)      = velocity(v.system, v.index)
position(v::AtomView)      = position(v.system, v.index)
atomic_mass(v::AtomView)   = atomic_mass(v.system, v.index)
atomic_symbol(v::AtomView) = atomic_symbol(v.system, v.index)
atomic_number(v::AtomView) = atomic_number(v.system, v.index)
element(v::AtomView)       = elements[atomic_symbol(v)]
n_dimensions(v::AtomView)  = n_dimensions(v.system)


text:
using Documenter
import Pkg
using Klimakoffer

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(Klimakoffer, :DocTestSetup, :(using Klimakoffer); recursive=true)

# Make documentation
makedocs(
    # Specify modules for which docstrings should be shown
    modules = [Klimakoffer],
    # Set sitename to Klimakoffer
    sitename="Klimakoffer.jl",
    # Provide additional formatting options
    format = Documenter.HTML(
        # Disable pretty URLs during manual testing
        prettyurls = get(ENV, "CI", nothing) == "true",
        # Explicitly add favicon as asset
        # assets = ["assets/favicon.ico"],
        # Set canonical URL to GitHub pages URL
        canonical = "https://klimakoffer.github.io/Klimakoffer.jl/dev"
    ),
    # Explicitly specify documentation structure
    pages = [
        "Home" => "index.md",
        "Reference" => "reference.md",
        "License" => "license.md"
    ],
    strict = true # to make the GitHub action fail when doctests fail, see https://github.com/neuropsychology/Psycho.jl/issues/34
)

deploydocs(
    repo = "github.com/klimakoffer/Klimakoffer.jl",
    devbranch = "main",
    push_preview = true
)


text:
using .JSort
using PyCall
import PyPlot; plt = PyPlot;
using Polynomials

#TODO Can maybe improved by fitting the actual exponential decay
function aligntime(spectra, edges; referenceindex=1,
                   prominence=50, distance=290, plot=false,
                   prompt=400, wavelength=45)
    signal = pyimport("scipy.signal")
    coefficients = Vector{Float64}[]
    aligned = Vector{eltype(spectra[1])}[]
    allpeaks = Vector{Float64}[]

    shift = first(edges)
    gain = step(edges)
    indextot = Poly([shift, gain])
    ttoindex = Poly([-shift/gain, 1/gain])
    distance = ttoindex(distance)

    for spectrum in spectra
        # Find the main peak
        peak = findpeaks(spectrum, numpeaks=1)[1]

        # Search for all smaller peaks to the right for the main peak
        region = (peak+40):length(spectrum)
        peaks, prop = signal.find_peaks(spectrum[region],
                                        prominence=prominence,
                                        distance=distance)
        peaks .+= first(region) # shift back
        peaks = [peak, peaks...]
        push!(allpeaks, peaks)
    end

    for (chn, peaks) in allpeaks |> enumerate
        # Last peak is often ugly
        N = length(peaks) - 1
        N = if N > 10 10 else N end
        peaks = peaks[1:N]
        peaks .= indextot.(peaks)
        # Main peak at 400 ms, next follows with distance 45 ms
        fact = [prompt, [prompt + wavelength*i for i in 1:N-1]...]
        @show fact
        @show peaks
        coeffs = leastsquares(peaks, fact)
        push!(coefficients, coeffs)
        if plot
            plt.scatter(peaks, repeat([chn+0.5], length(peaks)), color="r", marker="v")
        end
    end

    # flag = false
    # for (i, (peaks, spectrum)) in zip(allpeaks, spectra) |> enumerate
    #     if i == referenceindex
    #         push!(coefficients, [0.0, 1.0])
    #         push!(aligned, spectrum)
    #         continue
    #     end
    #     # Use as many peaks as possible for alignment
    #     refpeaks = allpeaks[referenceindex]
    #     N = min(length(peaks), length(refpeaks))
    #     # The last peak is often ugly. Remove it
    #     N = N-1
    #     refpeaks = refpeaks[1:N]
    #     peaks = peaks[1:N]
    #     coeffs = leastsquares(peaks, refpeaks)
    #     #push!(aligned, gainshift(spectrum, coeffs...))

    #     # Convert to time basis
    #     refpeaks .= indextot.(refpeaks)
    #     peaks .= indextot.(peaks)
    #     coeffs = leastsquares(peaks, refpeaks)
    #     push!(coefficients, coeffs)
    #     if plot
    #         plt.scatter(peaks, repeat([i+0.5], length(peaks)), color="r", marker="v")
    #         if !flag
    #             plt.scatter(refpeaks, repeat([referenceindex+0.5], length(refpeaks)), color="r", marker="v")
    #             flag = true
    #         end
    #     end
    # end
    coefficients
end

#=
# t: Uncorrected time
# e: Calibrated energy [keV]
# e_si: Calibrated SiRi back energy [keV]
=#
function chirp(time, e, e_si)
    corr[0]
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
function build_from_custom_attitude(::Type{ThrustDirectionAndAttitudeProvider}, arg0::AttitudeProvider, arg1::Vector3D)
    return jcall(ThrustDirectionAndAttitudeProvider, "buildFromCustomAttitude", ThrustDirectionAndAttitudeProvider, (AttitudeProvider, Vector3D), arg0, arg1)
end

function build_from_direction_in_frame(::Type{ThrustDirectionAndAttitudeProvider}, arg0::Frame, arg1::ThrustDirectionProvider, arg2::Vector3D)
    return jcall(ThrustDirectionAndAttitudeProvider, "buildFromDirectionInFrame", ThrustDirectionAndAttitudeProvider, (Frame, ThrustDirectionProvider, Vector3D), arg0, arg1, arg2)
end

function build_from_direction_in_lof(::Type{ThrustDirectionAndAttitudeProvider}, arg0::LOFType, arg1::ThrustDirectionProvider, arg2::Vector3D)
    return jcall(ThrustDirectionAndAttitudeProvider, "buildFromDirectionInLOF", ThrustDirectionAndAttitudeProvider, (LOFType, ThrustDirectionProvider, Vector3D), arg0, arg1, arg2)
end

function build_from_fixed_direction_in_satellite_frame(::Type{ThrustDirectionAndAttitudeProvider}, arg0::Vector3D)
    return jcall(ThrustDirectionAndAttitudeProvider, "buildFromFixedDirectionInSatelliteFrame", ThrustDirectionAndAttitudeProvider, (Vector3D,), arg0)
end

function get_attitude(obj::ThrustDirectionAndAttitudeProvider, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::ThrustDirectionAndAttitudeProvider, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_maneuver_attitude_provider(obj::ThrustDirectionAndAttitudeProvider)
    return jcall(obj, "getManeuverAttitudeProvider", AttitudeProvider, ())
end

function get_thruster_axis_in_satellite_frame(obj::ThrustDirectionAndAttitudeProvider)
    return jcall(obj, "getThrusterAxisInSatelliteFrame", Vector3D, ())
end


text:
"""
    _direct_dependencies(notebook::Notebook) -> String

Return the direct dependencies for a `notebook`.
"""
function _direct_dependencies(notebook::Notebook)::String
    deps = [last(pair) for pair in dependencies(notebook.nbpkg_ctx)]
    filter!(p -> p.is_direct_dep, deps)
    # Ignore stdlib modules.
    filter!(p -> !isnothing(p.version), deps)
    list = ["$(p.name) $(p.version)" for p in deps]
    sort!(list)
    return join(list, '\n')
end

function _direct_dependencies(notebook::AbstractString)
    return _direct_dependencies(_load_notebook(notebook))
end

"""
    _context(notebook::Notebook) -> String

Return build context, such as the Julia version and package versions, for `notebook`.
"""
function _context(notebook::Notebook)::String
    deps = _direct_dependencies(notebook)
    return """

Version

        Built with Julia $VERSION and

        $deps

        """
end


text:
using BaseTestMulti
using Base.Test


text:
JSON.lower(d::DataSource) = non_nothing_dict(d)

"""
    Table(name::String)

Table datasource.
"""
struct Table <: DataSource
    type::String
    name::String
    Table(name) = new("table", name)
end

"""
    Lookup(name::String)

Lookup datasource.
"""
struct Lookup <: DataSource
    type::String
    lookup::String
    Lookup(name) = new("lookup", name)
end

"""
    Unioned(names::Vector{String})

Union of table datasources.
"""
struct Unioned <: DataSource
    type::String
    dataSources::Vector{String}
    Unioned(names) = new("union", names)
end

"""
    Inline(column_names::Vector{String}, rows::Vector{Vector{T}})

Inline datasource.
"""
struct Inline <: DataSource
    type::String
    columnNames::Vector{String}
    rows::Vector{Vector{Any}}
    Inline(column_names, rows) = new("inline", column_names, rows)
end

"""
    QuerySource(query::Query)

Query datasource.
"""
struct QuerySource <: DataSource
    type::String
    query::Query
    QuerySource(query) = new("query", query)
end

"""
    Join(left::Union{Table, Lookup, Inline, QuerySource, Join},
        right::Union{Lookup, Inline, QuerySource}, rightPrefix::String,
        condition::String, joinType::String)

Join datasource.
"""
struct Join <: DataSource
    type::String
    left::Union{Table, Lookup, Inline, QuerySource, Join}
    right::Union{Lookup, Inline, QuerySource}
    rightPrefix::String
    condition::String
    joinType::String
    function Join(left, right, rightPrefix, condition, joinType)
        uppercase(joinType) ∈ ["INNER", "LEFT"] || error("Invalid joinType")
        new("join", left, right, rightPrefix, condition, uppercase(joinType))
    end
end


text:
function noboringzeros(n)
    while n ≠ 0 && n % 10 == 0
        n ÷= 10
    end
    return n
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using .CSTParser
using .CSTParser: EXPR, Call, Location, LocExpr, exprloc, charrange

expridx(x, ii) = (@assert isempty(ii); x)
expridx(x::Expr, ii) = isempty(ii) ? x : expridx(x.args[ii[1]], ii[2:end])
expridx(x::LocExpr, ii) = expridx(x.expr, ii)

function precedence_level(cst::EXPR, loc::Location)
  parent = cst[CSTParser.parent(loc)]
  if parent isa Union{CSTParser.BinaryOpCall,CSTParser.UnaryOpCall,CSTParser.ChainOpCall}
    Base.operator_precedence(Expr(parent).args[1])
  elseif parent isa Union{CSTParser.BinarySyntaxOpCall,CSTParser.UnarySyntaxOpCall}
    Base.operator_precedence(Expr(parent).head)
  else
    0
  end
end

# Get the range at loc, including trailing trivia from the previous node.
function separator_range(cst, loc)
  i = loc.ii[end]
  loc = CSTParser.parent(loc)
  start = charrange(cst, loc[i-1])[1][1]+1
  stop = charrange(cst, loc[i])[1][end]
  return start:stop
end

function separator(cst, loc, x::EXPR{Call}, i)
  length(x.args) == 3 && return ""
  length(x.args) == 4 && return ", "
  separator_range(cst, loc[max(i-1,4)])
end

function separator(cst, loc, x::EXPR{CSTParser.MacroCall}, i)
  if get(x.args, 2, nothing) isa CSTParser.PUNCTUATION # @foo(a, b)
    length(x.args) == 3 && return ""
    length(x.args) == 4 && return ", "
    separator_range(cst, loc[max(i-1,4)])
  else                                   # @foo a b
    length(x.args) == 1 && return " "
    outer, inner = charrange(cst, loc[i-1])
    inner[end]+1:outer[end]
  end
end

function separator(cst, loc, x::EXPR{CSTParser.TupleH}, i)
  brackets = x.args[1] isa CSTParser.PUNCTUATION
  length(x.args) == 2brackets && return ""
  length(x.args) == 1+2brackets && return ", "
  separator_range(cst, loc[max(i-1,2+brackets)])
end

ex = CSTParser.parse("a, b")
separator(ex, CSTParser.Location(), ex, 1)

ex = CSTParser.parse("(a,)")
separator(ex, CSTParser.Location(), ex, 2)

function separator(cst, loc, x::EXPR{CSTParser.Block}, i)
  out, in = charrange(cst, loc[max(i-1,1)])
  in[end]+1:out[end]
end

function separator(cst, loc, x::CSTParser.BinaryOpCall, i)
  separator_range(cst, loc[2])
end

function separator(cst::EXPR, loc::Location)
  parent = CSTParser.parent(loc)
  separator(cst, parent, cst[parent], loc.ii[end])
end

struct SourceFile
  path::String
  text::String
  cst::EXPR
  ast::LocExpr
end

function SourceFile(path::String, text = String(read(path)))
  cst = CSTParser.parse(text, true)
  SourceFile(path, text, cst, LocExpr(cst))
end

function replacement(src::SourceFile, p::Replace)
  loc = exprloc(src.ast, p.idx)
  if loc == nothing
    @warn "No location found for $(repr(CSTParser.striploc(expridx(src.ast, p.idx))))"
    return
  end
  prec = precedence_level(src.cst, loc)
  _, span = charrange(src.cst, loc)
  span => sprint(Base.show_unquoted, p.new, 0, prec)
end

function replacement(src::SourceFile, p::Insert)
  append = p.idx[end] > length(expridx(src.ast, p.idx[1:end-1]).args)
  append && (p.idx[end] -= 1)
  loc = exprloc(src.ast, p.idx)
  # TODO handle cases like this more generally
  src.cst[CSTParser.parent(loc)] isa EXPR{Call} && (loc.ii[end] = max(loc.ii[end], 2))
  _, span = charrange(src.cst, loc)
  point = append ? span[end] : span[1]-1
  sep = separator(src.cst, loc)
  sep isa AbstractRange && (sep = src.text[sep])
  (1:0).+point => sprint() do io
    append && write(io, sep)
    Base.show_unquoted(io, p.tree, 0, 0)
    append || write(io, sep)
  end
end

function replacement(src::SourceFile, p::Delete)
  loc = exprloc(src.ast, p.idx)
  span, _ = charrange(src.cst, loc)
  sep = separator(src.cst, loc)
  sep isa AbstractRange || (sep = span)
  span = span[1] > sep[1] ? (sep[1]:span[end]) : (span[1]:sep[end])
  span => ""
end

replacement(src::SourceFile, p::Patch) =
  filter(x -> x != nothing, [replacement(src, p) for p in p.ps])

function patch(io::IO, src::SourceFile, rs)
  rs = sort(rs, by=x->first(x[1]))
  offset = 1
  for r in rs
    print(io, src.text[offset:prevind(src.text, first(r[1]))])
    print(io, r[2])
    offset = last(r[1])+1
  end
  print(io, src.text[offset:end])
end

patch(io::IO, src::SourceFile, p::Patch) = patch(io, src, replacement(src, p))

patch(src::SourceFile, p) = sprint(patch, src, p)

function patch!(src::SourceFile, p)
  s = read(src.path)
  open(src.path, "w") do io
    try
      patch(io, src, p)
    catch _
      seek(io, 0)
      write(io, s)
      rethrow()
    end
  end
end

function sourcemap(f, src::SourceFile)
  expr = CSTParser.striploc(src.ast)
  ex = striplines(f(expr))
  patch(src, diff(expr, ex))
end

function sourcemap(f, path::AbstractString)
  isdir(path) && return sourcemap_dir(f, path)
  isfile(path) || error("No file at $path")
  s = SourceFile(path)
  expr = CSTParser.striploc(s.ast)
  ex = striplines(f(expr))
  patch!(s, diff(expr, ex))
  return
end

function sourcemap_dir(f, path)
  @sync for (dir, _, fs) in walkdir(path), file in fs
    @async if endswith(file, ".jl")
      @info "SourceWalk: processing $file"
      sourcemap(f, joinpath(dir, file))
    end
  end
end

sourcewalk(f, file) = sourcemap(x -> postwalk(f, x), file)

textmap(f, text) = sourcemap(f, SourceFile("", text))

textwalk(f, text) = sourcewalk(f, SourceFile("", text))


text:
#=
Model: FOM (1006 variables, 1 input)
=#
using Reachability, MAT

fom(o::Pair{Symbol, <:Any}...) = fom(Options(Dict{Symbol, Any}(o)))

function fom(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "fom.mat")
    A = float(read(file, "A")) # this matrix has Int entries

    # initial set
    X0 = Hyperrectangle(zeros(1006), [fill(0.0001, 400); zeros(606)])

    # input set
    B = read(file, "B")
    U = BallInf([0.0], 1.0)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # property: y < 185
    property = SafeStatesProperty(HalfSpace(
        read(matopen(@relpath "out.mat"), "M")[1, :], 185.))

    # =======================
    # Problem default options
    # =======================
    partition_2D = [(2*i-1:2*i) for i in 1:503] # 2D blocks
    partition_6D = vcat([(6*i-5:6*i) for i in 1:167], [1003:1006]) # 6D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [1],
                                  :partition => partition_2D,
                                  :plot_vars => [0, 1],
                                  :lazy_expm_discretize => true,
                                  :assume_sparse => true)
        # :projection_matrix => sparse(read(matopen(@relpath "out.mat"), "M")),
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => 1:1006,
                                  :partition => partition_6D,
                                  :property => property,
                                  :lazy_inputs_interval => -1,
                                  :lazy_expm_discretize => true,
                                  :assume_sparse => true)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using FileIO
path = dirname(@__FILE__)

custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

file = "$path/../reference/posterior.jld2"
data = Matrix{Float64}(load(file, "data")')
lh_expected = load(file, "lh_expected")
post_expected = load(file, "post_expected")

@testset "Check likelihood and posterior calculations" begin
    lh = likelihood(m, data)
    @test lh_expected ≈ lh

    post = DSGE.posterior(m, data)
    @test post_expected ≈ post

    x = map(α->α.value, m.parameters)
    post_at_start = DSGE.posterior!(m, x, data)
    @test post_expected ≈ post_at_start

    # Ensure if we are not evaluating at start vector, then we do not get the reference
    # posterior
    global y = x .+ 0.01
    post_not_at_start = DSGE.posterior!(m, y, data)
    ϵ = 1.0
    @test abs(post_at_start - post_not_at_start) > ϵ
end


text:
@testset "Finite difference timing" begin
    fd_input_tuple(method::Optim.FirstOrderOptimizer, prob) = ((MVP.objective(prob),),)
    fd_input_tuple(method::Optim.SecondOrderOptimizer, prob) = ((MVP.objective(prob),), (MVP.objective(prob), MVP.gradient(prob)))

    function run_optim_fd_tests(method;
                                problems = ("Extended Rosenbrock", "Large Polynomial", "Powell",
                                            "Paraboloid Diagonal", "Penalty Function I",),
                                show_name = false, show_trace = false,
                                show_time = false, show_res = false)

        # Loop over unconstrained problems
        for name in problems
            prob = MVP.UnconstrainedProblems.examples[name]
            show_name && printstyled("Problem: ", name, "\n", color=:green)
            options = Optim.Options(allow_f_increases=true, show_trace = show_trace)
            for (i, input) in enumerate(fd_input_tuple(method, prob))
                # Loop over appropriate input combinations of f, g!, and h!
                results = Optim.optimize(input..., prob.initial_x, method, options)

                debug_printing && printstyled("f-calls: $(Optim.f_calls(results))\n", color=:red)
                show_res && display(results)

                show_time && @time Optim.optimize(input..., prob.initial_x, method, options)

                @test Optim.converged(results)
                @test Optim.minimum(results) < prob.minimum + sqrt(eps(typeof(prob.minimum)))
                @test norm(Optim.minimizer(results) - prob.solutions) < 1e-2
            end
        end
    end

    @testset "Timing with LBFGS" begin
        debug_printing && printstyled("#####################\nSolver: L-BFGS\n", color=:blue)
        run_optim_fd_tests(LBFGS(), show_name=debug_printing, show_time = debug_printing)
    end

    @testset "Timing with Newton" begin
        debug_printing && printstyled("#####################\nSolver: Newton\n", color=:blue)
        run_optim_fd_tests(Newton(), show_name=debug_printing, show_time = debug_printing)
    end
end


text:
load("utils/namespace.jl")

showall(scanvars_let(:(let arg1, arg2::Int, arg3=1, arg4::Int=2, arg5=(gl=1)
    a
    b::Int
    c=1
    d::Int = 2
    local e
    local f::Int
    local g=1
    local h::Int = 2

    x, y::Float = 1, 2
    begin
        z=3
    end
    let q=1
        local f::String="f"
        r::Int = 3
    end

    k = :k
    local k::Symbol

    function f(x)
        local e
    end
    f(x,y) = (local z=x*y)

    abstract S
    abstract T <: S
    type U
        x::Uint
        y::Uint
    end
    type V <: T
        z::String
    end
end)))

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
# exing indices from a PointCloud
# http://www.pointclouds.org/documentation/tutorials/extract_indices.php

using PCLCommon
using PCLIO
using PCLFilters
using PCLSegmentation
using PCLSampleConsensus
using PCLVisualization
using Cxx

PT = PointXYZ

pcd_file = Pkg.dir("PCLIO", "test", "data", "table_scene_lms400.pcd")

cloud = PointCloud{PT}(pcd_file)
cloud_filtered = PointCloud{PT}()

vgf = VoxelGrid{PT}()
setInputCloud(vgf, cloud)
setLeafSize(vgf, 0.01, 0.01, 0.01)
filter(vgf, cloud_filtered)

println("PointCloud before filtering: ", length(cloud), " data points")
println("PointCloud after filtering: ", length(cloud_filtered), " data points")

coefficients = ModelCoefficients()
inliers = PointIndices()

seg = SACSegmentation{PT}()
setOptimizeCoefficients(seg, true)
setModelType(seg, SACMODEL_PLANE)
setMethodType(seg, SAC_RANSAC)
setMaxIterations(seg, 1000)
setDistanceThreshold(seg, 0.01)

ex = ExtractIndices{PT}()
nr_points = length(cloud_filtered)

planes = []
coefs = []
while length(cloud_filtered) > 0.3 * nr_points
    setInputCloud(seg, cloud_filtered)
    segment(seg, inliers, coefficients)

    coefvalues = icxx"$(coefficients.handle)->values;"
    indices = icxx"$(inliers.handle)->indices;"

    push!(coefs, coefvalues)

    if length(indices) == 0
        error("Could not estimate a planar model for the given dataset.")
    end

    cloud_p = PointCloud{PT}()
    cloud_f = PointCloud{PT}()

    setInputCloud(ex, cloud_filtered)
    setIndices(ex, inliers)
    setNegative(ex, false)
    filter(ex, cloud_p)
    println("PointCloud representing the planar component:",
        length(cloud_p), " data points")

    setNegative(ex, true)
    filter(ex, cloud_f)
    cloud_filtered = cloud_f

    push!(planes, cloud_p)
end

neg_plane = cloud_filtered

if isdefined(:vis) && vis
    viewer = PCLVisualizer("pcl visualizer")
    addCoordinateSystem(viewer, 1.2, 0,0,0)
    @show length(planes)
    for i in 1:length(planes)
        color = i == 1 ? (255,0,0) : i == 2 ? (0,255,0) : tuple(rand(UInt, 3)...)
        handler = PointCloudColorHandlerCustom(planes[i], color...)
        addPointCloud(viewer, planes[i], handler, id="plane $i")
    end
    white_hanlder = PointCloudColorHandlerCustom(neg_plane, 255,255,255)
    addPointCloud(viewer, neg_plane, white_hanlder, id="neg plane")
    spin(viewer)
    close(viewer)
    viewer = 0; gc(); gc();
end


text:
using Interpolations
using kRPC
using CSV
using LinearAlgebra
using DataFrames
include("../aerodynamics.jl")
conn = kRPCConnect("Aerodata")
using kRPC.Remote.SpaceCenter
totable = 1
outp_a = 1
outp_b = 1
vessel = ActiveVessel()
flight = Flight(vessel)
orbit = Orbit(vessel)
body = Body(orbit)
body_frame = ReferenceFrame(body)
control = Control(vessel)

aerodat = Aerodynamics.load_aerodata("lift_drag.csv",1.0,1.0,1.0)

function compute_reference_frame(body, lat, lon, alt::Float64)
    landing_position = SurfacePosition(body, lat, lon, body_frame)
    q_long = (
        0.0,
        sin(-lon * 0.5 * pi / 180),
        0.0,
        cos(-lon * 0.5 * pi / 180))
    q_lat = (
        0.0,
        0.0,
        sin(lat * 0.5 * pi / 180),
        cos(lat * 0.5 * pi / 180))
    cr1 = CreateRelative(body_frame, position=landing_position, rotation=q_long)
    cr2 = CreateRelative(cr1,rotation=q_lat)
    return CreateRelative(cr2,position=(alt,0.0,0.0))
end
rf = compute_reference_frame(body, -0.09720792138114814, -74.55767985150328, 0.0)
gfl = Flight(vessel, referenceFrame=rf)
sos = Float64(SpeedOfSound(gfl))

add_multiple_streams([kRPC.Remote.SpaceCenter.Delayed.AerodynamicForce(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Direction(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Velocity(gfl),
	kRPC.Remote.SpaceCenter.Delayed.Position(vessel, rf)], function (force,direction,velocity,posn)
	dir,vel = convert(Vector{Float64}, [direction...]), convert(Vector{Float64}, [velocity...])
	aerf,dragf,liftf = Aerodynamics.aero_force(aerodat,dir,vel,sos)
	aoa = acosd(dot(dir,vel)/norm(vel))
	mach = norm(vel)/sos
	liftd = cross(cross(dir, vel),vel)
	liftd = liftd/norm(liftd)
    println("lift ratio: $(liftf/dot(liftd,[force...])) ratio:$(norm(aerf)/norm(force)) dp:$(dot(aerf,force)/(norm(aerf) * norm(force)))")
end)

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
const SLICELENGTH = 10000

const MISSINGLEVEL = "."

const HEADLENGTH = 20

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
# random graph generating functions

using Distributions

function randomWalkMultiGraph(;n_edges::Int=100, alpha_prob::AbstractFloat=0.5,length_distribution::DiscreteUnivariateDistribution=Poisson(1), sizeBias::Bool=false)

  g = zeros(Int64,n_edges,2)
  g[1,:] = [1 2] # first edge
  n = 1
  nv = 2
  deg = [1;1]

  min_offset = 1 - minimum(length_distribution)

  coinDist = Bernoulli(alpha_prob)

  for i = 2:n_edges

    coin = rand(coinDist)
    vweights = (sizeBias ? deg : ones(Float64,nv) )

    stv = wsample(1:nv, vweights)  # sample start vertex

    if(Bool(coin)) # add a vertex
      nv += 1
      g[i,:] = [stv,nv]
      push!(deg,1)

    else # random walk edge
      K = rand(length_distribution) + min_offset
      edv = (K > 0 ? randomwalk(g[1:i-1,:], stv, K)[end] : stv)

      g[i,:] = [stv,edv]
      deg[stv] += 1
      deg[edv] += 1

    end

  end
  return g

end

function randomWalkSimpleGraph(;n_edges::Int=100, alpha_prob::AbstractFloat=0.5,length_distribution::DiscreteUnivariateDistribution=Poisson(1), sizeBias::Bool=false)
  # generates a simple graph from the random walk model with new vertex probability β
  # and random walk length distribution `length_distribution`
  g = zeros(Int64,n_edges,2)
  g[1,:] = [1 2] # first edge
  n = 1
  nv = 2
  deg = [1;1]

  min_offset = 1 - minimum(length_distribution)

  coinDist = Bernoulli(alpha_prob)

  for i = 2:n_edges

    coin = rand(coinDist)
    vweights = (sizeBias ? deg : ones(Float64,nv) )

    stv = wsample(1:nv, vweights)  # sample start vertex

    if(Bool(coin)) # add a vertex
      nv += 1
      g[i,:] = [stv,nv]
      push!(deg,1)

    else # random walk edge
      K = rand(length_distribution) + min_offset
      edv = (K > 0 ? randomwalk(g[1:i-1,:], stv, K)[end] : stv)

      if edv==stv || has_edge(g[1:i-1,:], stv, edv) # new vertex
        nv += 1
        g[i,:] = [stv,nv]
        push!(deg,1)

      else
        g[i,:] = [stv,edv]
        deg[stv] += 1
        deg[edv] += 1

      end

    end

  end
  return g

end

function has_edge(edgelist::Array{Int64,2},edge::Vector{Int64})
"""
  Checks for the directed edge in edgelist
"""
  return any( [ edge == edgelist[i,:] for i = 1:size(edgelist,1) ] )
end

function has_edge(edgelist::Array{Int64,2},v1::Int64,v2::Int64)
"""
  Checks for the undirected edge [v1, v2] in edgelist
"""
  return any( [ ([v1,v2] == edgelist[i,:] || [v2,v1] == edgelist[i,:]) for i = 1:size(edgelist,1) ] )
end

function randomwalk(edgelist::Array{Int64,2},start_vertex::Int64,len::Int64)

  A = edgelist2adj(edgelist)
  nv = size(A,2)
  vc = start_vertex
  p = zeros(Int64,len+1)
  p[1] = vc
  for k = 1:len
    vn = wsample(1:nv,full(A[:,vc]))
    vc = vn
    p[k+1] = vc
  end
  return p
end


text:
@doc Markdown.doc"""
    rand!(a::nf_elem, U::AbstractArray) -> nf_elem
Inplace, set the coefficients of $a$ to random elements in $U$.
$a$ is returned.
"""
function rand!(a::nf_elem, U::AbstractArray)
  for i=1:degree(parent(a))
    _num_setcoeff!(a, i-1, rand(U))
  end
  return a
end

@doc Markdown.doc"""
    rand(K::AnticNumberField, U::AbstractArray) -> nf_elem
Find an element in $K$ where the coefficients are selceted at random in $U$.
"""
function rand(K::AnticNumberField, U::AbstractArray)
  a = K()
  return rand!(a, U)
end

@doc Markdown.doc"""
    rand!(A::Generic.Mat{nf_elem}, U::AbstractArray) -> Generic.Mat{nf_elem}
Inplace, replace each element in $A$ by an element where the coefficients are
sected at random in $U$.
Returns $A$.
"""
function rand!(A::Generic.Mat{nf_elem}, U::AbstractArray)
  for i=1:nrows(A)
    for j=1:ncols(A)
      rand!(A[i,j], U)
    end
  end
  return A
end

@doc Markdown.doc"""
    rand(A::Generic.MatSpace{nf_elem}, U::AbstractArray) -> Generic.Mat{nf_elem}
Create a random matrix in $A$ where the coefficients are selected from $U$.
"""
function rand(A::Generic.MatSpace{nf_elem}, U::AbstractArray)
  return rand!(A(), U)
end

@doc Markdown.doc"""
    modular_lift(ap::Array{fq_nmod_mat, 1}, me::modular_env) -> Array
Given an array of matrices as computed by \code{modular_proj},
compute a global pre-image using some efficient CRT.
"""
function modular_lift(ap::Array{fq_nmod_mat, 1}, me::modular_env)
  A = zero_matrix(me.K, nrows(ap[1]), ncols(ap[1]))
  for i=1:nrows(A)
    for j=1:ncols(A)
      A[i,j] = modular_lift([ap[k][i,j] for k=1:length(ap)], me)
    end
  end
  return A
end

@doc Markdown.doc"""
    mod!(A::Generic.Mat{nf_elem}, m::fmpz)
Inplace: reduce all entries of $A$ modulo $m$, into the positive residue system.
"""
function mod!(A::Generic.Mat{nf_elem}, m::fmpz)
  for i=1:nrows(A)
    for j=1:ncols(A)
      _mod!(A[i,j], m)
    end
  end
end

@doc Markdown.doc"""
    mod_sym!!(A::Generic.Mat{nf_elem}, m::fmpz)
Inplace: reduce all entries of $A$ modulo $m$, into the symmetric residue system.
"""
function mod_sym!(A::Generic.Mat{nf_elem}, m::fmpz)
  for i=1:nrows(A)
    for j=1:ncols(A)
      mod_sym!(A[i,j], m)
    end
  end
end

function small_coeff(a::nf_elem, B::fmpz, i::Int)
  z = fmpz()
  Nemo.num_coeff!(z, a, i)
  return cmpabs(z, B) <= 0
end

@doc Markdown.doc"""
    rational_reconstruction(A::Generic.Mat{nf_elem}, M::fmpz) -> Bool, Generic.Mat{nf_elem}
Apply \code{rational_reconstruction} to each entry of $M$.
"""
function rational_reconstruction2(A::Generic.Mat{nf_elem}, M::fmpz)
  B = similar(A)
  sM = root(M, 2)
  d = one(A[1,1])
  di = inv(d)
  for i=1:nrows(A)
    for j=1:ncols(A)
      a = A[i,j]*d
      mod_sym!(a, M)
      if all(i->small_coeff(a, sM, i), 1:a.elem_length)
        B[i,j] = a*di
      else
        n, dn = algebraic_reconstruction(a, M)
        d*=dn
        if any(i->!small_coeff(d, sM, i), 1:a.elem_length)
          println("early $i $j abort")
          return false, B
        end
        di*=inv(dn)
        B[i,j] = n*di
      end
    end
  end
  println("final den: $d")
  return true, B//d
end

function rational_reconstruction(A::Generic.Mat{nf_elem}, M::fmpz)
  B = similar(A)
  for i=1:nrows(A)
    for j=1:ncols(A)
      fl, B[i,j] = rational_reconstruction(A[i,j], M)
      if !fl
        return false, B
      end
    end
  end
  return true, B
end

function algebraic_reconstruction(a::nf_elem, M::fmpz)
  K = parent(a)
  n = degree(K)
  Znn = MatrixSpace(FlintZZ, n, n, false)
#  L = [ Znn(1) representation_matrix_q(a)[1] ; Znn(0) Znn(M)]
  L = vcat(hcat(Znn(1), representation_matrix_q(a)[1]), hcat(Znn(0),Znn(M)))
  lll!(L)
  d = Nemo.elem_from_mat_row(K, sub(L, 1:1, 1:n), 1, fmpz(1))
  n = Nemo.elem_from_mat_row(K, sub(L, 1:1, n+1:2*n), 1, fmpz(1))
  return n,d
  return true, n//d
end

function algebraic_reconstruction(a::nf_elem, M::NfAbsOrdIdl)
  K = parent(a)
  n = degree(K)
  Znn = MatrixSpace(FlintZZ, n, n, false)
  L = [ Znn(1) representation_matrix_q(a)[1] ; Znn(0) basis_matrix(M, copy = false)]
  lll!(L)
  d = Nemo.elem_from_mat_row(K, sub(L, 1:1, 1:n), 1, fmpz(1))
  n = Nemo.elem_from_mat_row(K, sub(L, 1:1, n+1:2*n), 1, fmpz(1))
  return n,d
  return true, n//d
end

@doc Markdown.doc"""
    divexact!(A::Generic.Mat{nf_elem}, p::fmpz)
Inplace: divide each entry by $p$.
"""
function divexact!(A::Generic.Mat{nf_elem}, p::fmpz)
  for i=1:nrows(A)
    for j=1:ncols(A)
      A[i,j] = A[i,j]//p
    end
  end
end

#TODO/ To experiment:
# - vector reconstruction ala Storjohan
# - reconstruction with algebraic denominators
# - vector reconstruction with algebraic denominators
# - Talk to Bill: fq_nmod_mat is missing in Nemo, the inv is dreadfully slow
# - extend to non-unique solutions
# - make Aip*D mult faster, A*y as well?
#
function solve_dixon(A::Generic.Mat{nf_elem}, B::Generic.Mat{nf_elem})
  p = next_prime(p_start)
  K = base_ring(A)

  me = modular_init(K, p)
  ap = modular_proj(A, me)
  @time ap = [inv(x) for x= ap]
  Aip = modular_lift(ap, me)
  sol = 0*B
  D = B
  pp = fmpz(1)
  last_SOL = false
  nd = 0
  while true
    nd += 1
    y = Aip*D
    mod!(y, fmpz(p))
    sol += y*pp
    pp *= p
    fl, SOL = rational_reconstruction(sol, pp)
#    t = A*sol-B
#    mod!(t, pp)
#    @assert iszero(t)
    if fl
      if last_SOL == SOL && A*SOL == B
        println("used $nd $p-adic digits")
        return SOL
      else
        last_SOL = SOL
      end
    end
    D = D - A*y
    divexact!(D, fmpz(p))
#    if nbits(pp) > 10000 # a safety device to avoid infinite loops
#      error("not work")
#    end
  end
end


text:
############################################################################################
# Define abstract types
"""
$(TYPEDEF)
Available Particle Filter Kernels.
"""
abstract type ParticleKernel <: AbstractKernel end

"""
$(TYPEDEF)
Super type for various particle referencing techniques.
"""
abstract type ParticleReferencing end

"""
$(TYPEDEF)
Super type for various particle referencing techniques.
"""
abstract type AbstractParticles end

############################################################################################
# Import sub-folder
include("Core/Core.jl")
include("Kernels/Kernels.jl")

############################################################################################
# Export
export update!,
    init,
    init!,
    propose,
    propose!,
    propagate!,
    predict,
    ParticleKernel,
    ParticleReferencing,
    AbstractParticles


text:
using JLD2
@testset "shkdcmp" begin
    m = E7A.model
    expected = Workspace(load("data/shkdcmp_E7A.jld2"))

    empty!(m.sstate.constraints)
    @steadystate m lc = 1
    @steadystate m linv = 1
    clear_sstate!(m)
    sssolve!(m)
    @test 0 == check_sstate(m)

    p = Plan(m, 2021Q1:2035Q4)
    exog_data = steadystatedata(m, p)
    exog_data[2021Q1, m.shocks] .+= 0.1
    result = shockdecomp(m, p, exog_data; fctype = fcnatural)
    @test @compare result expected quiet
end


text:
module TestPushfirst

include("preamble.jl")

@testset begin
    @test pushfirst!!([0.0], 1.0) == [1.0, 0.0]
    @test pushfirst!!([0], 1.0) == [1.0, 0.0]
    @test pushfirst!!((0,), 1) === (1, 0)
    @test pushfirst!!(SVector(0), 1) == SVector(1, 0)
    @test_broken pushfirst!!(SVector(0), 1) === SVector(1, 0)
end

end  # module


text:
function power_posterior(cm::ConditionalModel)
    m = cm.model

    pr = prior(m, observed(cm)...)
    lik = likelihood(m, observed(cm)...)

    avs = argvals(cm)
    obs = observations(cm)
    function f(t, x)
        logdensity(pr(avs), x) + t * logdensity(lik(merge(avs, x)) | obs)
    end

    return f
end

# EXAMPLE
# m = @model begin
#     x ~ Normal()
#     y ~ Normal(x,1) |> iid(3)
#     return y
# end;

# prior(m, :y)

# likelihood(m, :y)

# y = rand(m())

# f = power_posterior(m() | (;y))

# f(1.0, (x=0.2,))
# f(0.0, (x=0.2,))


text:
using Documenter, AWSTools, FilePathsBase

makedocs(
    modules=[AWSTools],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets=[
            "assets/invenia.css",
        ],
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaCloud/AWSTools.jl/blob/{commit}{path}#L{line}",
    sitename="AWSTools.jl",
    authors="Invenia Technical Computing",
    checkdocs = :exports,
    strict = true,
)