text:
#=
https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/
=#

#
# libbson
#

function bson_oid_init(oid_ref::Ref{BSONObjectId}, context::Ptr{Cvoid})
    ccall((:bson_oid_init, libbson), Cvoid, (Ref{BSONObjectId}, Ptr{Cvoid}), oid_ref, context)
end

function bson_oid_init_from_string(oid_ref::Ref{BSONObjectId}, oid_string::String)
    ccall((:bson_oid_init_from_string, libbson), Cvoid, (Ref{BSONObjectId}, Cstring), oid_ref, oid_string)
end

function bson_oid_to_string(oid::BSONObjectId)
    buffer_len = 25
    buffer = zeros(UInt8, buffer_len)
    ccall((:bson_oid_to_string, libbson), Cvoid, (Ref{BSONObjectId}, Ref{UInt8}), Ref(oid), Ref(buffer, 1))
    @assert buffer[end] == 0
    return String(buffer[1:end-1])
end

function bson_oid_compare(oid1::BSONObjectId, oid2::BSONObjectId)
    ccall((:bson_oid_compare, libbson), Cint, (Ref{BSONObjectId}, Ref{BSONObjectId}), Ref(oid1), Ref(oid2))
end

function bson_oid_get_time_t(oid::BSONObjectId)
    ccall((:bson_oid_get_time_t, libbson), Clong, (Ref{BSONObjectId},), Ref(oid))
end

function bson_oid_is_valid(str::String)
    str_length = Csize_t(length(str))
    ccall((:bson_oid_is_valid, libbson), Bool, (Cstring, Csize_t), str, str_length)
end

function bson_append_oid(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::BSONObjectId)
    ccall((:bson_append_oid, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Ref{BSONObjectId}), bson_document, key, key_length, value)
end

function bson_append_int32(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Int32)
    ccall((:bson_append_int32, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Cint), bson_document, key, key_length, value)
end

function bson_append_int64(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Int64)
    ccall((:bson_append_int64, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Clonglong), bson_document, key, key_length, value)
end

function bson_append_utf8(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::String, len::Int)
    ccall((:bson_append_utf8, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Cstring, Cint), bson_document, key, key_length, value, len)
end

function bson_append_bool(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Bool)
    ccall((:bson_append_bool, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Bool), bson_document, key, key_length, value)
end

function bson_append_double(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Float64)
    ccall((:bson_append_double, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Cdouble), bson_document, key, key_length, value)
end

function bson_append_date_time(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Int64)
    ccall((:bson_append_date_time, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Clonglong), bson_document, key, key_length, value)
end

function bson_append_document(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Ptr{Cvoid})
    ccall((:bson_append_document, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}), bson_document, key, key_length, value)
end

function bson_append_array(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::Ptr{Cvoid})
    ccall((:bson_append_array, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}), bson_document, key, key_length, value)
end

function bson_append_binary(bson_document::Ptr{Cvoid}, key::String, key_length::Int, subtype::BSONSubType, value::Vector{UInt8}, val_length::UInt32)
    ccall((:bson_append_binary, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, BSONSubType, Ptr{Cvoid}, Culong), bson_document, key, key_length, subtype, value, val_length)
end

function bson_append_code(bson_document::Ptr{Cvoid}, key::String, key_length::Int, value::String)
    ccall((:bson_append_code, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint, Cstring), bson_document, key, key_length, value)
end

#bool bson_append_null (bson_t *bson, const char *key, int key_length);
function bson_append_null(bson_document::Ptr{Cvoid}, key::String, key_length::Int)
    ccall((:bson_append_null, libbson), Bool, (Ptr{Cvoid}, Cstring, Cint), bson_document, key, key_length)
end

function bson_new()
    ccall((:bson_new, libbson), Ptr{Cvoid}, ())
end

function bson_new_from_json(data::String, len::Int=-1)
    ccall((:bson_new_from_json, libbson), Ptr{Cvoid}, (Ptr{UInt8}, Cssize_t, Ptr{Cvoid}), data, len, C_NULL)
end

function bson_destroy(bson_document::Ptr{Cvoid})
    ccall((:bson_destroy, libbson), Cvoid, (Ptr{Cvoid},), bson_document)
end

function bson_as_canonical_extended_json(bson_document::Ptr{Cvoid})
    ccall((:bson_as_canonical_extended_json, libbson), Cstring, (Ptr{Cvoid}, Ptr{Cvoid}), bson_document, C_NULL)
end

function bson_as_relaxed_extended_json(bson_document::Ptr{Cvoid})
    ccall((:bson_as_relaxed_extended_json, libbson), Cstring, (Ptr{Cvoid}, Ptr{Cvoid}), bson_document, C_NULL)
end

function bson_copy(bson_document::Ptr{Cvoid})
    ccall((:bson_copy, libbson), Ptr{Cvoid}, (Ptr{Cvoid},), bson_document)
end

function bson_has_field(bson_document::Ptr{Cvoid}, key::String)
    ccall((:bson_has_field, libbson), Bool, (Ptr{Cvoid}, Cstring), bson_document, key)
end

function bson_iter_init(iter_ref::Ref{BSONIter}, bson_document::Ptr{Cvoid})
    ccall((:bson_iter_init, libbson), Bool, (Ref{BSONIter}, Ptr{Cvoid}), iter_ref, bson_document)
end

function bson_iter_next(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_next, libbson), Bool, (Ref{BSONIter},), iter_ref)
end

function bson_iter_find(iter_ref::Ref{BSONIter}, key::String)
    ccall((:bson_iter_find, libbson), Bool, (Ref{BSONIter}, Cstring), iter_ref, key)
end

function bson_iter_init_find(iter_ref::Ref{BSONIter}, bson_document::Ptr{Cvoid}, key::String)
    ccall((:bson_iter_init_find, libbson), Bool, (Ref{BSONIter}, Ptr{Cvoid}, Cstring), iter_ref, bson_document, key)
end

function bson_iter_type(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_type, libbson), BSONType, (Ref{BSONIter},), iter_ref)
end

function bson_iter_key(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_key, libbson), Cstring, (Ref{BSONIter},), iter_ref)
end

function bson_iter_int32(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_int32, libbson), Int32, (Ref{BSONIter},), iter_ref)
end

function bson_iter_int64(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_int64, libbson), Int64, (Ref{BSONIter},), iter_ref)
end

function bson_iter_utf8(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_utf8, libbson), Cstring, (Ref{BSONIter}, Ptr{Cvoid}), iter_ref, C_NULL)
end

function bson_iter_bool(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_bool, libbson), Bool, (Ref{BSONIter},), iter_ref)
end

function bson_iter_double(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_double, libbson), Float64, (Ref{BSONIter},), iter_ref)
end

function bson_iter_oid(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_oid, libbson), Ptr{BSONObjectId}, (Ref{BSONIter},), iter_ref)
end

function bson_iter_date_time(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_date_time, libbson), Clonglong, (Ref{BSONIter},), iter_ref)
end

function bson_iter_recurse(iter_ref::Ref{BSONIter}, child_iter_ref::Ref{BSONIter})
    ccall((:bson_iter_recurse, libbson), Bool, (Ref{BSONIter}, Ref{BSONIter}), iter_ref, child_iter_ref)
end

function bson_iter_code(iter_ref::Ref{BSONIter})
    ccall((:bson_iter_code, libbson), Cstring, (Ref{BSONIter}, Ptr{UInt32}), iter_ref, C_NULL)
end

function bson_iter_binary(iter_ref::Ref{BSONIter}, lengthPtr::Array{UInt32}, dataPtr::Array{Ptr{UInt8}})
    bsonsubtype = BSON_SUBTYPE_BINARY
    ccall((:bson_iter_binary, libbson), Cvoid, (Ref{BSONIter}, Ref{BSONSubType}, Ptr{UInt32}, Ptr{Ptr{UInt8}}),iter_ref, bsonsubtype, lengthPtr, dataPtr)
    nothing
end

function bson_free(mem::Ptr{Cvoid})
    ccall((:bson_free, libbson), Cvoid, (Ptr{Cvoid},), mem)
end

function bson_reader_destroy(bson_reader_handle::Ptr{Cvoid})
    ccall((:bson_reader_destroy, libbson), Cvoid, (Ptr{Cvoid},), bson_reader_handle)
end

function bson_reader_new_from_file(filepath::AbstractString, bson_error::BSONError)
    ccall((:bson_reader_new_from_file, libbson), Ptr{Cvoid}, (Cstring, Ref{BSONError}), filepath, Ref(bson_error))
end

function bson_reader_new_from_data(data::Ptr{UInt8}, data_length::Integer)
    ccall((:bson_reader_new_from_data, libbson), Ptr{Cvoid}, (Ptr{UInt8}, Csize_t), data, data_length)
end

function bson_reader_read(bson_reader_handle::Ptr{Cvoid}, reached_eof_ref::Ref{Bool})
    ccall((:bson_reader_read, libbson), Ptr{Cvoid}, (Ptr{Cvoid}, Ref{Bool}), bson_reader_handle, reached_eof_ref)
end

function bson_writer_destroy(bson_writer_handle::Ptr{Cvoid})
    ccall((:bson_writer_destroy, libbson), Cvoid, (Ptr{Cvoid},), bson_writer_handle)
end

function bson_writer_begin(bson_writer_handle::Ptr{Cvoid}, bson_document_handle_ref::Ref{Ptr{Cvoid}})
    ccall((:bson_writer_begin, libbson), Bool, (Ptr{Cvoid}, Ref{Ptr{Cvoid}}), bson_writer_handle, bson_document_handle_ref)
end

function bson_writer_end(bson_writer_handle::Ptr{Cvoid})
    ccall((:bson_writer_end, libbson), Cvoid, (Ptr{Cvoid},), bson_writer_handle)
end

function bson_writer_new(buffer_handle_ref::Ref{Ptr{UInt8}}, buffer_length_ref::Ref{Csize_t}, offset::Csize_t, realloc_func::Ptr{Cvoid}, realloc_func_ctx::Ptr{Cvoid})
    ccall((:bson_writer_new, libbson), Ptr{Cvoid}, (Ref{Ptr{UInt8}}, Ref{Csize_t}, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}), buffer_handle_ref, buffer_length_ref, offset, realloc_func, realloc_func_ctx)
end

function bson_writer_get_length(bson_writer_handle::Ptr{Cvoid})
    ccall((:bson_writer_get_length, libbson), Csize_t, (Ptr{Cvoid},), bson_writer_handle)
end

function bson_copy_to(src_bson_handle::Ptr{Cvoid}, dst_bson_handle::Ptr{Cvoid})
    ccall((:bson_copy_to, libbson), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), src_bson_handle, dst_bson_handle)
end

function bson_copy_to_excluding_noinit(src_bson_handle::Ptr{Cvoid}, dst_bson_handle::Ptr{Cvoid})

    # this hack will create a key that is not present in src_bson
    # since bson_copy_to_excluding_noinit requires at least one `exclude` arg.
    function exclude_key(bson_handle::Ptr{Cvoid})
        new_exclude_key() = "___" * string(Int(rand(UInt16)))

        exclude = new_exclude_key()
        bson = BSON(bson_handle, enable_finalizer=false) # disable finalizer

        while haskey(bson, exclude)
            exclude = new_exclude_key()
        end

        return exclude
    end

    exclude = exclude_key(src_bson_handle)

    ccall((:bson_copy_to_excluding_noinit, libbson), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring), src_bson_handle, dst_bson_handle, exclude, C_NULL)
end

function bson_copy_to_excluding_noinit(src_bson_handle::Ptr{Cvoid}, dst_bson_handle::Ptr{Cvoid}, exclude::AbstractString)
    ccall((:bson_copy_to_excluding_noinit, libbson), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring), src_bson_handle, dst_bson_handle, exclude, C_NULL)
end

#
# libmongoc
#

function mongoc_init()
    ccall((:mongoc_init, libmongoc), Cvoid, ())
end

function mongoc_uri_new_with_error(uri_string::String, bson_error::BSONError)
    ccall((:mongoc_uri_new_with_error, libmongoc), Ptr{Cvoid}, (Cstring, Ref{BSONError}), uri_string, Ref(bson_error))
end

function mongoc_uri_destroy(uri_handle::Ptr{Cvoid})
    ccall((:mongoc_uri_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), uri_handle)
end

function mongoc_client_new_from_uri(uri_handle::Ptr{Cvoid})
    ccall((:mongoc_client_new_from_uri, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid},), uri_handle)
end

function mongoc_client_destroy(client_handle::Ptr{Cvoid})
    ccall((:mongoc_client_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), client_handle)
end

function mongoc_client_set_appname(client_handle::Ptr{Cvoid}, appname::String)
    ccall((:mongoc_client_set_appname, libmongoc), Bool, (Ptr{Cvoid}, Cstring), client_handle, appname)
end

function mongoc_client_get_database(client_handle::Ptr{Cvoid}, db_name::String)
    ccall((:mongoc_client_get_database, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Cstring), client_handle, db_name)
end

function mongoc_client_find_databases_with_opts(client_handle::Ptr{Cvoid}, bson_opts::Ptr{Cvoid})
    ccall((:mongoc_client_find_databases_with_opts, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), client_handle, bson_opts)
end

function mongoc_client_start_session(client_handle::Ptr{Cvoid}, session_options_handle::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_client_start_session, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), client_handle, session_options_handle, Ref(bson_error))
end

function mongoc_client_session_destroy(session_handle::Ptr{Cvoid})
    ccall((:mongoc_client_session_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), session_handle)
end

function mongoc_client_session_append(session_handle::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_client_session_append, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), session_handle, bson_opts, Ref(bson_error))
end

function mongoc_client_session_start_transaction(session_handle::Ptr{Cvoid}, transaction_options_handle::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_client_session_start_transaction, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), session_handle, transaction_options_handle, Ref(bson_error))
end

function mongoc_client_session_abort_transaction(session_handle::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_client_session_abort_transaction, libmongoc), Bool, (Ptr{Cvoid}, Ref{BSONError}), session_handle, Ref(bson_error))
end

function mongoc_client_session_commit_transaction(session_handle::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_client_session_commit_transaction, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), session_handle, bson_reply, Ref(bson_error))
end

function mongoc_client_session_in_transaction(session_handle::Ptr{Cvoid})
    ccall((:mongoc_client_session_in_transaction, libmongoc), Bool, (Ptr{Cvoid},), session_handle)
end

function mongoc_session_opts_new()
    ccall((:mongoc_session_opts_new, libmongoc), Ptr{Cvoid}, ())
end

function mongoc_session_opts_destroy(session_options_handle::Ptr{Cvoid})
    ccall((:mongoc_session_opts_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), session_options_handle)
end

function mongoc_session_opts_get_causal_consistency(session_options_handle::Ptr{Cvoid})
    ccall((:mongoc_session_opts_get_causal_consistency, libmongoc), Bool, (Ptr{Cvoid},), session_options_handle)
end

function mongoc_session_opts_set_causal_consistency(session_options_handle::Ptr{Cvoid}, casual_consistency::Bool)
    ccall((:mongoc_session_opts_set_causal_consistency, libmongoc), Cvoid, (Ptr{Cvoid}, Bool), session_options_handle, casual_consistency)
end

function mongoc_database_destroy(database_handle::Ptr{Cvoid})
    ccall((:mongoc_database_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), database_handle)
end

function mongoc_database_get_collection(database_handle::Ptr{Cvoid}, collection_name::String)
    ccall((:mongoc_database_get_collection, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Cstring), database_handle, collection_name)
end

function mongoc_database_find_collections_with_opts(database_handle::Ptr{Cvoid}, bson_opts::Ptr{Cvoid})
    ccall((:mongoc_database_find_collections_with_opts, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), database_handle, bson_opts)
end

function mongoc_database_command_simple(database_handle::Ptr{Cvoid}, bson_command::Ptr{Cvoid}, read_prefs::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_database_command_simple, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), database_handle, bson_command, read_prefs, bson_reply, Ref(bson_error))
end

function mongoc_database_add_user(database_handle::Ptr{Cvoid}, username::String, password::String, bson_roles::Ptr{Cvoid}, bson_custom_data::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_database_add_user, libmongoc), Bool, (Ptr{Cvoid}, Cstring, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), database_handle, username, password, bson_roles, bson_custom_data, Ref(bson_error))
end

function mongoc_database_remove_user(database_handle::Ptr{Cvoid}, username::String, bson_error::BSONError)
    ccall((:mongoc_database_remove_user, libmongoc), Bool, (Ptr{Cvoid}, Cstring, Ref{BSONError}), database_handle, username, Ref(bson_error))
end

function mongoc_collection_command_simple(collection_handle::Ptr{Cvoid}, bson_command::Ptr{Cvoid}, read_prefs::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_command_simple, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_command, read_prefs, bson_reply, Ref(bson_error))
end

function mongoc_collection_destroy(collection_handle::Ptr{Cvoid})
    ccall((:mongoc_collection_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), collection_handle)
end

function mongoc_collection_insert_one(collection_handle::Ptr{Cvoid}, bson_document::Ptr{Cvoid}, bson_options::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_insert_one, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_document, bson_options, bson_reply, Ref(bson_error))
end

function mongoc_collection_find_with_opts(collection_handle::Ptr{Cvoid}, bson_filter::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, read_prefs::Ptr{Cvoid})
    ccall((:mongoc_collection_find_with_opts, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), collection_handle, bson_filter, bson_opts, read_prefs)
end

function mongoc_collection_count_documents(collection_handle::Ptr{Cvoid}, bson_filter::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, read_prefs::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_count_documents, libmongoc), Int64, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_filter, bson_opts, read_prefs, bson_reply, Ref(bson_error))
end

function mongoc_collection_create_bulk_operation_with_opts(collection_handle::Ptr{Cvoid}, bson_opts::Ptr{Cvoid})
    ccall((:mongoc_collection_create_bulk_operation_with_opts, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), collection_handle, bson_opts)
end

function mongoc_collection_delete_one(collection_handle::Ptr{Cvoid}, bson_selector::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_delete_one, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_selector, bson_opts, bson_reply, Ref(bson_error))
end

function mongoc_collection_delete_many(collection_handle::Ptr{Cvoid}, bson_selector::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_delete_many, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_selector, bson_opts, bson_reply, Ref(bson_error))
end

function mongoc_collection_update_one(collection_handle::Ptr{Cvoid}, bson_selector::Ptr{Cvoid}, bson_update::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_update_one, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_selector, bson_update, bson_opts, bson_reply, Ref(bson_error))
end

function mongoc_collection_update_many(collection_handle::Ptr{Cvoid}, bson_selector::Ptr{Cvoid}, bson_update::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_collection_update_many, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), collection_handle, bson_selector, bson_update, bson_opts, bson_reply, Ref(bson_error))
end

function mongoc_collection_aggregate(collection_handle::Ptr{Cvoid}, flags::QueryFlags, bson_pipeline::Ptr{Cvoid}, bson_opts::Ptr{Cvoid}, read_prefs::Ptr{Cvoid})
    ccall((:mongoc_collection_aggregate, libmongoc), Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), collection_handle, flags, bson_pipeline, bson_opts, read_prefs)
end

function mongoc_cursor_destroy(cursor_handle::Ptr{Cvoid})
    ccall((:mongoc_cursor_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), cursor_handle)
end

function mongoc_cursor_next(cursor_handle::Ptr{Cvoid}, bson_document_ref::Ref{Ptr{Cvoid}})
    ccall((:mongoc_cursor_next, libmongoc), Bool, (Ptr{Cvoid}, Ref{Ptr{Cvoid}}), cursor_handle, bson_document_ref)
end

function mongoc_cursor_set_limit(cursor_handle::Ptr{Cvoid}, limit::Int)
    ccall((:mongoc_cursor_set_limit, libmongoc), Bool, (Ptr{Cvoid}, Int64), cursor_handle, limit)
end

function mongoc_cursor_error(cursor_handle::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_cursor_error, libmongoc), Bool, (Ptr{Cvoid}, Ref{BSONError}), cursor_handle, Ref(bson_error))
end

function mongoc_bulk_operation_destroy(bulk_operation_handle::Ptr{Cvoid})
    ccall((:mongoc_bulk_operation_destroy, libmongoc), Cvoid, (Ptr{Cvoid},), bulk_operation_handle)
end

function mongoc_bulk_operation_insert_with_opts(bulk_operation_handle::Ptr{Cvoid}, bson_document::Ptr{Cvoid}, bson_options::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_bulk_operation_insert_with_opts, libmongoc), Bool, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), bulk_operation_handle, bson_document, bson_options, Ref(bson_error))
end

function mongoc_bulk_operation_execute(bulk_operation_handle::Ptr{Cvoid}, bson_reply::Ptr{Cvoid}, bson_error::BSONError)
    ccall((:mongoc_bulk_operation_execute, libmongoc), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{BSONError}), bulk_operation_handle, bson_reply, Ref(bson_error))
end


text:
export DualCGPEvolution, DualCGPGAEvo, NSGA2Evo, generation

import Cambrian.populate, Cambrian.evaluate, Cambrian.log_gen, Cambrian.save_gen
using Statistics
using Formatting

mutable struct NSGA2Evo{T} <: Cambrian.AbstractEvolution
    config::NamedTuple
    logid::String
    logger::CambrianLogger
	resdir::String
    population::Vector{T}
    fitness::Function
    gen::Int64
    atari_games::Vector{Game}
end

function NSGA2Evo(
    config::NamedTuple,
    resdir::String,
    fitness::Function,
    init_population::Function,
    rom_name::String
)
    logid = config.id
    log_path = joinpath(resdir, logid, "logs/logs.csv")
    logger = CambrianLogger(log_path)
    population = init_population(config)
    atari_games = [Game(rom_name, 0) for _ in 1:length(population)]
    NSGA2Evo(config, logid, logger, resdir, population, fitness, 0, atari_games)
end

populate(e::NSGA2Evo{T}) where T = IICGP.nsga2_populate(e)
evaluate(e::NSGA2Evo{T}) where T = IICGP.fitness_evaluate(e, e.fitness)

function get_n_elites(population::Vector{T}) where T
	count(ind -> ind.is_elite, population)
end

function log_gen(
    e::NSGA2Evo{T},
    max_fitness::Vector{Float64},
	min_fitness::Vector{Float64}
) where T
	sep = ";"
    if e.gen == 1
        #=
		with_logger(e.logger) do
            @info Formatting.format(
                string("generation", sep, "rank", sep, "fitness", sep,
					   "normalized_fitness", sep, "reached_frames", sep,
					   "dna_id")
            )
        end
		=#
		f = open(joinpath(e.resdir, e.logid, "logs/logs.csv"), "w+")
        write(f, string("gen_number", sep, "rank", sep, "fitness", sep,
			"normalized_fitness", sep, "reached_frames", sep, "dna_id\n"))
        close(f)
    end
	enco_path = joinpath(e.resdir, e.logid, Formatting.format("gens/encoder_{1:04d}", e.gen))
    cont_path = joinpath(e.resdir, e.logid, Formatting.format("gens/controller_{1:04d}", e.gen))
    mkpath(enco_path)
    mkpath(cont_path)
    for i in eachindex(e.population)
        dna_id = Formatting.format("{1:04d}", i)
		# Log results
		raw_fitness = e.population[i].fitness .* (max_fitness .- min_fitness) .+ min_fitness
        #=
		with_logger(e.logger) do
            @info Formatting.format(
                string(e.gen, sep, e.population[i].rank, sep, raw_fitness, sep,
					   e.population[i].fitness, sep,
					   e.population[i].reached_frames, sep, dna_id)
            )
        end
		=#
		f = open(joinpath(e.resdir, e.logid, "logs/logs.csv"), "a+")
        write(f, Formatting.format(string(e.gen, sep, e.population[i].rank, sep,
			raw_fitness, sep, e.population[i].fitness, sep,
			e.population[i].reached_frames, sep, dna_id, "\n")))
        close(f)
		# Log individuals
		enco = IPCGPInd(e.config.e_config, e.population[i].e_chromosome)
		enco.fitness .= e.population[i].fitness
		f = open(string(enco_path, "/", dna_id, ".dna"), "w+")
        write(f, string(enco))
        close(f)
        cont = CGPInd(e.config.c_config, e.population[i].c_chromosome)
		cont.fitness .= e.population[i].fitness
        f = open(string(cont_path, "/", dna_id, ".dna"), "w+")
        write(f, string(cont))
        close(f)
    end
    flush(e.logger.stream)
end

function generation(
	e::NSGA2Evo{T},
    max_fitness::Vector{Float64},
	min_fitness::Vector{Float64}
) where T
    #==#
	# Sort all individuals according to Pareto efficiency
    fast_non_dominated_sort!(e)
    # Select elites as most Pareto efficient solutions
    sort!(e.population, by = ind -> ind.rank)
    current_rank = 1
    i_start = 1
    i_end = findlast(ind -> ind.rank == current_rank, e.population)
    while get_n_elites(e.population) < e.config.n_elite
        n_candidates = i_end - i_start + 1
        if !(get_n_elites(e.population) + n_candidates > e.config.n_elite)
            # Flag all the individuals of this rank
			@inbounds for i in i_start:i_end
				e.population[i].is_elite = true
			end
        else
            # Compute crowding distances and flag less crowded individuals
            assign_crowding_distances!(e, current_rank)
            sort!(e.population, by = ind -> ind.crowding_distance, rev=true)
  			@inbounds for i in 1:(e.config.n_elite - get_n_elites(e.population))
  				e.population[i].is_elite = true
  			end
        end
        current_rank += 1
        i_start = i_end + 1
        i_end = findlast(ind -> ind.rank == current_rank, e.population)
    end
	# Log results
	if ((e.config.log_gen > 0) && (e.gen == 1 || mod(e.gen, e.config.log_gen) == 0))
        log_gen(e, max_fitness, min_fitness)
    end
	# Remove non-elite individuals
	filter!(ind -> ind.is_elite, e.population)
	# Reset flags
	for ind in e.population
        ind.is_elite = false
    end
	#==#
    #=
    # Sort all individuals according to Pareto efficiency
    fast_non_dominated_sort!(e)
    # Select elites as most Pareto efficient solutions
    new_population = Vector{T}()
    sort!(e.population, by = ind -> ind.rank)
    current_rank = 1
    i_start = 1
    i_end = findlast(ind -> ind.rank == current_rank, e.population)
    while length(new_population) < e.config.n_elite
        n_candidates = i_end - i_start + 1
        if !(length(new_population) + n_candidates > e.config.n_elite)
            # Add all the individuals of this rank
            push!(new_population,
                  [NSGA2IndCopy(ind) for ind in e.population[i_start:i_end]]...)
        else
            # Compute crowding distances and add less crowded individuals
            assign_crowding_distances!(e, current_rank)
            sort!(e.population, by = ind -> ind.crowding_distance, rev=true)
            push!(new_population, [NSGA2IndCopy(ind) for ind in
                  e.population[1:e.config.n_elite-length(new_population)]]...)
        end
        current_rank += 1
        i_start = i_end + 1
        i_end = findlast(ind -> ind.rank == current_rank, e.population)
    end
    e.population = new_population
    #new_population
    =#
end

"""
    assign_crowding_distances!(e::NSGA2Evo, rank::Int64)

Assign crowding distances as defined in [1] for individuals with rank
corresponding to the input rank.

[1] Deb, Kalyanmoy, et al. "A Fast and Elitist Multiobjective Genetic
Algorithm: NSGA-II." IEEE transactions on evolutionary computation 6.2
(2002): 182-197.
"""
function assign_crowding_distances!(e::NSGA2Evo, rank::Int64)
    subpop = view(e.population, findall(ind -> ind.rank == rank, e.population))
    for ind in subpop
        ind.crowding_distance = 0
    end
    l = length(subpop)
    @inbounds for i in 1:e.config.d_fitness
        sort!(subpop, by = ind -> ind.fitness[i])
        if subpop[1].fitness[i] != subpop[end].fitness[i]
            subpop[1].crowding_distance = Inf
            subpop[end].crowding_distance = Inf
            cd_norm = subpop[end].fitness[i] - subpop[1].fitness[i]
            for j in 2:l-1
                cd_dist = subpop[j+1].fitness[i] - subpop[j-1].fitness[i]
                subpop[j].crowding_distance += cd_dist / cd_norm
            end
        end
    end
end

"""
    fast_non_dominated_sort!(e::NSGA2Evo)

Fast non-dominated sorting algorithm from the NSGA2 paper [1].

[1] Deb, Kalyanmoy, et al. "A Fast and Elitist Multiobjective Genetic
Algorithm: NSGA-II." IEEE transactions on evolutionary computation 6.2
(2002): 182-197.
"""
function fast_non_dominated_sort!(e::NSGA2Evo)
    # Re-initialize individuals statistics
    for ind in e.population
        ind.rank = 0
        ind.domination_count = 0
        empty!(ind.domination_list)
        ind.crowding_distance = 0.0
    end
    # Set domination lists, counts and 1st ranks
    @inbounds for i in 1:length(e.population)
        for j in i+1:length(e.population)
            if dominates(e.population[i], e.population[j])
                push!(e.population[i].domination_list, j)
                e.population[j].domination_count += 1
            elseif dominates(e.population[j], e.population[i])
                push!(e.population[j].domination_list, i)
                e.population[i].domination_count += 1
            end
        end
        if e.population[i].domination_count == 0
            e.population[i].rank = 1
        end
    end
    # Set ranks higher than 2
    current_rank = 2
    while any([ind.rank == current_rank - 1 for ind in e.population])
        for ind in e.population
            if ind.rank == current_rank - 1
                @inbounds for index in ind.domination_list
                    e.population[index].domination_count -= 1
                    if e.population[index].domination_count == 0
                        e.population[index].rank = current_rank
                    end
                end
            end
        end
        current_rank += 1
    end
end

"""
    init_chomosomes(cfg::NamedTuple)

Initialize a symbolic population of individual with random chromosomes of sizes
corresponding to the input configuration file.
"""
function init_sympop(cfg::NamedTuple)::Vector{SymInd}
    population = Vector{SymInd}(undef, cfg.n_population)
    for i in 1:cfg.n_population
        chromosome = rand_CGPchromosome(cfg)
        population[i] = SymInd(chromosome, i, cfg.type)
    end
    population
end

mutable struct DualCGPGAEvo <: Cambrian.AbstractEvolution
    config::NamedTuple
    logid::String
    resdir::String
    encoder_config::NamedTuple
    encoder_sympop::Vector{SymInd}
    encoder_logger::CambrianLogger
    controller_config::NamedTuple
    controller_sympop::Vector{SymInd}
    controller_logger::CambrianLogger
    fitness::Function
    fitness_matrix::Matrix{Float64}
    elites_matrix::Matrix{Bool}
    eval_matrix::Matrix{Bool}
    eval_mutant::Bool
    n_eval::Int64
    n_elite::Int64
    tournament_size::Int64
    gen::Int64
end

populate(e::IICGP.DualCGPGAEvo) = IICGP.ga_populate(e)
evaluate(e::IICGP.DualCGPGAEvo) = IICGP.fitness_evaluate(e, e.fitness)

"""
    DualCGPGAEvo(
        main_config::Dict,
        encoder_config::NamedTuple,
        controller_config::NamedTuple,
        fitness::Function,
        logid::String,
        resdir::String
    )

Dual CGP evolution framework.
"""
function DualCGPGAEvo(
    main_config::Dict,
    encoder_config::NamedTuple,
    controller_config::NamedTuple,
    fitness::Function,
    logid::String,
    resdir::String
)
    encoder_logfile = joinpath(resdir, "logs", logid, "encoder.csv")
    controller_logfile = joinpath(resdir, "logs", logid, "controller.csv")
    encoder_logger = CambrianLogger(encoder_logfile)
    controller_logger = CambrianLogger(controller_logfile)
    encoder_sympop = init_sympop(encoder_config)
    controller_sympop = init_sympop(controller_config)
    @assert encoder_config.d_fitness == controller_config.d_fitness
    config = (
        d_fitness=encoder_config.d_fitness,
        seed=encoder_config.seed,
        n_gen=max(encoder_config.n_gen, controller_config.n_gen),
        log_gen=min(encoder_config.log_gen, controller_config.log_gen),
        save_gen=min(encoder_config.save_gen, controller_config.save_gen)
    )
    mat_size = (encoder_config.n_population, controller_config.n_population)
    fitness_matrix = -Inf * ones(mat_size...)
    elites_matrix = falses(mat_size...)
    eval_matrix = falses(mat_size...)
    n_eval = main_config["n_eval"]
    n_elite = main_config["n_elite"]
    eval_mutant = main_config["eval_mutant"]
    tournament_size = main_config["tournament_size"]
    DualCGPGAEvo(
        config, logid, resdir,
        encoder_config, encoder_sympop, encoder_logger,
        controller_config, controller_sympop, controller_logger,
        fitness, fitness_matrix, elites_matrix, eval_matrix, eval_mutant,
        n_eval, n_elite, tournament_size, 0
    )
end

mutable struct DualCGPEvolution{T} <: Cambrian.AbstractEvolution
    config::NamedTuple
    logid::String
    resdir::String
    encoder_config::NamedTuple
    encoder_population::Array{T}
    encoder_logger::CambrianLogger
    controller_config::NamedTuple
    controller_population::Array{T}
    controller_logger::CambrianLogger
    fitness::Function
    gen::Int
end

populate(e::IICGP.DualCGPEvolution) = IICGP.oneplus_populate(e)
evaluate(e::IICGP.DualCGPEvolution) = IICGP.fitness_evaluate(e, e.fitness)

"""
    function DualCGPEvolution(
        encoder_config::NamedTuple,
        controller_config::NamedTuple,
        fitness::Function
        resdir::String;;
        bootstrap::Bool=false,
        kwargs...
    )

Dual CGP evolution framework.
"""
function DualCGPEvolution(
    encoder_config::NamedTuple,
    controller_config::NamedTuple,
    fitness::Function,
    resdir::String;
    bootstrap::Bool=false,
    kwargs...
)
    kwargs_dict = Dict(kwargs)

    # Logger
    if haskey(kwargs_dict, :logid)
        logid = kwargs_dict[:logid]
    else
        logid = encoder_config.id
    end
    encoder_logfile = joinpath(resdir, "logs", logid, "encoder.csv")
    controller_logfile = joinpath(resdir, "logs", logid, "controller.csv")
    encoder_logger = CambrianLogger(encoder_logfile)
    controller_logger = CambrianLogger(controller_logfile)

    # CGP populations
    if haskey(kwargs_dict, :encoder_init_function)
        encoder_population = Cambrian.initialize(CGPInd, encoder_config,
            init_function=kwargs_dict[:encoder_init_function])
    else
        encoder_population = Cambrian.initialize(CGPInd, encoder_config)
    end
    if haskey(kwargs_dict, :controller_init_function)
        controller_population = Cambrian.initialize(CGPInd, controller_config,
            init_function=kwargs_dict[:controller_init_function])
    else
        controller_population = Cambrian.initialize(CGPInd, controller_config)
    end

    # Global evolution config
    @assert encoder_config.d_fitness == controller_config.d_fitness
    config = (
        d_fitness=encoder_config.d_fitness,
        seed=encoder_config.seed,
        n_gen=max(encoder_config.n_gen, controller_config.n_gen),
        log_gen=min(encoder_config.log_gen, controller_config.log_gen),
        save_gen=min(encoder_config.save_gen, controller_config.save_gen)
    )

    # Create and return DualCGPEvolution
    DualCGPEvolution(
        config, logid, resdir,
        encoder_config, encoder_population, encoder_logger,
        controller_config, controller_population, controller_logger,
        fitness, 0
    )
end

DualEvo = Union{DualCGPEvolution, DualCGPGAEvo}

function getpop(e::DualEvo; copypop::Bool=false)
    if typeof(e) <: DualCGPGAEvo
        if copypop
            epop = [SymIndCopy(ind) for ind in e.encoder_sympop]
            cpop = [SymIndCopy(ind) for ind in e.controller_sympop]
        else
            epop = e.encoder_sympop
            cpop = e.controller_sympop
        end
    elseif typeof(e) <: DualCGPEvolution
        epop = e.encoder_population
        cpop = e.controller_population
    else
        throw(TypeError("Type $(typeof(e)) not supported in loggen."))
    end
    return epop, cpop
end

function log_gen_from_pop_logger(
    d_fitness::Int64,
    pop::AbstractArray,
    logger::CambrianLogger,
    gen::Int64
)
    for d in 1:d_fitness
        maxs = map(i->i.fitness[d], pop)
        with_logger(logger) do
            @info Formatting.format("{1:04d},{2:e},{3:e},{4:e}",
                                    gen, maximum(maxs), mean(maxs), std(maxs))
        end
    end
    flush(logger.stream)
end

"""
    function log_gen(e::DualCGPEvolution)

Log a generation within a dual CGP evolution framework, including max, mean,
and std of each fitness dimension.
"""
function log_gen(e::DualEvo)
    epop, cpop = getpop(e)
    log_gen_from_pop_logger(e.config.d_fitness, epop, e.encoder_logger, e.gen)
    log_gen_from_pop_logger(e.config.d_fitness, cpop, e.controller_logger, e.gen)
end

function save_gen_at(path::String, pop::AbstractArray)
    sort!(pop)
    for i in eachindex(pop)
        f = open(Formatting.format("{1}/{2:04d}.dna", path, i), "w+")
        write(f, string(pop[i]))
        close(f)
    end
end

"""
    function save_gen(e::DualEvo)

Save the population of a dual CGP evolution framework in `gens/` contained in
the results directory.
"""
function save_gen(e::DualEvo)
    encoder_path = joinpath(e.resdir, Formatting.format("gens/{1}/encoder_{2:04d}", e.logid, e.gen))
    controller_path = joinpath(e.resdir, Formatting.format("gens/{1}/controller_{2:04d}", e.logid, e.gen))
    mkpath(encoder_path)
    mkpath(controller_path)
    epop, cpop = getpop(e; copypop=true)
    save_gen_at(encoder_path, epop)
    save_gen_at(controller_path, cpop)
end

ext: jl


text:
#######################################
#
# FTest
#
#######################################

"""
    ftest(s::SSR)
    ftest(s::SSR; kwargs...)

Run f-test statistics on steady state response measurement.

# Arguments

* `freq_of_interest`: Frequency to analyse for presence of a response
* `side_freq`: Adjacent frequencies above and below the frequency of interest used to quantify the noise
* `ID`: Participant ID for storage in resulting dataframe
* `spill_bins`: Number of bins each side of the FFT bin of interes to ignore in noise computation
* `results_key`: Dictionary key name to store results in `s.processing`

# Examples
```julia
s = read_SSR(fname)
s.modulationrate = 33.2u"Hz"
s = ftest(s)
println(s.processing['statistics'])
```

# Reference
Hofmann, M., Wouters, J. Improved Electrically Evoked Auditory Steady-State Response Thresholds in Humans. JARO 13, 573–589 (2012). https://doi.org/10.1007/s10162-012-0321-8

Luke, Robert, and Jan Wouters. "Kalman filter based estimation of auditory steady state response parameters." IEEE Transactions on Neural Systems and Rehabilitation Engineering 25.3 (2016): 196-204.
"""
function ftest(
    s::SSR;
    freq_of_interest::Union{Real,AbstractArray} = modulationrate(s),
    side_freq::Number = 0.5,
    ID::AbstractString = "",
    spill_bins::Int = 2,
    results_key::AbstractString = "statistics",
    kwargs...,
)

    # Do calculation here once, instead of in each low level call
    spectrum = Neuroimaging._ftest_spectrum(s.processing["sweeps"])
    spectrum = compensate_for_filter(s.processing, spectrum, samplingrate(s))
    frequencies = range(0, stop = 1, length = Int(size(spectrum, 1))) * samplingrate(s) / 2

    for freq in freq_of_interest

        snrDb, phase, signal, noise, statistic =
            ftest(spectrum, frequencies, freq, side_freq, spill_bins)

        result = DataFrame(
            ID = vec(repeat([ID], length(channelnames(s)), 1)),
            Channel = copy(channelnames(s)),
            ModulationRate = copy(modulationrate(s)),
            AnalysisType = vec(repeat(["F-test"], length(channelnames(s)))),
            AnalysisFrequency = vec(repeat([freq], length(channelnames(s)))),
            SignalAmplitude = vec(sqrt.(signal)),
            SignalPhase = vec(phase),
            NoiseAmplitude = vec(sqrt.(noise)),
            SNRdB = vec(snrDb),
            Statistic = vec(statistic),
        )

        result = add_dataframe_static_rows(result, kwargs)

        if haskey(s.processing, results_key)
            s.processing[results_key] = vcat(s.processing[results_key], result)
        else
            s.processing[results_key] = result
        end

    end

    return s
end

#######################################
#
# Helper functions
#
#######################################

"""
Save the processing step information from SSR type.

By default saves information in the key `statistics`,
but this can be modified by the user.
"""
# Save ftest results to file
function save_results(
    a::SSR;
    name_extension::AbstractString = "",
    results_key::AbstractString = "statistics",
    kwargs...,
)

    file_name = string(a.file_name, name_extension, ".csv")

    # Rename to save space
    results = a.processing

    # Index of keys to be exported
    result_idx = find_keys_containing(results, results_key)

    if length(result_idx) > 0

        to_save = get(results, collect(keys(results))[result_idx[1]], 0)

        CSV.write(file_name, to_save)
    end

    @info("File saved to $file_name")

    return a
end


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

name = "HybridLambda"
version = v"0.6.3-beta"

# Collection of sources required to build HybridLambda
sources = [
    "https://github.com/hybridLambda/hybrid-Lambda.git" =>
    "804160711e36143f89667cab405977b2f0812ce6",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd hybrid-Lambda/src/
sed -i.bak 's/g++/${CXX}/' Makefile
make
mkdir ${prefix}/bin
mv hybrid-Lambda ${prefix}/bin/

"""

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

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "hybrid-Lambda", :hybridlambda)
]

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

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: cecileane/HybridLambdaBuilder

max_issues_repo_name: cecileane/HybridLambdaBuilder

max_forks_repo_name: cecileane/HybridLambdaBuilder


text:
"""
    nondimensional_exchange_functions(
        m::AtmosModel{FT},
        entr::EntrainmentDetrainment,
        state::Vars,
        aux::Vars,
        t::Real,
        ts,
        env,
        i,
    ) where {FT}

Returns the nondimensional entrainment and detrainment
functions following Cohen et al. (JAMES, 2020), given:
 - `m`, an `AtmosModel`
 - `entr`, an `EntrainmentDetrainment` model
 - `state`, state variables
 - `aux`, auxiliary variables
 - `t`, the time
 - `ts`, NamedTuple of thermodynamic states
 - `env`, NamedTuple of environment variables
 - `i`, the updraft index
"""
function nondimensional_exchange_functions(
    m::AtmosModel{FT},
    entr::EntrainmentDetrainment,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    env,
    i,
) where {FT}

    # Alias convention:
    gm = state
    up = state.turbconv.updraft
    up_aux = aux.turbconv.updraft
    en_aux = aux.turbconv.environment

    # precompute vars
    w_min = entr.w_min
    N_up = n_updrafts(m.turbconv)
    ρinv = 1 / gm.ρ
    a_up_i = up[i].ρa * ρinv
    w_up_i = up[i].ρaw / up[i].ρa

    # thermodynamic variables
    RH_up = relative_humidity(ts.up[i])
    RH_en = relative_humidity(ts.en)

    Δw = filter_w(w_up_i - env.w, w_min)
    Δb = up_aux[i].buoyancy - en_aux.buoyancy

    c_δ = sign(condensate(ts.en) + condensate(ts.up[i])) * entr.c_δ

    # compute dry and moist aux functions
    μ_ij = (entr.χ - a_up_i / (a_up_i + env.a)) * Δb / Δw
    D_ε = entr.c_ε / (1 + exp(-μ_ij / entr.μ_0))
    M_ε = c_δ * (max((RH_en^entr.β - RH_up^entr.β), 0))^(1 / entr.β)
    D_δ = entr.c_ε / (1 + exp(μ_ij / entr.μ_0))
    M_δ = c_δ * (max((RH_up^entr.β - RH_en^entr.β), 0))^(1 / entr.β)
    return D_ε, D_δ, M_δ, M_ε
end;

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Test

# @testset Collision: Sphere-sphere begin

# Parameters
origin = Origin{Float64}()
pbody = Sphere(0.5, 1.0)
cbody = Sphere(0.5, 1.0)
joint = JointConstraint(Fixed(origin, pbody))

bodies = [pbody, cbody]
joints = [joint]

collision = SphereSphereCollision{Float64,2,3,6}(
        szeros(3),
        szeros(3),
        pbody.shape.r,
        cbody.shape.r)
friction_parameterization = SA{Float64}[
    1.0  0.0
    0.0  1.0
]
body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

contacts = [ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)]

mech = Mechanism(origin, bodies, joints, contacts,
            gravity=1.0 * -9.81,
            timestep=0.1)

mech.bodies[1].state.x2 = [0.0, 0.0, 0.0]
mech.bodies[2].state.x2 = [0.0, 0.0, 2.0]

d = distance(mech.contacts[1].model.collision,
    mech.bodies[1].state.x2, mech.bodies[1].state.q2,
    mech.bodies[2].state.x2, mech.bodies[2].state.q2,)

@test abs(d - 1.0) < 1.0e-6

cp = contact_point(:parent, mech.contacts[1].model.collision,
    mech.bodies[1].state.x2, mech.bodies[1].state.q2,
    mech.bodies[2].state.x2, mech.bodies[2].state.q2,)

@test norm(cp - [0.0; 0.0; 0.5], Inf) < 1.0e-6

cc = contact_point(:child, mech.contacts[1].model.collision,
    mech.bodies[1].state.x2, mech.bodies[1].state.q2,
    mech.bodies[2].state.x2, mech.bodies[2].state.q2,)

@test norm(cc - [0.0; 0.0; 1.5], Inf) < 1.0e-6

cn = contact_normal(mech.contacts[1].model.collision,
    mech.bodies[1].state.x2, mech.bodies[1].state.q2,
    mech.bodies[2].state.x2, mech.bodies[2].state.q2,)

@test norm(cn - [0.0 0.0 -1.0], Inf) < 1.0e-6

ct = contact_tangent(mech.contacts[1].model.collision,
    mech.bodies[1].state.x2, mech.bodies[1].state.q2,
    mech.bodies[2].state.x2, mech.bodies[2].state.q2,)

@test norm(ct - [0.0 1.0 0.0; -1.0 0.0 0.0], Inf) < 1.0e-6

am = adjacency_matrix(mech.joints, mech.bodies, mech.contacts)

@test norm(am - [0 1 0 0; 1 0 0 1; 0 0 0 1; 0 1 1 0], Inf) < 1.0e-6

storage = simulate!(mech, 2.0,
    verbose=false,
    record=true)

vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis)


text:
function checkInput()

  if haskey(ENV, "PDESOLVER_BUNDLE_DEPS")

    if !haskey(ENV, "PDESOLVER_PKGDIR")
      error("PDESOLVER_PKGDIR environment variable must be specified when bundling dependencies")
    end

    pkgdir = ENV["PDESOLVER_PKGDIR"]
    if !isdir(pkgdir)
      error("PDESOLVER_PKGDIR $pkgdir must exist")
    end

    if length(readdir(pkgdir)) != 0
      error("PDESOLVER_PKGDIR $pkgdir must be empty directory")
    end

  end

  if haskey(ENV, "PDESOLVER_UNBUNDLE_DEPS")

    if haskey(ENV, "PDESOLVER_PKGDIR")
      error("Environment variable PDESOLVER_PKGDIR should not be set when unbundling dependencies (PDESOLVER_UNBUNDLE_DEPS).  Copy the contents of the directory directory PDESOLVER_PKGDIR to the Julia package directory (eg. ~/.julia/v0.6) on the new machine")
    end
  end

  return nothing
end

"""
  This function determines if a git identifier is a commit or something else
  (ie. a branch or tag)

  **Inputs**

   * pkg_dir: path to directory containing the repo
   * hash: git object identifier

  **Outputs**

   * bool, true if hash represents a commit, false otherwise
"""
function is_commit(pkg_dir::AbstractString, hash::AbstractString)

  start_dir = pwd()
  if pkg_dir != ""
    cd(pkg_dir)
  end

  iscommit = false
  try
    vals = readstring(`git show-ref \| grep $hash`)
    # make sure this is an exact match
    name = split(vals, '/')[end]
    if name[1:end-1] == hash
      iscommit = false
    else
      iscommit = true
    end
  catch x
    iscommit = true
  end

  cd(start_dir)
  return iscommit
end

"""
  Checks out a git object of some kind.  If the object is a tag or a branch,
  it is checkout out as normal.  If it is a commit, a new branch is created
  called detatched_from_commit, where commit is the commit hash

  **Inputs**

   * pkg_dir: directory where the git repo is located
   * hash: git object identifier
"""
function set_hash(pkg_dir::AbstractString, hash::AbstractString)
# this function checks out a particular git hash of a particular package

  new_branchname = "detached_from_$hash"

  start_dir = pwd()
  cd(pkg_dir)
#  branch_name = readstring(`git rev-parse --abbrev-ref HEAD`)
#  if branch_name[1:end-1] != new_branchname
#    run(`git checkout -b pdesolver_version $hash`)
#  end
  iscommit = is_commit("", hash)
  run(`git checkout $hash`)

  if iscommit
    run(`git checkout -b $new_branchname HEAD`)
  end
  cd(start_dir)

  println("finished setting hash")
end

"""
  This function installs a specified package.  Unless overridded, if the package is
  already installed, it will not be modified.

  **Inputs**

  * dir: directory in which to install package
  * pkg_name: name of the package (not including .jl)
  * git_url: URL of git repository
  * git_commit: anything that git checkout can use to identify what to check out
  * pkg_dict: dictionary of installed package
  * f: an IO object to write info to
  * force: force package to install even if already present, (keyword arg), default
           false
"""
function install_pkg(dir::AbstractString, pkg_name::AbstractString, git_url::AbstractString, git_commit::AbstractString, pkg_dict, f; force=false)
# pkg_name = name of package
# git_url: value to be passed to Pkg.clone()
# git_commit: commit has to check out
# pkg_dict: output of Pkg.installed, passed as an argument because it takes so long to get
# f: file handle to write progress to
# force: install the package regardless of other flags, default false

  already_installed = haskey(pkg_dict, pkg_name)
  pkg_path = joinpath(dir, pkg_name)
  force_specific = haskey(ENV, "PDESOLVER_FORCE_DEP_INSTALL_$pkg_name")
  if !already_installed || force_specific || FORCE_INSTALL_ALL || force
    println(f, "Installing package $pkg_name")
    try
      if !isdir(pkg_path)  # if package not already present
        start_dir = pwd()
        cd(dir)
        run(`git clone $git_url`)
        name_ext = string(pkg_name, ".jl")
        run(`mv -v ./$name_ext ./$pkg_name`)
        set_hash(pkg_path, git_commit)
        cd(start_dir)
      end
#      else
#        reset_repo(pkg_name)
#      end
      println(f, "Resolving dependencies")
      Pkg.resolve()
      println(f, "Building package")
      println(f, "  Note: if this step fails, there is no way to report it in this log")
      Pkg.build(pkg_name)
      Libc.flush_cstdio()
      println(f, "  Installation appears to have completed sucessfully")
    catch x
      Libc.flush_cstdio()
      println(f, "Error installing package $pkg_name")
      println(f, "Error is $x")
    end

  else
    println(f, "Skipping installation of package $pkg_name")
    println(f, "already installed: ", already_installed, ", specifically forced: ", force_specific,
               ", force general: ", FORCE_INSTALL_ALL, ", force override: ", force)
  end

  flush(f)
end

function reset_repo(pkg_name::AbstractString)
  start_dir = pwd()
  cd(Pkg.dir(pkg_name))
  run(`git reset --hard`)  # remove changes to files under version control
  run(`git clean -x -f -d`) # remove any files not under version control
end

"""
  This function downloads packages to a specified directory, in preparation for
  copying them to a machine without internet access.

  **Inputs**

   * dir: the directory to put the packages into
   * pkg_name: the name of the package
   * git_url: the git URL of the package
   * git_commit: the name of a commit/tag/branch to checkout
   * f: an IO to write logging messages to
"""
function bundle_pkg(dir::AbstractString, pkg_name::AbstractString, git_url::AbstractString, git_commit::AbstractString, f::IO)
# clone all packages into a specified directory, without installing any of them
  println(f, "Bundling package $pkg_name")

  start_dir = pwd()
  cd(dir)

  println(f, "now in directory ", pwd())
  # check if package already exists
  if isdir(pkg_name)
    println(f, "package $pkg_name is already present, skipping...")
    return
  end

  try
    run(`git clone $git_url`)
    println(f, "ls of pwd() = ", readstring(`ls`))
    name_ext = string(pkg_name, ".jl")
    run(`mv -v ./$name_ext ./$pkg_name`)
    set_hash("./$pkg_name", git_commit)
    println(f, "set hash to ", git_commit)

    if isdir("./$pkg_name/deps")
      cd("./$pkg_name/deps")

      # perform any needed downloads
      if isfile("./download.sh")
        run(`./download.sh`)
      end
    end

    println(f, "Bundling of package $pkg_name appears to have completed successfully")
  catch x
    println(f, "Error bundling package $pkg_name")
    println(f, "error is $x")
  end

  flush(f)
  cd(start_dir)
end

"""
  Unbundle packages and install them, without accessing the internet.
  This is the inverse function to `bundle_pkg()`.

  For this to work, the contents of the directory created by
  PDESOLVER_BUNDLE_DEPS should be copied to the new machine and placed in the
  Julia package directory (ie. ~/.julia/v0.6 or JULIA_PKGDIR).  When done,
  the contents of this directory should look like

  ~/.julia/v0.6:
    - PDESolver
    - ODLCommonTools
    - PumiInterface
    - etc...

  Note that the order in which the dependencies is installed is important.
  They must be installed in an order such that all the dependencies of a
  package are installed before the package itself

  The order in which packages are installed is determined by the order in
  which they appear listed in `build.jl`.

  **Inputs**

   * dir: the julia package directory
   * pkg_name: the name of the package
   * f: an IO object to write logging messages to
"""
function unbundle_pkg(dir::AbstractString, pkg_name::AbstractString, f::IO)
# install a package by directly calling the build.jl script
# this bypasses the REQUIRE file, and therefore any internet access
# the dependencies must be installed in order for this to work

  println(f, "unbundling package $pkg_name")
  start_dir = pwd()
  pkgdir = joinpath(dir, pkg_name)
  cd(pkgdir)
  if isdir("./deps")
    cd("./deps")

    if isfile("./build.jl")
      try
        buildpath = joinpath(dir, pkg_name, "deps", "build.jl")
        include(buildpath)
        println(f, "unbundling package $pkg_name appears to have completed successfully")
      catch x
        println(f, "Error unbundling package $pkg_name")
        println(f, "Error is $x")
      end

    else
      println(f, "No build.jl file for this package, skipping...")

    end   # end if isfile()

  else
    println(f, "No deps directory for this package, skipping...")
  end  # end isdir()

  flush(f)
end


text:
using TransformVariables, Parameters, Statistics, StatsFuns, Optim
using NLSolversBase

function makeLoss(model)

    t = getTransform(model)

    fpre = @eval $(logdensity(model))
    f(par, data) = Base.invokelatest(fpre, par, data)

    loss(x, data) = -f(t(x), data)

    (loss=loss, t=t)
end

export getMAP
function getMAP(m :: Model ;kwargs...)
    @unpack loss, t = makeLoss(m)
    d = dimension(t)

    kwargs = Dict(kwargs)
    # @show kwargs
    data = get(kwargs, :data, NamedTuple{}())
    init = get(kwargs, :init, t(zeros(d)))

    f(x) = loss(x, data)

    opt = optimize(f, inverse(t)(init),method=LBFGS())
    t(Optim.minimizer(opt))
end


text:
using Test
using Dates
using UnetSockets

# start fjåge

println("Starting 2-node simulation...")
run(`bash sim.sh start`)
sleep(5)

# tests

println("Starting tests...")
try

  @testset "socket" begin
    s = UnetSocket("localhost", 1101)
    @test getlocaladdress(s) == 232
    @test host(s, "A") == 232
    @test host(s, "B") == 31
    @test !isclosed(s)
    close(s)
    @test isclosed(s)
  end

  @testset "gateway" begin
    s = UnetSocket("localhost", 1101)
    gw = getgateway(s)
    @test gw isa Fjage.Gateway
    shell = agentforservice(gw, Services.SHELL)
    @test shell isa AgentID
    @test shell.language == "Groovy"
    close(s)
  end

  @testset "agents" begin
    s = UnetSocket("localhost", 1101)
    node = agent(s, "node")
    @test node isa AgentID
    @test node.address == 232
    @test node.nodeName == "A"
    phy = agentforservice(s, Services.PHYSICAL)
    @test phy isa AgentID
    @test phy.name == "phy"
    @test phy.MTU > 0
    close(s)
  end

  @testset "bind/connect" begin
    s = UnetSocket("localhost", 1101)
    @test getlocalprotocol(s) == -1
    @test !isbound(s)
    @test bind(s, 42)
    @test isbound(s)
    @test getlocalprotocol(s) == 42
    unbind(s)
    @test !isbound(s)
    @test getlocalprotocol(s) == -1
    @test getremoteaddress(s) == -1
    @test getremoteprotocol(s) == 0
    @test !isconnected(s)
    @test connect(s, 31, 0)
    @test getremoteaddress(s) == 31
    @test getremoteprotocol(s) == 0
    @test isconnected(s)
    disconnect(s)
    @test !isconnected(s)
    @test getremoteaddress(s) == -1
    @test getremoteprotocol(s) == 0
    close(s)
  end

  @testset "timeouts" begin
    s = UnetSocket("localhost", 1101)
    @test bind(s, 0)
    @test gettimeout(s) == -1
    settimeout(s, 1000)
    @test gettimeout(s) == 1000
    t1 = now().instant
    @test receive(s) == nothing
    dt = now().instant - t1
    @test dt.value > 1000
    settimeout(s, 0)
    @test gettimeout(s) == 0
    t1 = now().instant
    @test receive(s) == nothing
    dt = now().instant - t1
    @test dt.value < 500
    close(s)
  end

  @testset "cancel" begin
    s = UnetSocket("localhost", 1101)
    @test bind(s, 0)
    @test gettimeout(s) == -1
    t1 = now().instant
    @async begin
      sleep(2)
      cancel(s)
    end
    @test receive(s) == nothing
    dt = now().instant - t1
    @test dt.value > 1500
    close(s)
  end

  @testset "communication" begin
    s1 = UnetSocket("localhost", 1101)
    s2 = UnetSocket("localhost", 1102)
    @test bind(s2, Protocol.USER)
    settimeout(s2, 1000)
    @test !send(s1, [1,2,3])
    @test send(s1, [1,2,3], 31)
    @test receive(s2) == nothing
    @test send(s1, [1,2,3], 31, Protocol.USER)
    ntf = receive(s2)
    @test ntf isa DatagramNtf
    @test ntf.data == [1,2,3]
    @test connect(s1, 31, Protocol.USER)
    @test send(s1, [1,2,3])
    ntf = receive(s2)
    @test ntf isa DatagramNtf
    @test ntf.data == [1,2,3]
    @test send(s1, [1,2,3], 31, 0)
    @test receive(s2) == nothing
    @test send(s1, [1,2,3], 27, Protocol.USER)
    @test receive(s2) == nothing
    @test send(s1, [1,2,3])
    ntf = receive(s2)
    @test ntf isa DatagramNtf
    @test ntf.data == [1,2,3]
    disconnect(s1)
    @test !send(s1, [1,2,3])
    @test send(s1, [1,2,3], 31, Protocol.USER)
    ntf = receive(s2)
    @test ntf isa DatagramNtf
    @test ntf.data == [1,2,3]
    close(s1)
    close(s2)
  end

finally

# stop simulator

  println("Stopping 2-node simulation...")
  p = run(`bash sim.sh stop`; wait=false)
  sleep(1)
  success(p) || println("Could not stop!")

end


text:
using ITensors
using ITensors.NDTensors
using Random
using Test
import ITensors: In, Out, Neither

@testset "Index" begin
  @testset "Index with dim" begin
    i = Index(2)
    @test id(i) != 0
    @test dim(i) == 2
    @test dir(i) == Neither
    @test plev(i) == 0
    @test tags(i) == TagSet("")
  end
  @testset "Index with all args" begin
    i = Index(1, 2, In, "Link", 1)
    @test id(i) == 1
    @test dim(i) == 2
    @test dir(i) == In
    @test plev(i) == 1
    @test tags(i) == TagSet("Link")
    j = copy(i)
    @test id(j) == 1
    @test dim(j) == 2
    @test dir(j) == In
    @test plev(j) == 1
    @test tags(j) == TagSet("Link")
    @test j == i
  end
  @testset "prime" begin
    i = Index(2)
    @test plev(i) == 0
    i2 = prime(i, 2)
    @test plev(i2) == 2
    i1 = i'
    @test plev(i1) == 1
    i2 = i''
    @test plev(i2) == 2
    i3 = i'''
    @test plev(i3) == 3
    i6 = i^6
    @test plev(i6) == 6
    i0 = noprime(i)
    @test plev(i0) == 0
  end
  @testset "IndexVal" begin
    i = Index(2)
    @test_deprecated i[1]
    @test val(i => 1) == 1
    @test ind(i => 1) == i
    @test isindequal(i, i => 2)
    @test isindequal(i => 2, i)
    @test plev(i' => 2) == 1
    @test val(i' => 2) == 2
    @test plev(prime(i => 2, 4)) == 4

    @test plev(i => 2) == 0
    @test plev(i' => 2) == 1
    @test prime(i => 2) == (i' => 2)
    iv = i => 2
    ĩv = sim(i => 2)
    @test ind(iv) ≠ ind(ĩv)
    @test val(iv) == val(ĩv)
  end
  @testset "Iteration" begin
    i = Index(3)

    c = 1
    for iv in eachindval(i)
      @test iv == (i => c)
      c += 1
    end

    c = 1
    for n in eachval(i)
      @test n == c
      c += 1
    end
  end
  @testset "Broadcasting" begin
    N = 3
    i = Index(2)
    ps = (n - 1 for n in 1:4)
    is = prime.(i, ps)
    @test is[1] == i
    @test is[2] == i'
    @test is[3] == i''
    ts = ("i$n" for n in 1:4)
    is = settags.(i, ts)
    @test is[1] == addtags(i, "i1")
    @test is[2] == addtags(i, "i2")
    @test is[3] == addtags(i, "i3")
  end
  @testset "Index ID random seed" begin
    Random.seed!(index_id_rng(), 1234)
    i = Index(2)
    j = Index(2)
    Random.seed!(index_id_rng(), 1234)
    ĩ = Index(2)
    j̃ = Index(2)
    Random.seed!(index_id_rng(), 123)
    ĩ′ = Index(2)
    j̃′ = Index(2)
    @test id(i) == id(ĩ)
    @test id(j) == id(j̃)
    @test id(i) ≠ id(ĩ′)
    @test id(j) ≠ id(j̃′)

    Random.seed!(index_id_rng(), 1234)
    Random.seed!(1234)
    i = Index(2)
    j = Index(2)
    Random.seed!(1234)
    ĩ = Index(2)
    j̃ = Index(2)
    Random.seed!(123)
    ĩ′ = Index(2)
    j̃′ = Index(2)
    @test id(i) ≠ id(ĩ)
    @test id(j) ≠ id(j̃)
    @test id(i) ≠ id(ĩ′)
    @test id(j) ≠ id(j̃′)

    Random.seed!(index_id_rng(), 1234)
    Random.seed!(1234)
    i = Index(2)
    j = Index(2)
    A = randomITensor(i, j)
    Random.seed!(1234)
    ĩ = Index(2)
    j̃ = Index(2)
    Ã = randomITensor(ĩ, j̃)
    Random.seed!(123)
    ĩ′ = Index(2)
    j̃′ = Index(2)
    Ã′ = randomITensor(ĩ′, j̃′)
    @test id(i) ≠ id(ĩ)
    @test id(j) ≠ id(j̃)
    @test id(i) ≠ id(ĩ′)
    @test id(j) ≠ id(j̃′)
    @test all(tensor(A) .== tensor(Ã))
    @test all(tensor(A) .≠ tensor(Ã′))

    Random.seed!(index_id_rng(), 1234)
    Random.seed!(1234)
    i = Index(2)
    j = Index(2)
    A = randomITensor(i, j)
    Random.seed!(index_id_rng(), 1234)
    Random.seed!(1234)
    ĩ = Index(2)
    j̃ = Index(2)
    Ã = randomITensor(ĩ, j̃)
    Random.seed!(index_id_rng(), 1234)
    Random.seed!(123)
    ĩ′ = Index(2)
    j̃′ = Index(2)
    Ã′ = randomITensor(ĩ′, j̃′)
    @test id(i) == id(ĩ)
    @test id(j) == id(j̃)
    @test id(i) == id(ĩ′)
    @test id(j) == id(j̃′)
    @test all(tensor(A) .== tensor(Ã))
    @test all(tensor(A) .≠ tensor(Ã′))
  end
  @testset "directsum" begin
    i = Index(2, "i")
    j = Index(3, "j")
    ij = ITensors.directsum(i, j; tags="test")
    @test dim(ij) == 5
    @test hastags(ij, "test")
  end
end

nothing


text:
export FiniteException

###
### Exceptions
###
struct KeywordArgumentException <: Exception
    key::Symbol
    given::Any
    msg::String
end
function KeywordArgumentException(key, given)
    KeywordArgumentException(key, given, "")
end
function Base.showerror(io::IO, E::KeywordArgumentException)
    print(io, "Invalid argument $given for $key. ", msg)
end

struct FiniteException <: Exception
    dim::Int
end
function Base.showerror(io::IO, E::FiniteException)
    print(
        io,
        "FiniteException: The solution set of the given system has at least dimension " *
        "$(E.dim) > 0. Consider intersecting with an (affine) subspace of codimension " *
        "$(E.dim) to reduce to (possibly) finitely many solutions.",
    )
end

@noinline function unsupported_kwargs(kwargs)
    if !(isempty(kwargs))
        msg = join(["$k = $v" for (k, v) in pairs(kwargs)], ", ")
        @warn "Ingored unsupported keyword arguments: $msg"
    end
end

"""
    rand_unitary_matrix(n::Int, T=ComplexF64)

Samples a `n × n` unitary Matrix uniformly from the space of all unitary n × n matrices.

See https://arxiv.org/abs/math-ph/0609050 for a derivation.
"""
function rand_unitary_matrix(n::Int, T::Type = ComplexF64)
    Z = randn(T, n, n) ./ sqrt(2)
    Q, R = LA.qr(Z)
    Λ = LA.diagm(0 => [R[i, i] / abs(R[i, i]) for i = 1:n])
    Q * Λ
end

fast_abs(z::Complex) = sqrt(abs2(z))
fast_abs(x::Real) = abs(x)

"Like `min(a,b)`` but ignoring any `NaN` values."
nanmin(a, b) = isnan(a) ? b : (isnan(b) ? a : min(a, b))

"""
    unpack(a::Union{Nothing, T}, b::T)

Returns `a` if it is not `nothing`, otherwise `b`.
"""
unpack(a::Union{Nothing,T}, b::T) where {T} = a === nothing ? b : a

plural(singularstr, n) = n == 1 ? singularstr : singularstr * "s"

"""
     print_fieldnames(io::IO, obj)

 A better default printing for structs.
 """
function print_fieldnames(io::IO, obj)
    println(io, typeof(obj), ":")
    for name in fieldnames(typeof(obj))
        if getfield(obj, name) !== nothing
            val = getfield(obj, name)
            print(io, " • ", name, " → ")
            if val isa AbstractFloat
                println(io, round(val; sigdigits = 5))
            else
                println(io, val)
            end
        end
    end
end

mutable struct SegmentStepper
    start::ComplexF64
    target::ComplexF64
    abs_Δ::Float64
    forward::Bool
    # current
    s::Float64
    # proposed
    s′::Float64
end
SegmentStepper(start::Number, target::Number) =
    SegmentStepper(ComplexF64(start), ComplexF64(target))
SegmentStepper(start::ComplexF64, target::ComplexF64) =
    init!(SegmentStepper(start, target, 0.0, true, 0.0, 0.0), start, target)

init!(S::SegmentStepper, start::Number, target::Number) =
    init!(S, ComplexF64(start), ComplexF64(target))
function init!(S::SegmentStepper, start::ComplexF64, target::ComplexF64)
    S.start = start
    S.target = target
    S.abs_Δ = abs(target - start)
    S.forward = abs(start) < abs(target)
    S.s = S.s′ = S.forward ? 0.0 : S.abs_Δ
    S
end

is_done(S::SegmentStepper) = S.forward ? S.s == S.abs_Δ : S.s == 0.0

step_success!(S::SegmentStepper) = (S.s = S.s′; S)
function propose_step!(S::SegmentStepper, Δs::Real)
    if S.forward
        S.s′ = min(S.s + Δs, S.abs_Δ)
    else
        S.s′ = max(S.s - Δs, 0.0)
    end
    S
end
dist_to_target(S::SegmentStepper) = S.forward ? S.abs_Δ - S.s : S.s

function Base.getproperty(S::SegmentStepper, sym::Symbol)
    if sym == :Δs
        s = getfield(S, :s)
        s′ = getfield(S, :s′)
        Δs = getfield(S, :forward) ? s′ - s : s - s′
        return Δs
    elseif sym === :t
        s = getfield(S, :s)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        return _t_helper(start, target, s, abs_Δ, forward)
    elseif sym == :t′
        s′ = getfield(S, :s′)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        return _t_helper(start, target, s′, abs_Δ, forward)
    elseif sym == :Δt
        s = getfield(S, :s)
        s′ = getfield(S, :s′)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        if forward
            Δt = ((s′ - s) / abs_Δ) * (target - start)
        else
            Δt = ((s - s′) / abs_Δ) * (target - start)
        end
        return Δt
    else # fallback to getfield
        return getfield(S, sym)
    end
end
function _t_helper(start, target, s, Δ, forward)
    if forward
        if s == 0.0
            return start
        elseif s == Δ
            return target
        else
            return start + (s / Δ) * (target - start)
        end
    else
        if s == Δ
            return start
        elseif s == 0.0
            return target
        else
            return target + (s / Δ) * (start - target)
        end
    end
end

function Base.show(io::IO, val::SegmentStepper)
    print(io, "SegmentStepper:")
    for field in [:start, :target, :t, :Δt]
        print(io, "\n • ", field, " → ", getproperty(val, field))
    end
end

"""
    nthroot(x, n)

Compute the `n`-th root of `x`.
"""
function nthroot(x::Real, N::Integer)
    if N == 4
        √(√(x))
    elseif N == 2
        √(x)
    elseif N == 3
        cbrt(x)
    elseif N == 1
        x
    elseif N == 0
        one(x)
    else
        x^(1 / N)
    end
end

import Base: MPFR
import Base.MPFR: MPFRRoundingMode

function set!(x::BigFloat, y::BigFloat, r::MPFRRoundingMode = MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_set, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
        x,
        y,
        r,
    )
    x
end

function set!(x::BigFloat, d::Float64, r::MPFRRoundingMode = MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_set_d, :libmpfr),
        Int32,
        (Ref{BigFloat}, Float64, MPFRRoundingMode),
        x,
        d,
        r,
    )
    x
end

# Basic arithmetic without promotion
for (fJ, fC) in ((:add!, :add), (:mul!, :mul))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(
                ($(string(:mpfr_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
                z,
                x,
                y,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        function ($fJ)(z::BigFloat, x::BigFloat, c::Float64)
            ccall(
                ($(string(:mpfr_, fC, :_d)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Cdouble, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::Float64, x::BigFloat) = ($fJ)(x, c)

        function ($fJ)(z::BigFloat, x::BigFloat, c::Int64)
            ccall(
                ($(string(:mpfr_, fC, :_si)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Cdouble, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::Int64, x::BigFloat) = ($fJ)(x, c)

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(
                ($(string(:mpfr_, fC, :_z)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigInt}, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::BigInt, x::BigFloat) = ($fJ)(x, c)
    end
end

for (fJ, fC) in ((:sub!, :sub), (:div!, :div))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(
                ($(string(:mpfr_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
                z,
                x,
                y,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        # Float32/Float64
        function ($fJ)(z::BigFloat, x::BigFloat, c::Float64)
            ccall(
                ($(string(:mpfr_, fC, :_d)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Cdouble, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        function ($fJ)(z::BigFloat, c::Float64, x::BigFloat)
            ccall(
                ($(string(:mpfr_, :d_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Cdouble, Ref{BigFloat}, MPFRRoundingMode),
                z,
                c,
                x,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(
                ($(string(:mpfr_, fC, :_z)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigInt}, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        # no :mpfr_z_div function
    end
end

function sub!(z::BigFloat, c::BigInt, x::BigFloat)
    ccall(
        (:mpfr_z_sub, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigInt}, Ref{BigFloat}, MPFRRoundingMode),
        z,
        c,
        x,
        MPFR.ROUNDING_MODE[],
    )
    return z
end

function rem!(z::BigFloat, x::BigFloat, y::BigFloat, ::RoundingMode{:Nearest})
    ccall(
        (:mpfr_remainder, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
        z,
        x,
        y,
        MPFR.ROUNDING_MODE[],
    )
    return z
end


text:
import Base.copy

mutable struct tabgraph
	nv::Int
	ne::Int
	adj::Array{Bool,2}
	pds::Array{Real,2}
	pred::Array{Int,2}
	tabgraph(nbs::Int) = new(nbs,0,zeros(Bool,(nbs,nbs)),Array{Real}(undef,nbs,nbs),Array{Int}(undef,nbs,nbs))
end

function copy(G::tabgraph)
	F=tabgraph(G.nv)
	F.adj=copy(G.adj)
	F.ne=G.ne
	F.pds=copy(G.pds)
	F.pred=copy(G.pred)
	return F
end

function relie!(G,s,d,p)

	if G.adj[s,d] != 1
		G.ne += 1
	end
	G.adj[s,d] = 1
	G.pds[s,d] = p
end

function aff(G)
	println("Nombre de sommets : ", G.nv)
	println("Nombre d'arrêtes  : ", G.ne)

	for  i in 1:G.nv
		for j in 1:G.nv
			if G.adj[i,j]nb_tirage = rand(1:G.nv)
	preds = zeros(nb_tirage)
	suivs = zeros(nb_tirage)

				println("(",i,",",j,",",G.pds[i,j],")")
			end
		end
	end
end

function alea!(G)

	for i in 1:G.nv
		for j in 1:G.nv
			if rand(1:4) == 1
				relie!(G,i,j,rand())
			end
		end
	end
end

function voisins(G,s)
	ret = Array{Integer}()
	for j in 1:
		if G.adj[s,j]
			push!(ret, j)
		end
	end
	return ret
end

function lire!(G,nomf)

	f = open(nomf)
	while l = readline() != EOF
		arette = parse(l," ")
		relie!(G,arette[1],arette[2],arette[3])
	end

end

function ecrire!(G)

end


text:
using LinearAlgebraicRepresentation, SparseArrays
Lar = LinearAlgebraicRepresentation;
using IntervalTrees,LinearAlgebra
# using Revise, OhMyREPL

function settestpoints2d(W,EV,FV,f, copEV,copFE) # W by rows
	e = findnz(copFE[f,:])[1][1] # first edge of face f
	v1,v2 = findnz(copEV[e,:])[1] # two (global) verts incident on in
    t = W[v2,:] - W[v1,:]
    n = [-t[2],t[1]]
	p0 = (W[v1,:] + W[v2,:]) ./ 2
	ϵ = 1.0e-4
	ptest1 = p0 + ϵ * n
	ptest2 = p0 - ϵ * n
	return ptest1, ptest2
end

function getinternalpoint2d(W,EV,FV, f, copEV,copFE) # W by rows
	#edges for v1=FV[1][1]
	ptest1, ptest2 = Lar.settestpoints2d(W,EV,FV, f, copEV,copFE)
    edges = [findnz(copEV[e,:])[1] for e in findnz(copFE[f,:])[1]]
    V = convert(Lar.Points,W')
    classify = Lar.pointInPolygonClassification(V,edges)
    if classify(ptest1) == "p_in"
        return ptest1
    elseif classify(ptest2) == "p_in"
        return ptest2
    else
        error("classifying inner point in face $f")
    end
end

function chainbasis2polygons(V,copEV,copFE)
	FE = [findnz(copFE[k,:])[1] for k=1:copFE.m]
	EV = [findnz(copEV[k,:])[1] for k=1:copEV.m]

	FEs = Array{Int64,1}[]
	EVs = Array{Array{Int64,1},1}[]
	FVs = Array{Int64,1}[]
	for f=1:copFE.m
		push!( FEs, collect(Set(cat([e for e in FE[f]]))) )
		# edges in EVs are aggregated by face, in order to answer point-classifications
		push!( EVs, [EV[e] for e in FE[f]]  )
		push!( FVs, collect(Set(cat([EV[e] for e in FE[f]])))  )
	end
	polygons = collect(zip(EVs,FVs,FEs)) # (EV,FV,FFE)s for polygon
	W = convert(Lar.Points,V')
	return W,polygons,FE
end

function internalpoints2d(W,copEV,copFE) # W by rows
	# transform each 2-cell in a solid (via Lar model)
	#-------------------------------------------------------------------------------
	U,pols,FE = Lar.chainbasis2polygons(W,copEV,copFE) # V by rows
	# compute, for each `pol` (3-cell) in `pols`, one `internalpoint`.
	#-------------------------------------------------------------------------------
	internalpoints = []
	for f=1:length(pols)
		(EV,FV,FE) = pols[f]
		#GL.VIEW([ GL.GLFrame, GL.GLLines(V,EV) ]);
		internalpoint = Lar.getinternalpoint2d(W,EV,FV, f, copEV,copFE)
		push!(internalpoints,internalpoint)
	end
	return internalpoints
end

"""
	testinternalpoint2d(V::Lar.Points, EV::Lar.Cells, FV::Lar.Cells)

"""
function testinternalpoint2d(listOfModels)
	function testinternalpoint0(testpoint)
		intersectedfaces = Int64[]
		depot = []
		# actual containment test of ray point in faces within depot
		for (k,model) in enumerate(listOfModels)
            verts,edges = model
			classify = Lar.pointInPolygonClassification(verts,edges)
			inOut = classify(testpoint)
			if inOut == "p_in"
				push!(intersectedfaces,k)
			end
		end
		return intersectedfaces
	end
	return testinternalpoint0
end

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

function bool2d(assembly)
	# input of affine assembly
	#-------------------------------------------------------------------------------
	# V,EV = Lar.struct2lar(assembly) #TODO proper different method
	V,EV = Lar.struct2lar(assembly)
	cop_EW = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
	W = convert(Lar.Points, V');
	# generate the 3D space arrangement
	#-------------------------------------------------------------------------------
	W, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
	#V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)
	innerpoints = Lar.internalpoints2d(W,copEV,copFE[1:end,:])
	# associate internal points to 2-cells
	#-------------------------------------------------------------------------------
	listOfModels = Lar.evalStruct(assembly)
	inputfacenumbers = [length(listOfModels[k][2]) for k=1:length(listOfModels)]
	# test input data for containment of reference points
	#-------------------------------------------------------------------------------
	boolmatrix = BitArray(undef, length(innerpoints), length(listOfModels))
    containmenttest = Lar.testinternalpoint2d(listOfModels)
	for (k,point) in enumerate(innerpoints) # k runs on columns
		cells = containmenttest(point) # contents of columns
		#println(k," ",faces)
		for l in cells
			boolmatrix[k,l] = 1
		end
	end
	return W, copEV, copFE, boolmatrix
end

"""
	boolops( assembly::Lar.Struct, op::Symbol )

User interface to 2d Boolean ops, where `op` symbol ``in`` {`:|`, `:&`, `:-`}.
Return a pair `(V,EV)`.
# Example
```julia
assembly = Lar.Struct([wall, openings])
V,EV = boolops(assembly, :-)
```
"""
function boolops(assembly::Lar.Struct, op::Symbol)
    W, copEV, copFE, boolmatrix = Lar.bool2d(assembly)
    boolvars = [boolmatrix[:,k] for k=1:size(boolmatrix,2)]
    if eval(op) == -
        solution = boolvars[1] .& .!(.|(boolvars[2:end]...))
    else
        operator = eval(op)
        solution = operator.(boolvars...)
    end
    chain1d = sparse(copFE') * Int.(solution)
    EV = Lar.cop2lar(copEV)
    EVop = [ev for (k,ev) in enumerate(EV) if abs(chain1d[k])==1 ]
    V = convert(Lar.Points,W')
    return V,EVop
end


text:
module PrettyNumber

using Printf
export prettify, prettyprint

############################
#                          #
#      Default options     #
#                          #
############################

const max_int_len = 5
const max_whole_part_len = 3
const max_decimals = 3

############################
#                          #
#         Unicode          #
#                          #
############################

const supers = [  '⁰'; '¹'; '²'
                ; '³'; '⁴'; '⁵'
                ; '⁶'; '⁷'; '⁸'
                ; '⁹']

const superminus = '⁻'

############################
#                          #
#        Functions         #
#                          #
############################

exponent(x::Number) = floor(Int64,log10(abs(x)))
mantissa(x::Number) = x / exp10(exponent(x))

# Superindex
"Return the unicode version of x as a superindex."
function integer_to_sup(x::Integer)
    str = string(abs(x))
    L = length(str)
    sups = Array{Char}(undef, L)
    for i in 1:L
        c = str[i]
        n = Int64(c) - 48
        sups[i] = supers[n+1]
    end
    result = String(sups)
    if x > 0
        return result
    else
        return "$superminus$result"
    end
end

"""
Pretty print in scientific notation, with `l` decimal places. if `hideexp`,
the exponent will not be shown if it's 0.
Returns a `String`.
"""
function scientific(x::Number; l=:full, hideexp=true)
    m = mantissa(abs(x)) |> string
    # Sign will be added later.
    if l == :full
        mstr = m[1:end]
    elseif l == 0
        mstr = m[1]
    elseif l ≥ 1
        # - 2 because the first number and the point.
        if l > length(m) - 2 # needs padding
            pad = fill('0',l-length(m)+2) |> join
            mstr = string(m,pad)
        else
            mstr = string(m)[1:l+2]
        end
    else
        error("l ∉ [0,∞). Use :full for all places")
    end
    expstr = x |> exponent |> integer_to_sup
    if x < 0
        mstr = string('-',mstr)
    end
    if hideexp && exponent(x) == 0
        return string(mstr)
    else
        return string(mstr,"⋅10",expstr)
    end
end

"""
Number of significant digits.
1.420 has 4, for example, and 0.0540 has 3.
"""
function sig_digits(x::AbstractString)
    # remove sign and exponent
    x = replace(x,r"e.*"=>"")
    x = replace(x,"+"=>"")
    x = replace(x,"-"=>"")
    # remove leading zeroes and possible dot
    x = replace(x,r"^0*\.?0*"=>"")
    # remove point so we dont count it as character
    x = replace(x,"."=>"")
    return length(x)
end
sig_digits(x::Number) = x |> string |> sig_digits

function pretty(x::Number)
    if isinteger(x) && exponent(x) < max_int_len
        char_size = exponent(x) + 1
        if char_size ≤ max_int_len
            return @sprintf("%d",x)
        else # Scientific notation
            return scientific(x,l=0)
        end
    else # It's a float
        whole = floor(x)
        decimal = x - floor(x)
        whole_part_len = @sprintf("%d",abs(whole)) |> length
        decimal_part_len = sig_digits(x) - whole_part_len
        # Small enough, don't touch it.
        if whole_part_len ≤ max_whole_part_len && round(x,digits=max_decimals) ≠ 0
            if decimal_part_len ≤ max_decimals
                return string(x)
            else #truncate decimals
                return round(x,digits=max_decimals) |> string
            end
        elseif x / 10.0^exponent(x) ≈ round(mantissa(x))
            # Just an "int". Does not need decimals.
            # An example is 5e-5, is just 5⋅10⁻⁵.
            superindex = x |> exponent |> integer_to_sup
            whole = round(x / 10.0^exponent(x))
            return @sprintf("%d⋅10%s",mantissa(x) ,superindex)
        else # Too big and not an pseudoint. To scientific.
            dec_places = min(sig_digits(x)-1, max_decimals)
            return scientific(x,l=dec_places)
        end
    end
end

"""
    prettify(str)

Change all the numbers in str by nice looking ones.
"""
function prettify(str::AbstractString)
    # Search for numbers in the string, and
    # substitute them.
    num_r = r"[+-]?\d{1,}\.?\d*e?[+-]?\d*"
    # Regex explanation:
    # ±?  num(1,∞)  .?  num(0,∞)  e? ±?  num(0,∞)
    # +     2       .   4254      e  -   20
    # +2.4254e-20 (e.g.)
    result = str # will be overwriten
    for match_result in eachmatch(num_r,str)
        numstr = match_result.match
        value = parse(Float64,numstr)
        result = replace(result, numstr => pretty(value))
    end
    return result
end

prettyprint(str) = str |> prettify |> println

end # module


text:
# Various mathematical operations
using Distributed,LinearAlgebra

#=
&,And : Logical & of A and B
=#

(&)(A::Assoc,B::Assoc) = and(A::Assoc,B::Assoc)
function and(A::Assoc, B::Assoc)
    #First, create the row and col of the intersection
    ABrow = intersect(A.row,B.row)
    ABcol = intersect(A.col,B.col)
    #Filling the sparse matrix with

    AA = spzeros(size(ABrow,1), size(ABcol,1))
    rowMapping = searchsortedmapping(ABrow,A.row)
    colMapping = searchsortedmapping(ABcol,A.col)
    AA = LinearAlgebra.fillstored!(copy(A.A[rowMapping,colMapping]),1)
    AA = round.(Int64,AA)

    BB = spzeros(size(ABrow,1), size(ABcol,1))
    rowMapping = searchsortedmapping(ABrow,B.row)
    colMapping = searchsortedmapping(ABcol,B.col)
    BB = LinearAlgebra.fillstored!(copy(B.A[rowMapping,colMapping]),1)
    BB = round.(Int64,BB)

    ABA = AA .& BB
    ABA = ABA * 1.0

    return Assoc(ABrow,ABcol,promote([1.0],A.val)[1],ABA)
end

#=
 +,plus : matrix add for Assoc.

 If Assoc has String in value mapping, the Assoc will be logical instead.
=#

+(A::Assoc,B::Assoc) = plus(A::Assoc,B::Assoc)
function plus(A::Assoc,B::Assoc)
    # Check if A or B is empty
    if isempty(A)
       return B;
    end
    if isempty(B)
       return A;
    end

    #Check A,B, if string => Logical
    At = A
    Bt = B

    # There can be numberic values in the "vals" this will not work in this case
    if(A.val != [1.0])
        At = logical(A)
    end
    if(A.val != [1.0])
        Bt = logical(B)
    end

    ## A+B operation
    ABrow = sortedunion(At.row,Bt.row)
    ABcol = sortedunion(At.col,Bt.col)
    if (size(ABrow,1) == 0 || size(ABcol,1) == 0)
        return Assoc([1],[1],0,(+))
    end

    ABA =  spzeros(length(ABrow),length(ABcol))
    Arow = searchsortedmapping(At.row,ABrow)
    Acol = searchsortedmapping(At.col,ABcol)
    Brow = searchsortedmapping(Bt.row,ABrow)
    Bcol = searchsortedmapping(Bt.col,ABcol)
    ABA[Arow,Acol] += At.A
    ABA[Brow,Bcol] += Bt.A
    AB = Assoc(ABrow,ABcol, Array{Union{AbstractString,Number}}([1.0]), ABA)

    return condense(AB) #Incase of negation
end

# Minus for Associative Arrays
-(A::Assoc,B::Assoc) = minus(A::Assoc,B::Assoc)
function minus(A::Assoc,B::Assoc)
    #Check A,B, if string => Logical
    At = A
    Bt = B
    if(A.val != [1.0])
        At = logical(A)
    end
    if(B.val !=[1.0])
        Bt = logical(B)
    end
    ## A*B operation
    ABrow = sortedunion(At.row,Bt.row)
    ABcol = sortedunion(At.col,Bt.col)
    if (size(ABrow,1) == 0 || size(ABcol,1) == 0)
        return Assoc([1],[1],0,(+))
    end

    ABA =  spzeros(length(ABrow),length(ABcol))
    Arow = searchsortedmapping(At.row,ABrow)
    Acol = searchsortedmapping(At.col,ABcol)
    Brow = searchsortedmapping(Bt.row,ABrow)
    Bcol = searchsortedmapping(Bt.col,ABcol)
    ABA[Arow,Acol] += At.A
    ABA[Brow,Bcol] -= Bt.A
    AB = Assoc(ABrow,ABcol, Array{Union{AbstractString,Number}}([1.0]), ABA)

    return condense(AB) #Incase of negation
end

#=
* : matrix multiply between two Assoc.
If there are Strings in Assoc Value mapping, the Assoc will be reduced to logical.
=#

function *(A::Assoc,B::Assoc)
    #Check A,B, if string => Logical
    At = A
    Bt = B
    if(!isa(A.val[1],Number))
        At = logical(A)
    end
    if(!isa(B.val[1],Number))
        Bt = logical(B)
    end
    ## A*B operation
    #Improve intersect with unique sorted sets
    ABintersect = sortedintersect(At.col,Bt.row)

#    if (size(ABintersect,1) == 0)
#        return Assoc([1],[1],0,(+))
#    end

#    AintMap = searchsortedmapping(ABintersect,At.col)
#    BintMap = searchsortedmapping(ABintersect,Bt.row)

    Aref = @spawn searchsortedmapping(ABintersect,At.col)
    Bref = @spawn searchsortedmapping(ABintersect,Bt.row)
    AintMap = fetch(Aref)
    BintMap = fetch(Bref)
#    AintMap, BintMap = sortedintersectmapping(At.col,Bt.row)

#    AintMap,BintMap = map(x -> searchsortedmapping(x[1],x[2]),[(ABintersect,At.col) (ABintersect, Bt.row)])
    Aref = @spawn At.A[:,AintMap]
    Bref = @spawn Bt.A[BintMap,:]
    AA = fetch(Aref)
    BB = fetch(Bref)

    ABA = AA*BB

    AB = Assoc(At.row,Bt.col,Array{Union{AbstractString,Number}}([1.0]),ABA)

    return condense(AB)
end

#=
 transpose : return the transpose of Given Assoc
=#

function transpose(A::Assoc)
    return Assoc(A.col,A.row,A.val,Transpose(A.A))
end

function adjoint(A::Assoc)
    return Assoc(A.col,A.row,A.val,Adjoint(A.A))
end

#abs: Absolutee value of an inputted numerical associative array.
function abs(A::Assoc)

    if isa(A.val[1],AbstractString)
        error("Attempt to Abs a string valued Associative Array, only numerical Associative Array is accepted")
    end

    AT = deepcopy(A)
    AT.A = abs(AT.A)
    return AT
end

# Same as A'*A, perhaps faster
function sqIn(A::Assoc)
    AtA = deepcopy(A)
    if ! isa(A.val[1], Number)
        AtA = logical(A)
    end

    AA = getadj(AtA)
    AAtAA = AA' * AA

    #=
    AtA.A = AAtAA
    AtA.row = AtA.col
    return AtA
=#
    return Assoc(copy(AtA.col),copy(AtA.col),copy(AtA.val),AAtAA)
end

# Same as A*A', perhaps faster
function sqOut(A::Assoc)
    AtA = deepcopy(A)
    if ! isa(A.val[1], Number)
        AtA = logical(A)
    end

    AA = getadj(AtA)
    AAtAA = AA * AA';

#=
    AtA.A = AAtAA;
    AtA.col = AtA.row
    return AtA
=#
    return Assoc(copy(AtA.row),copy(AtA.row),copy(AtA.val),AAtAA)
end

function sum(A::Assoc, dim::Int64)

    if(A.val != [1.0])
        A = logical(A)
    end

    if(dim == 1)
        return condense(Assoc(promote([1.0],A.row)[1],A.col, promote([1.0],A.val)[1],sparse(sum(A.A, dims = 1))))
    elseif (dim == 2)

        return condense(Assoc(A.row,promote([1.0],A.col)[1], promote([1.0],A.val)[1],sparse(sum(A.A, dims = 2))))
    end
end

function sum(A::Assoc)

    if(A.val != [1.0])
        A = logical(A)
    end

    sum(A.A)

end

# Matrix multiply, but resulting values are inner dimension keys concatenated together
function CatKeyMul(A::Assoc,B::Assoc)
    if isa(getcol(A)[1],AbstractString) && isa(getrow(B)[1],AbstractString)
        AB = sortedintersect(A.col,B.row)
        A = A[:,AB]
        B = B[AB,:]
        rrr,ccc,vvv = findnz(getadj(A*B))
        ABVal = Array{Union{AbstractString,Number}}(undef,length(rrr))
        for i in 1:length(rrr)
            r = rrr[i]
            c = ccc[i]
            ABvalList = sortedintersect(getcol(A[r,:]),getrow(B[:,c]))
            if length(ABvalList) > 0
                val = join(ABvalList,";")*";"
                ABVal[i] = val
            end
        end
        return Assoc(getrow(A)[rrr],getcol(B)[ccc],ABVal)
    else
        return A*B
    end

end

# Matrix multiply, but resulting values are previous values concatenated together
function CatValMul(A::Assoc,B::Assoc)
    if isa(getval(A)[1],AbstractString) && isa(getval(B)[1],AbstractString)
        AB = sortedintersect(A.col,B.row)
        A = A[:,AB]
        B = B[AB,:]
        rrr,ccc,vvv = findnz(getadj(A*B))
        ABVal = Array{Union{AbstractString,Number},1}(undef, length(rrr))
        for i in 1:length(rrr)
            r = rrr[i]
            c = ccc[i]
            ABIntersect = sortedintersect(getcol(A[r,:]),getrow(B[:,c]))
            ABValList = Array{Union{AbstractString,Number},1}(undef, 0)
            for x in ABIntersect
                push!(ABValList,getval(A[r,x*","])[1])
                push!(ABValList,getval(B[x*",",c])[1])
            end
            if length(ABValList) > 0
                val = join(ABValList,";")*";"
                ABVal[i] = val
            end
        end
        return Assoc(getrow(A)[rrr],getcol(B)[ccc],ABVal)
    else
        return A*B
    end

end

#=
OutDegree : Calculate the out-degree distribution of graph A
InDegree : Calculate the in-degree distribution of graph A
=#

using SparseArrays

function OutDegree(A)
    if isa(A,Assoc)
        A = A.A
    end
    dout = sum(A,dims = 2)
    dout_i = getindex.(findall(!iszero,dout),1)
    dout_v = dout[dout_i]
    ndout = sum(sparse(dout_i,dout_v,1),dims = 1)
    return ndout
end

function InDegree(A)
    if isa(A,Assoc)
        A = A.A
    end
    din = sum(A,dims = 1)
    din_i = getindex.(findall(!iszero,din),1)
    din_v = din[din_i]
    ndin = sum(sparse(din_i,din_v,1),dims = 1)
    return ndin
end


text:
using sem, Arrow, ModelingToolkit, LinearAlgebra, SparseArrays, DataFrames, Optim, LineSearches

cd("./test")
## Observed Data
dat = DataFrame(Arrow.Table("comparisons/reg_1.arrow"))
par = DataFrame(Arrow.Table("comparisons/reg_1_par.arrow"))
start_lav = DataFrame(Arrow.Table("comparisons/reg_1_start.arrow"))

# observed
semobserved = SemObsCommon(data = Matrix{Float64}(dat))

## Model definition
@ModelingToolkit.variables x[1:243]

#F
Ind = collect(1:25)
Jnd = collect(1:25)
V = fill(1,25)
F = sparse(Ind, Jnd, V, 25, 46)

#A
Ind = [collect(2:25); 5; collect(6:25); fill(27, 19); 1; 4]
Jnd = [fill(26, 24); 27; collect(27:46); collect(28:46); 26; 27]
V =[x[1:64]; 1.0; 1.0]
A = sparse(Ind, Jnd, V, 46, 46)

#S
Ind = [collect(1:6); 26; 27; collect(28:46)]
Jnd = [collect(1:6); 26; 27; collect(28:46)]
V = [x[65:72]; fill(1.0, 19)]
S = sparse(Ind, Jnd, V, 46, 46)
xind = 73
for i = 28:46
    for j = (i+1):46
       S[i,j] = x[xind]
       S[j,i] = x[xind]
       global xind += 1
    end
end

par_order = [
    collect(2:6);
    collect(11:29);
    8;9;
    collect(49:67);
    collect(106:124);
    collect(125:130);
    131; 132;
    collect(133:303)]

start_val = Vector{Float64}(start_lav.est[par_order])

# loss
loss = Loss([SemML(semobserved, [0.0], similar(start_val))])

# imply
imply = ImplySymbolic(A, S, F, x, start_val)

# diff
grad_ml = sem.∇SemML(A, S, F, x, start_val)

diff_ana =
    SemAnalyticDiff(
        LBFGS(
            m = 50,
            alphaguess = InitialHagerZhang(),
            linesearch = HagerZhang()),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_ml,))

model_ana = Sem(semobserved, imply, loss, diff_ana)

using BenchmarkTools

@benchmark solution_ana = sem_fit(model_ana)

all(abs.(solution_ana.minimizer .- par.est[par_order]) .< 0.05*abs.(par.est[par_order]))