text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    SDFileBond, SmilesBond, SmartsBond,
    setorder, setnotation, setstereo

struct SDFileBond <: Bond
    """Bond
    * Notation
        * Single bond
            * 0: u - v
            * 1: u ◀ v (Up-arrow)
            * 4: u ~ v (Up or down)
            * 6: u ◁ v (Down-arrow)
        * Double bond
            * 0: v ニ u (clockwise, default)
            * 1: u ニ v (counter-clockwise)
            * 2: u ＝ v (equal length, for terminal bond by default)
            * 3: u × v (Cis-Trans Unknown)
    """
    order::Int
    notation::Int
    stereo::Symbol
end

SDFileBond() = SDFileBond(1, 0, :unspecified)
SDFileBond(order) = SDFileBond(order, 0, :unspecified)
SDFileBond(order, notation) = SDFileBond(order, notation, :unspecified)

setorder(edge::SDFileBond, order
    ) = SDFileBond(order, edge.notation, edge.stereo)

setnotation(edge::SDFileBond, notation
    ) = SDFileBond(edge.order, notation, edge.stereo)

setstereo(edge::SDFileBond, cistrans
    ) = SDFileBond(edge.order, edge.notation, cistrans)

struct SmilesBond <: Bond
    order::Int
    isaromatic::Bool
    direction::Symbol
    stereo::Symbol
end

SmilesBond() = SmilesBond(1, false, :unspecified, :unspecified)
SmilesBond(order) = SmilesBond(order, false, :unspecified, :unspecified)
SmilesBond(order, isaromatic, direction
    ) = SmilesBond(order, isaromatic, direction, :unspecified)

setorder(edge::SmilesBond, order
    ) = SmilesBond(order, edge.isaromatic, edge.direction, edge.stereo)

setstereo(edge::SmilesBond, cistrans
    ) = SmilesBond(edge.order, edge.isaromatic, edge.direction, cistrans)

struct SmartsBond <: QueryBond
    query::Pair
end

SmartsBond() = SmartsBond(:any => true)


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

mutable struct IoK8sApiCoreV1Toleration <: SwaggerModel
    effect::Any # spec type: Union{ Nothing, String } # spec name: effect
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    operator::Any # spec type: Union{ Nothing, String } # spec name: operator
    tolerationSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: tolerationSeconds
    value::Any # spec type: Union{ Nothing, String } # spec name: value

    function IoK8sApiCoreV1Toleration(;effect=nothing, key=nothing, operator=nothing, tolerationSeconds=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Toleration, Symbol("effect"), effect)
        setfield!(o, Symbol("effect"), effect)
        validate_property(IoK8sApiCoreV1Toleration, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApiCoreV1Toleration, Symbol("operator"), operator)
        setfield!(o, Symbol("operator"), operator)
        validate_property(IoK8sApiCoreV1Toleration, Symbol("tolerationSeconds"), tolerationSeconds)
        setfield!(o, Symbol("tolerationSeconds"), tolerationSeconds)
        validate_property(IoK8sApiCoreV1Toleration, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiCoreV1Toleration

const _property_map_IoK8sApiCoreV1Toleration = Dict{Symbol,Symbol}(Symbol("effect")=>Symbol("effect"), Symbol("key")=>Symbol("key"), Symbol("operator")=>Symbol("operator"), Symbol("tolerationSeconds")=>Symbol("tolerationSeconds"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiCoreV1Toleration = Dict{Symbol,String}(Symbol("effect")=>"String", Symbol("key")=>"String", Symbol("operator")=>"String", Symbol("tolerationSeconds")=>"Int64", Symbol("value")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1Toleration }) = collect(keys(_property_map_IoK8sApiCoreV1Toleration))
Swagger.property_type(::Type{ IoK8sApiCoreV1Toleration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Toleration[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Toleration }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Toleration[property_name]

function check_required(o::IoK8sApiCoreV1Toleration)
    true
end

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

ext: jl


text:
include("CPP.jl")
include("DCPP.jl")
include("MCPP.jl")
include("MCARP.jl")

function starts_in(arc, vertex, arc_set)
    return(arc_set[arc, :start_node]==vertex)
end

function ends_in(arc, vertex, arc_set)
    return(arc_set[arc, :end_node]==vertex)
end

function find_problem_type(nb_vehicles, graph)
    """Finds of wich type the problem is."""

    if nb_vehicles == 1     # Chinese postman problem type
        if (nrow(graph["req_arcs"])==0) & (nrow(graph["noreq_arcs"])==0)
            problem_type = "CPP"
        elseif (nrow(graph["req_edges"])==0) & (nrow(graph["noreq_edges"])==0)
            problem_type = "DCPP"
        else
            problem_type = "MCPP"
        end
    elseif nb_vehicles > 1  # MCARP problem type
        problem_type = "MCARP"
    else
        error("Problem with vehicle number. Check the dataset's validity.")
    end

    return(problem_type)
end

function solve_with_appropriate_model(problem_type, data, graph)
    """Solve problem with the apropriated modelisation."""

    if problem_type == "MCARP"
        (tours, solving_time) = MCARP_solve(data, graph)
    elseif problem_type == "CPP"
        (tours, solving_time) = CPP_solve(data, graph)
    elseif problem_type == "DCPP"
        (tours, solving_time) = DCPP_solve(data, graph)
    elseif problem_type == "MCPP"
        (tours, solving_time) = MCPP_solve(data, graph)
    else
        error("Problem while choosing wich model to use.")
    end

    return(tours, solving_time)
end

function solve_single_depot(data, graph)
    """Solve problem in the case of a single-depot dataset."""

    problem_type = find_problem_type(data[1, :nb_vehicles], graph)
    println("Algorith used : ", problem_type)
    println("Solving problem ...")
    (tours, solving_time) = solve_with_appropriate_model(problem_type, data, graph)

    return(tours, solving_time)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
########################
# multi-threading task
########################

@inline function subspace_mul_generic_task!(S::AbstractVector, indices, U::AbstractMatrix, subspace, range, offset)
    D = StrideArrays.static_length(indices)
    y = StrideArray{T}(undef, (D, ))
    if D == 4
        for s in range
            k = subspace[s]
            subspace_mul_kernel_generic_4x4!(S, y, indices, U, k, offset)
        end
    else
        for s in range
            k = subspace[s]
            subspace_mul_kernel_generic!(S, y, indices, U, k, offset)
        end
    end
    return
end

@inline function subspace_mul_generic_task!(S::AbstractMatrix, indices, U::AbstractMatrix, subspace, range, offset)
    D = StrideArrays.static_length(indices)
    stride = size(S, 1)
    y = StrideArray{T}(undef, (D, ))
    if D == 4
        for (b, s) in range
            k = subspace[s]
            subspace_mul_kernel_generic_4x4!(S, y, indices, U, k, b, offset)
        end
    else
        for (b, s) in range
            k = subspace[s]
            subspace_mul_kernel_generic!(S, y, indices, U, k, b, offset)
        end
    end
    return
end

@inline function subspace_mul_complex_task!(Sr::AbstractMatrix{T}, indices, U, subspace, range::CartesianPartition, offset) where {T <: Base.HWReal}
    D = StrideArrays.static_length(indices)
    y = StrideArray{T}(undef, (D, )), StrideArray{T}(undef, (D, ))

    for s in range
        k = subspace[s]
        subspace_mul_kernel_complex!(Sr, y, indices, U, k, offset)
    end
    return
end

@inline function subspace_mul_complex_task!(Sr::AbstractArray{T, 3}, indices, U, subspace, range::CartesianPartition, offset) where {T <: Base.HWReal}
    D = ArrayInterface.static_length(indices)
    C = StrideArray{T}(undef, (D, StaticInt{8}())), StrideArray{T}(undef, (D, StaticInt{8}()))
    Bmax = size(Sr, 2)

    for (_b, s) in range
        k = subspace[s]
        subspace_mul_kernel_complex!(Sr, C, indices, U, k, _b, Bmax, offset)
    end
    return
end


text:
module SKCrossValidators

using PyCall

# standard included modules
using DataFrames
using Random
using ..AbsTypes
using ..Utils

import ..CrossValidators: crossvalidate
export crossvalidate

const metric_dict = Dict{String,PyObject}()
const SKM = PyNULL()

function __init__()
   copy!(SKM, pyimport_conda("sklearn.metrics","scikit-learn"))

   metric_dict["roc_auc_score"]                   = SKM.roc_auc_score
   metric_dict["accuracy_score"]                  = SKM.accuracy_score
   metric_dict["auc"]                             = SKM.auc
   metric_dict["average_precision_score"]         = SKM.average_precision_score
   metric_dict["balanced_accuracy_score"]         = SKM.balanced_accuracy_score
   metric_dict["brier_score_loss"]                = SKM.brier_score_loss
   metric_dict["classification_report"]           = SKM.classification_report
   metric_dict["cohen_kappa_score"]               = SKM.cohen_kappa_score
   metric_dict["confusion_matrix"]                = SKM.confusion_matrix
   metric_dict["f1_score"]                        = SKM.f1_score
   metric_dict["fbeta_score"]                     = SKM.fbeta_score
   metric_dict["hamming_loss"]                    = SKM.hamming_loss
   metric_dict["hinge_loss"]                      = SKM.hinge_loss
   metric_dict["log_loss"]                        = SKM.log_loss
   metric_dict["matthews_corrcoef"]               = SKM.matthews_corrcoef
   metric_dict["multilabel_confusion_matrix"]     = SKM.multilabel_confusion_matrix
   metric_dict["precision_recall_curve"]          = SKM.precision_recall_curve
   metric_dict["precision_recall_fscore_support"] = SKM.precision_recall_fscore_support
   metric_dict["precision_score"]                 = SKM.precision_score
   metric_dict["recall_score"]                    = SKM.recall_score
   metric_dict["roc_auc_score"]                   = SKM.roc_auc_score
   metric_dict["roc_curve"]                       = SKM.roc_curve
   metric_dict["jaccard_score"]                   = SKM.jaccard_score
   metric_dict["zero_one_loss"]                   = SKM.zero_one_loss
   # regression
   metric_dict["mean_squared_error"]              = SKM.mean_squared_error
   metric_dict["mean_squared_log_error"]          = SKM.mean_squared_log_error
   metric_dict["mean_absolute_error"]             = SKM.mean_absolute_error
   metric_dict["median_absolute_error"]           = SKM.median_absolute_error
   metric_dict["r2_score"]                        = SKM.r2_score
   metric_dict["max_error"]                       = SKM.max_error
   metric_dict["mean_poisson_deviance"]           = SKM.mean_poisson_deviance
   metric_dict["mean_gamma_deviance"]             = SKM.mean_gamma_deviance
   metric_dict["mean_tweedie_deviance"]           = SKM.mean_tweedie_deviance
   metric_dict["explained_variance_score"]        = SKM.explained_variance_score
   # custom
   metric_dict["make_scorer"]                     = SKM.make_scorer
end

function checkfun(sfunc::String)
    if !(sfunc in keys(metric_dict))
        println("$sfunc metric is not supported")
        println("metric: ",keys(metric_dict))
        error("Metric keyword error")
    end
end

"""
    crossvalidate(pl::Machine,X::DataFrame,Y::Vector,sfunc::String="balanced_accuracy_score",nfolds=10)

Runs K-fold cross-validation using balanced accuracy as the default. It supports the
following metrics for classification:
- accuracy_score
- balanced_accuracy_score
- cohen_kappa_score
- jaccard_score
- matthews_corrcoef
- hamming_loss
- zero_one_loss
- f1_score
- precision_score
- recall_score

and the following metrics for regression:
- mean_squared_error
- mean_squared_log_error
- median_absolute_error
- r2_score
- max_error
- explained_variance_score

Additionally, it supports custom metric functions by supplying a function that
meets the requirements of sklearn's make_scorer function in the place of the
supplied metrics.
"""
function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
                       sfunc::String; nfolds=10, verbose::Bool=true)
    checkfun(sfunc)
    pfunc = metric_dict[sfunc]
    metric(a,b) = pfunc(a,b)
    crossvalidate(pl,X,Y,metric,nfolds,verbose)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,sfunc::String,folds::Int)
   crossvalidate(pl,X,Y,sfunc,nfolds=folds)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,sfunc::String,verby::Bool)
   crossvalidate(pl,X,Y,sfunc,verbose=verby)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
                       sfunc::String, folds::Int,verby::Bool)
   crossvalidate(pl,X,Y,sfunc,nfolds=folds,verbose=verby)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
                       sfunc::String,averagetype::String,nfolds=10,verbose::Bool=true)
    checkfun(sfunc)
    pfunc = metric_dict[sfunc]
    metric(a,b) = pfunc(a,b,average=averagetype)
    crossvalidate(pl,X,Y,metric,nfolds,verbose)
end

function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
    sfunc::Function; nfolds=10, verbose::Bool=true)
    pfunc = metric_dict["make_scorer"]
    pfunc = pfunc(sfunc, greater_is_better=true)
    metric(a,b) = pfunc(a,b)
    crossvalidate(pl,X,Y,metric,nfolds,verbose)
end

end


text:
__precompile__()
module Requests

export URI, FileParam, headers, cookies, statuscode, post, requestfor
export get_streaming, post_streaming, write_chunked
export save
export set_proxy, set_https_proxy, get_request_settings

import Base: get, write

import Base.Filesystem: File

import URIParser: URI
import HttpCommon: Cookie

using Compat
using HttpParser
using HttpCommon
using URIParser
using MbedTLS
using Codecs
using JSON
using Libz

const CRLF = "\r\n"

include("parsing.jl")
include("multipart.jl")
include("streaming.jl")
include("mimetypes.jl")

function __init__()
    __init_parsing__()
    __init_streaming__()
    init_proxy()
end

type GlobalSettings
    http_proxy::Nullable{URI}
    https_proxy::Nullable{URI}
    max_redirects::Int
end
GlobalSettings() = GlobalSettings(Nullable{URI}(), Nullable{URI}(), 5)

function Base.show(io::IO, settings::GlobalSettings)
    println(io, "HTTP proxy: ", isnull(settings.http_proxy) ? "No proxy" : get(settings.http_proxy))
    println(io, "HTTPS proxy: ", isnull(settings.https_proxy) ? "No proxy" : get(settings.https_proxy))
    print(io, "Max redirects: ", settings.max_redirects)
end

const SETTINGS = GlobalSettings()

function init_proxy()
    if haskey(ENV, "http_proxy")
        try
            proxy = ENV["http_proxy"]
            set_proxy(proxy)
            info("Using proxy $proxy from http_proxy environment variable")
        catch err
            warn("Problem parsing http_proxy environment variable $(ENV["http_proxy"]); ignoring\n $err")
        end
    end
    if haskey(ENV, "https_proxy")
        try
            proxy = ENV["https_proxy"]
            set_https_proxy(proxy)
            info("Using HTTPS proxy $proxy from https_proxy environment variable")
        catch err
            warn("Problem parsing https_proxy environment variable $(ENV["https_proxy"]); ignoring\n $err")
        end
    end
end

set_proxy(proxy::URI) = set_proxy(Nullable(proxy))
set_proxy(proxy::Nullable{URI}) = SETTINGS.http_proxy = proxy
set_proxy(proxy) = set_proxy(URI(proxy))
set_https_proxy(proxy) = SETTINGS.https_proxy=Nullable(URI(proxy))

get_request_settings() = SETTINGS

## Convenience methods for extracting the payload of a response
for kind in [:Response, :Request]
    @eval bytes(r::$kind) = r.data
    @eval text(r::$kind) = Compat.String(bytes(r))
    @eval Base.String(r::$kind) = text(r)
    if VERSION < v"0.6.0-dev.2521"
        @eval function Base.bytestring(r::$kind)
            Base.depwarn("bytestring(r::$($kind)) is deprecated, use String(r) instead.", :bytestring)
            String(r)
        end
    end
    @eval Base.readstring(r::$kind) = text(r)
    if VERSION < v"0.6.0-dev.2521"
        @eval function Base.readall(r::$kind)
            Base.depwarn("readall(r::$($kind)) is deprecated, use readstring(r) instead.", :readall)
            readstring(r)
        end
    end
    @eval Base.read(r::$kind) = bytes(r)
    @eval json(r::$kind; kwargs...) = JSON.parse(text(r); kwargs...)

    ## Response getters to future-proof against changes to the Response type
    @eval headers(r::$kind) = r.headers
end

cookies(r::Response) = r.cookies
statuscode(r::Response) = r.status

function requestfor(r::Response)
    isnull(r.request) && error("No associated request for response")
    get(r.request)
end

history(r::Response) = r.history

# Stolen from https://github.com/dcjones/Gadfly.jl/blob/7fd56991e55b6617d37d7e3d0d69a310bdd36b05/src/Gadfly.jl#L1016
function open_file(filename)
    if is_apple()
        run(`open $(filename)`)
    elseif is_unix()
        run(`xdg-open $(filename)`)
    elseif is_windows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    end
end

function mimetype(r::Response)
    if haskey(headers(r), "Content-Type")
        ct = split(headers(r)["Content-Type"], ";")[1]
        return Nullable(ct)
    else
        return Nullable{String}()
    end
end

function contentdisposition(r::Response)
    if haskey(headers(r), "Content-Disposition")
        cd = split(headers(r)["Content-Disposition"], ";")
        if length(cd) ≥ 2
            filepart = split(cd[2], "=", limit=2)
            if length(filepart) == 2
                return  filepart[2] |> strip |> f->strip(f, '"') |> Nullable
            end
        end
    end
    return Nullable{String}()
end

"""
`save(r::Response, path=".")`

Saves the data in the response in the directory `path`. If the path is a directory,
then the filename is automatically chosen based on the response headers.

Returns the full pathname of the saved file.
"""
function save(r::Response, path=".")
    if !isdir(path)
        filename = path
    else
        maybe_basename = contentdisposition(r)
        if !isnull(maybe_basename)
            filename = joinpath(path, get(maybe_basename))
        else
            ext = "txt"
            maybe_mt = mimetype(r)
            if !isnull(maybe_mt)
                mt = get(maybe_mt)
                if haskey(MIMETYPES, mt)
                    ext = MIMETYPES[mt]
                else
                    if '/' ∉ mt
                        ext = mt
                    end
                end
            end
            basefile = Dates.format(now(), "mm-dd-yyyy-H-M")
            filename = joinpath(path, "$basefile.$ext")
        end
    end
    open(filename, "w") do file
        write(file, bytes(r))
    end
    filename
end

"""
`view(r::Response)`

View the data in the response with whatever application is associated with
its mimetype.
"""
function view(r::Response)
    path = save(r, mktempdir())
    open_file(path)
end

http_port(uri) = uri.port == 0 ? 80 : uri.port
https_port(uri) = uri.port == 0 ? 443 : uri.port

function default_request(method,resource,host,data,user_headers=Dict{Union{},Union{}}())
    headers = Dict(
        "User-Agent" => "Requests.jl/0.0.0",
        "Host" => host,
        "Accept" => "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8"
        )
    if !isempty(data)
        headers["Content-Length"] = dec(sizeof(data))
    end
    merge!(headers,user_headers)
    Request(method,resource,headers,data)
end

function resourcefor(uri::URI)
    r = uri.path
    if !isempty(uri.query)
        r = "$r?$(uri.query)"
    end
    r
end

function default_request(uri::URI,headers,data,method)
    resource = resourcefor(uri)
    if !isempty(uri.userinfo) && !haskey(headers,"Authorization")
        headers["Authorization"] = "Basic $(base64encode(uri.userinfo))"
    end
    host = uri.port == 0 ? uri.host : "$(uri.host):$(uri.port)"
    request = default_request(method,resource,host,data,headers)
    request.uri = uri
    return request
end

scheme(uri::URI) = isdefined(uri, :scheme) ? uri.scheme : uri.schema

function format_query_str(queryparams; uri = URI(""))
    query_str = isempty(uri.query) ? string() : string(uri.query, "&")

    for (k, v) in queryparams
        if isa(v, Array)
            query_str *= join(map(vi -> "$(URIParser.escape(string(k)))=$(URIParser.escape(string(vi)))", v), "&") * "&"
        else
            query_str *= "$(URIParser.escape(string(k)))=$(URIParser.escape(string(v)))&"
        end
    end
    chop(query_str) # remove the trailing &
end

timeout_in_sec(::Void) = Inf
timeout_in_sec(t::Dates.TimePeriod) = Dates.toms(t)/1000.
timeout_in_sec(t) = convert(Float64, t)

cookie_value(c::Cookie) = c.value
cookie_value(s) = s
function cookie_request_header(d::Dict)
    join(["$key=$(cookie_value(val))" for (key,val) in d], ';')
end
cookie_request_header(cookies::AbstractVector{Cookie}) =
    cookie_request_header(Dict([Pair(cookie.name, cookie.value) for cookie in cookies]))

const is_location = r"^location$"i

function get_redirect_uri(response)
    300 <= statuscode(response) < 400 || return Nullable{URI}()
    hdrs = headers(response)
    for (key, val) in hdrs
        if is_location(key)
            uri = URI(val)
            if isempty(uri.host)  # Redirect URL was given as a relative path
                request = requestfor(response)
                uri = URI(request.uri.scheme, request.uri.host, request.uri.port, uri.path, uri.query)
            end
            return Nullable(uri)
        end
    end
    return Nullable{URI}()
end

immutable RedirectException <: Exception
    max_redirects::Int
end

function Base.show(io::IO, err::RedirectException)
    print(io, "RedirectException: more than $(err.max_redirects) redirects attempted.")
end

macro check_body()
  has_body = esc(:has_body)
  write_body = esc(:write_body)
  quote
    $write_body || error("Incompatible arguments: write_body cannot be false if a data argument is provided.")
    $has_body && error("Multiple body options specified. Please only specify one")
    $has_body = true
  end
end

function do_request(uri::URI, verb; kwargs...)
    response_stream = do_stream_request(uri, verb; kwargs...)
    response = response_stream.response
    response.data = read(response_stream)
    if get(response.headers, "Content-Encoding", "") ∈ ("gzip", "deflate")
        if !isempty(response.data)
            response.data = response.data |> ZlibInflateInputStream |> read
        end
    end
    response
end

parse_request_data(data) = (data, "application/octet-stream")
parse_request_data(data::Associative) =
  (format_query_str(data), "application/x-www-form-urlencoded")

function do_stream_request(uri::URI, verb; headers = Dict{AbstractString, AbstractString}(),
                            cookies = nothing,
                            data = nothing,
                            json = nothing,
                            files = FileParam[],
                            timeout = nothing,
                            query::Dict = Dict(),
                            allow_redirects = true,
                            max_redirects = SETTINGS.max_redirects,
                            history = Response[],
                            tls_conf = TLS_VERIFY,
                            write_body = true,
                            gzip_data = false,
                            compressed = false,
                            proxy = SETTINGS.http_proxy,
                            https_proxy = SETTINGS.https_proxy
                            )

    query_str = format_query_str(query; uri = uri)
    newuri = URI(uri; query = query_str)
    timeout_sec = timeout_in_sec(timeout)

    body = ""
    has_body = false
    if json ≠ nothing
        @check_body
        if get(headers,"Content-Type","application/json") != "application/json"
            error("Tried to send json data with incompatible Content-Type")
        end
        headers["Content-Type"] = "application/json"
        body = JSON.json(json)
    end

    if data ≠ nothing
        @check_body
        body, default_content_type = parse_request_data(data)
        if "Content-Type" ∉ keys(headers)
            headers["Content-Type"] = default_content_type
        end
        if gzip_data
            headers["Content-Encoding"] = "gzip"
            body = Vector{UInt8}(body) |> ZlibDeflateInputStream |> read
        end
    end

    if compressed
        headers["Accept-Encoding"] = "gzip, deflate"
    end

    if cookies ≠ nothing
        headers["Cookie"] = cookie_request_header(cookies)
    end

    request = default_request(newuri, headers, body, verb)
    if isempty(files)
        response_stream = open_stream(request, tls_conf, timeout_sec, proxy, https_proxy)
        if write_body
            write(response_stream, request.data)
        end
    else
        @check_body
        verb == "POST" || error("Multipart file post only supported with POST")
        if haskey(headers,"Content-Type") && !beginswith(headers["Content-Type"],"multipart/form-data")
            error("Tried to send form data with invalid Content-Type. ")
        end
        multipart_settings = prepare_multipart_request!(request, files)
        response_stream = open_stream(request, tls_conf, timeout_sec, proxy, https_proxy)
        send_multipart(response_stream, multipart_settings, files)
    end
    main_task = current_task()
    @schedule begin
        try
            process_response(response_stream)
        catch err
            Base.throwto(main_task, err)
        end
        while response_stream.state < BodyDone
            wait(response_stream)
        end
        if response_stream.state != UpgradeConnection
            close(response_stream)
        end
    end
    if write_body
        while response_stream.state < HeadersDone
            wait(response_stream)
        end
        response_stream.response.history = history
        if allow_redirects && verb ≠ :head
            redirect_uri = get_redirect_uri(response_stream)
            if !isnull(redirect_uri)
                length(response_stream.response.history) > max_redirects &&
                    throw(RedirectException(max_redirects))
                push!(history, response_stream.response)
                return do_stream_request(get(redirect_uri), verb; headers=headers,
                     data=data, json=json, files=files, timeout=timeout,
                     allow_redirects=allow_redirects, max_redirects=max_redirects,
                     history=history, tls_conf=tls_conf, compressed=compressed)
            end
        end

    end
    return response_stream
end

for f in [:get, :post, :put, :delete, :head,
          :trace, :options, :patch, :connect]
    f_str = uppercase(string(f))
    f_stream = Symbol(string(f, "_streaming"))
    @eval begin
        function ($f)(uri::URI, data::AbstractString; headers::Dict=Dict())
            do_request(uri, $f_str; data=data, headers=headers)
        end
        function ($f_stream)(uri::URI, data::AbstractString; headers::Dict=Dict())
            do_stream_request(uri, $f_str; data=data, headers=headers)
        end

        ($f)(uri::AbstractString; args...) = ($f)(URI(uri); args...)
        ($f)(uri::URI; args...) = do_request(uri, $f_str; args...)

        ($f_stream)(uri::AbstractString; args...) = ($f_stream)(URI(uri); args...)
        ($f_stream)(uri::URI; args...) = do_stream_request(uri, $f_str; args...)
    end
end

end


text:
"""
    scaled_LHC_sampling_plan(search_range::Array{Tuple{Float64,Float64},1},num_samples,sampling_plan_opt_gens;trace=false)

Create a Latin Hypercube sampling plan scaled to the `search_range` where
`size(plan) = (num_dimensions,num_samples)`.
"""
function scaled_LHC_sampling_plan(search_range::Array{Tuple{Float64,Float64},1},num_samples,sampling_plan_opt_gens;trace=false)

    #Printing of creation
    trace && println("Creating optimized Latin Hypercube Sampling Plan ... \n")

    unscaled_plan = LHCoptim(num_samples,length(search_range),sampling_plan_opt_gens)[1]'

    plan = Array{eltype(eltype(search_range))}(undef,size(unscaled_plan))
    for (i,sr) in enumerate(search_range)
        plan[i,:] = scale(unscaled_plan[i,:],sr[1],sr[2])
    end
    return plan
end


text:
# Given the row and column indices of the lower left
# vertex, add the location where the contour level
# crosses the specified edge.
function interpolate(x, y, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
    xi, yi = ind
    @inbounds if edge == W
        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
        x_interp = x[xi]
    elseif edge == E
        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
        x_interp = x[xi + 1]
    elseif edge == N
        y_interp = y[yi + 1]
        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
    elseif edge == S
        y_interp = y[yi]
        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
    end

    return VT(x_interp, y_interp)
end

function interpolate(x::AbstractRange, y::AbstractRange, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
    xi, yi = ind
    @inbounds if edge == W
        y_interp = y[yi] + step(y) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
        x_interp = x[xi]
    elseif edge == E
        y_interp = y[yi] + step(y) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
        x_interp = x[xi + 1]
    elseif edge == N
        y_interp = y[yi + 1]
        x_interp = x[xi] + step(x) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
    elseif edge == S
        y_interp = y[yi]
        x_interp = x[xi] + step(x) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
    end

    return VT(x_interp, y_interp)
end

function interpolate(x::AbstractMatrix, y::AbstractMatrix, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
    xi, yi = ind
    @inbounds if edge == W
        Δ = [y[xi,  yi+1] - y[xi,  yi  ], x[xi,  yi+1] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi,  yi+1] - z[xi,  yi  ])
        y_interp = y[xi,yi] + Δ[1]
        x_interp = x[xi,yi] + Δ[2]
    elseif edge == E
        Δ = [y[xi+1,yi+1] - y[xi+1,yi  ], x[xi+1,yi+1] - x[xi+1,yi  ]].*(h - z[xi+1,yi  ])/(z[xi+1,yi+1] - z[xi+1,yi  ])
        y_interp = y[xi+1,yi] + Δ[1]
        x_interp = x[xi+1,yi] + Δ[2]
    elseif edge == N
        Δ = [y[xi+1,yi+1] - y[xi,  yi+1], x[xi+1,yi+1] - x[xi,  yi+1]].*(h - z[xi,  yi+1])/(z[xi+1,yi+1] - z[xi,  yi+1])
        y_interp = y[xi,yi+1] + Δ[1]
        x_interp = x[xi,yi+1] + Δ[2]
    elseif edge == S
        Δ = [y[xi+1,yi  ] - y[xi,  yi  ], x[xi+1,yi  ] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi+1,yi  ] - z[xi,  yi  ])
        y_interp = y[xi,yi] + Δ[1]
        x_interp = x[xi,yi] + Δ[2]
    end

    return VT(x_interp, y_interp)
end


text:
struct ResultSet
    columns::AbstractVector{AbstractVector}
    names::Vector{Symbol}
    description::String
    metadata::String
end

function featherread(filename::AbstractString; use_mmap=true)
    data = loadfile(filename, use_mmap=use_mmap)
    ctable = getctable(data)
    ncols = length(ctable.columns)
    colnames = [Symbol(col.name) for col in ctable.columns]
    coltypes = [juliatype(col) for col in ctable.columns]
    columns = ArrowVector[constructcolumn(coltypes[i], data, ctable.columns[i].metadata, ctable.columns[i]) for i in 1:ncols]
    return ResultSet(columns, colnames, ctable.description, ctable.metadata)
end

#=====================================================================================================
    new column construction stuff
=====================================================================================================#
Base.length(p::Metadata.PrimitiveArray) = p.length

startloc(p::Metadata.PrimitiveArray) = p.offset+1

Arrow.nullcount(p::Metadata.PrimitiveArray) = p.null_count

function bitmasklength(p::Metadata.PrimitiveArray)
    nullcount(p) == 0 ? 0 : padding(bytesforbits(length(p)))
end

function offsetslength(p::Metadata.PrimitiveArray)
    isprimitivetype(p.dtype) ? 0 : padding((length(p)+1)*sizeof(Int32))
end

valueslength(p::Metadata.PrimitiveArray) = p.total_bytes - offsetslength(p) - bitmasklength(p)

function offsetsloc(p::Metadata.PrimitiveArray)
    if isprimitivetype(p.dtype)
        throw(ErrorException("Trying to obtain offset values for primitive array."))
    end
    startloc(p) + bitmasklength(p)
end

# override default offset type
Locate.Offsets(col::Metadata.PrimitiveArray) = Locate.Offsets{Int32}(offsetsloc(col))

Locate.length(col::Metadata.PrimitiveArray) = length(col)
Locate.values(col::Metadata.PrimitiveArray) = startloc(col) + bitmasklength(col) + offsetslength(col)
# this is only relevant for lists, values type must be UInt8
Locate.valueslength(col::Metadata.PrimitiveArray) = valueslength(col)
Locate.bitmask(col::Metadata.PrimitiveArray) = startloc(col)

function constructcolumn(::Type{T}, data::Vector{UInt8}, meta::Metadata.CategoryMetadata, col::Metadata.Column) where T
    reftype = juliatype(col.values.dtype)
    DictEncoding{T}(locate(data, reftype, col.values), locate(data, T, col.metadata.levels))
end

function constructcolumn(::Type{Union{T,Missing}}, data::Vector{UInt8}, meta::Metadata.CategoryMetadata, col::Metadata.Column) where T
    reftype = Union{juliatype(col.values.dtype),Missing}
    DictEncoding{Union{T,Missing}}(locate(data, reftype, col.values), locate(data, T, col.metadata.levels))
end

function constructcolumn(::Type{T}, data::Vector{UInt8}, meta, col::Metadata.Column) where T
    locate(data, T, col.values)
end


text:
function LG.SimpleGraph(g::VStaticGraph{T,N,false}) where {T,N}

end

function LG.SimpleDiGraph(g::VStaticGraph{T,N,true}) where {T,N}

end

function LG.add_edge!(g::VStaticGraph{T,N}, s, d) where {T,N}
    LG.has_vertex(g, s) || return false
    LG.has_vertex(g, d) || return false
    flist = fadj(g)
    @inbounds list = flist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)

    g.ne += 1
    s == d && return true  # selfloop

    blist = badj(g)
    @inbounds list = blist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    return true  # edge successfully added
end

function LG.rem_edge!(g::VStaticGraph{T,N}, e) where {T,N}
    s = LG.src(e)
    d = LG.dst(e)
    LG.has_vertex(g, s) || return false
    LG.has_vertex(g, d) || return false
    flist = fadj(g)
    @inbounds list = flist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false   # edge not in graph

    deleteat!(list, index)

    g.ne -= 1
    s == d && return true  # selfloop

    blist = badj(g)
    @inbounds list = blist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    return true  # edge successfully deleted
end


text:
using Base.Test

include("../roc.jl")
@test Roc.tpNumNaive([1,11]      , [10,20,30,40], 3) == 1
@test Roc.tpNumNaive([1,11,13]   , [10,20,30,40], 3) == 2
@test Roc.tpNumNaive([1,11,14]   , [10,20,30,40], 3) == 1
@test Roc.tpNumNaive([1,11,14,30], [10,20,30,40], 3) == 2

@test Roc.tpRate([1,31   ], [10,20,30,40], 3) == 1/4
@test Roc.tpRate([1,31,44], [10,20,30,40], 3) == 1/4
@test Roc.tpRate([1,31,43], [10,20,30,40], 3) == 2/4

@test Roc.faRate([1,31   ], [10,20,30,40], 3, 100) == 1/(1+(100-4*3))
@test Roc.tpNum([1,11]      , [10,20,30,40], 3) == 1

#           T             T    T    T    T   F    F     T   F     T
detections=[0,           388, 395, 407, 418, 436, 448, 473, 487, 502]
changes   =[0, 100, 199, 300, 401, 501]

max_delay  = 30
println("... test_roc.jl > TP num. ", Roc.tpNum(detections, changes, max_delay))
println("... test_roc.jl > TP rate ", Roc.tpRate(detections, changes, max_delay))

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
linevalue(key::String, value::String) = "\"$key\":\"$value\""

function Base.contains(line::String, datafilter::RedditDataFilter)
    if isempty(datafilter.author) && isempty(datafilter.subreddit)
        return true
    else
        authorcheck = any(contains.(line, linevalue.("author", datafilter.author)))
        subredditcheck = any(contains.(line, linevalue.("subreddit", datafilter.subreddit)))
        return any([authorcheck, subredditcheck])
    end
end

function Base.contains(json::JSON3.Object, datafilter::RedditDataFilter)
    if isempty(datafilter.author) && isempty(datafilter.subreddit)
        return true
    else
        authorcheck = get(json, :author, nothing) ∈ datafilter.author
        subredditcheck = get(json, :subreddit, nothing) ∈ datafilter.subreddit
        return any([authorcheck, subredditcheck])
    end
end

function Base.filter(datafilter::RedditDataFilter, data::Dict{Symbol,Any})
    if isempty(datafilter.fields)
        return data
    else
        incols(x) = x ∈ Symbol.(datafilter.fields)
        return filter(d -> incols(d.first), data)
    end
end


text:
cnfl_residential = Residential(20, 550, -0.005, nothing,
                        [0.4, 0.3, 0.4, 0.48, 0.55, 0.63, 0.55, 0.65, 0.68, 0.7, 0.75, 0.81, 0.78, 0.72, 0.68, 0.64, 0.62, 0.65, 0.88, 0.95, 0.88, 0.75, 0.6, 0.45],
                        R_CNFL,
                        10,
                        0.0);

cnfl_commercial = CommIndus(20, 5000, 0.0, nothing,
                        [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.875, 0.8125, 0.6875, 1.0, 0.9375, 0.875, 0.8125, 0.8125, 0.8125, 0.75, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25],
                        CI_CNFL,
                        10,
                        0.10);

cnfl_industrial = CommIndus(20, 8050, 0.0, nothing,
                        [0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.666666667, 0.666666667, 0.666666667, 0.666666667, 1.0, 1.0, 1.0, 0.666666667, 0.66666666, 0.666666667, 0.666666667, 0.666666667, 1.0, 1.0, 0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037],
                        CI_CNFL,
                        10,
                        0.10);

cnfl_mediumvoltage = TMT(20, 20550, -0.005, nothing,
                        [0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.666666667, 0.666666667, 0.666666667, 0.766666667, 1.0, 1.0, 1.0, 0.766666667, 0.66666666, 0.666666667, 0.666666667, 0.766666667, 1.0, 1.0, 1.0, 0.77037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037],
                        TMT_CNFL,
                        10,
                        0.10);

ice_residential = Residential(20, 550, -0.005, nothing,
                        [0.4, 0.3, 0.4, 0.48, 0.55, 0.63, 0.55, 0.65, 0.68, 0.7, 0.75, 0.81, 0.78, 0.72, 0.68, 0.64, 0.62, 0.65, 0.88, 0.95, 0.88, 0.75, 0.6, 0.45],
                        R_ICE,
                        10,
                        0.0);

ice_commercial = CommIndus(20, 5000, 0.0, nothing,
                        [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.875, 0.8125, 0.6875, 1.0, 0.9375, 0.875, 0.8125, 0.8125, 0.8125, 0.75, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25],
                        CI_ICE,
                        10,
                        0.10);

ice_industrial = CommIndus(20, 8050, 0.0, nothing,
                        [0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.666666667, 0.666666667, 0.666666667, 0.666666667, 1.0, 1.0, 1.0, 0.666666667, 0.66666666, 0.666666667, 0.666666667, 0.666666667, 1.0, 1.0, 0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037],
                        CI_ICE,
                        10,
                        0.10);

ice_mediumvoltage = TMT(20, 20550, -0.005, nothing,
                        [0.37037037, 0.37037037, 0.37037037, 0.37037037, 0.666666667, 0.666666667, 0.666666667, 0.766666667, 1.0, 1.0, 1.0, 0.766666667, 0.66666666, 0.666666667, 0.666666667, 0.766666667, 1.0, 1.0, 1.0, 0.77037037, 0.37037037, 0.37037037, 0.37037037, 0.37037037],
                        TMT_ICE,
                        10,
                        0.10);

get_pmax(cnfl_residential);
get_pmax(cnfl_commercial);
get_pmax(cnfl_industrial);
get_pmax(cnfl_mediumvoltage);

get_pmax(ice_residential);
get_pmax(ice_commercial);
get_pmax(ice_industrial);
get_pmax(ice_mediumvoltage);

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
using LinearAlgebra
using Random
using ForwardDiff
using DataFrames
using GLM

function marquardt(f, x₀; ε=1e-6, maxiter=1000, λ=10.0^4, α=2)
    x = x₀
    fx = f(x)
    converged = false
    i = 0
    while i < maxiter
        i += 1
        ∇f = ForwardDiff.gradient(f, x)
        if norm(∇f) ≤ ε
            converged = true
            break
        end
        ∇²f = ForwardDiff.hessian(f, x)
        x′ = x - (∇²f + λ*I) \ ∇f
        fx′ = f(x′)
        if fx′ < fx
            λ *= 0.5
            fx = fx′
            x = x′
        else
            λ *= 2.0
        end
    end
    (x=x, fx=fx, converged=converged, iters=i)
end

rosenbrock(x) = sum((1-x[i])^2 + 100(x[i+1]-x[i]^2)^2 for i in 1:length(x)-1)

Random.seed!(1234)
x = rand(20);
marquardt(rosenbrock, x)

Random.seed!(123)
N = 1000
df = DataFrame(rand(N, 4))
df.y = sum.(eachrow(df)) + randn(N);
lm(@formula(y~x1+x2+x3+x4), df)

# This is not very efficient, but is good enough in our case
function ssq(x)
    m = Matrix(df)[:, startswith.(string.(names(df)), "x")]
    norm(df.y - fill(x[1], nrow(df)) - m * x[2:end])
end
marquardt(ssq, rand(5))

max_stars_repo_path: 4. Optimization/optimziation.jl

max_issues_repo_path: 4. Optimization/optimziation.jl

max_forks_repo_path: 4. Optimization/optimziation.jl


text:
include("bipartite.jl")

"""
    all_different(variables::Vector{Variable})

Create a BasicConstraint which will later be used by `all_different(com, constraint)`. \n
Can be used i.e by `add_constraint!(com, CS.all_different(variables))`.
"""
function all_different(variables::Vector{Variable})
    constraint = BasicConstraint()
    constraint.fct = all_different
    constraint.indices = Int[v.idx for v in variables]
    constraint.hash = constraint_hash(constraint)
    return constraint
end

"""
    all_different(com::CS.CoM, constraint::BasicConstraint; logs = true)

Tries to reduce the search space by the all_different constraint.
Fixes values and then sets com.changed to true for the corresponding index.
Return a ConstraintOutput object and throws a warning if infeasible and `logs` is set
"""
function all_different(com::CS.CoM, constraint::BasicConstraint; logs = true)
    indices = constraint.indices
    pvals = constraint.pvals
    nindices = length(indices)

    search_space = com.search_space
    fixed_vals, unfixed_indices = fixed_vs_unfixed(search_space, indices)

    fixed_vals_set = Set(fixed_vals)
    # check if one value is used more than once
    if length(fixed_vals_set) < length(fixed_vals)
        logs && @warn "The problem is infeasible"
        return false
    end

    bfixed = true
    current_fixed_vals = fixed_vals

    while bfixed
        new_fixed_vals = Int[]
        bfixed = false
        for i in 1:length(unfixed_indices)
            pi = unfixed_indices[i]
            ind = indices[pi]
            @views c_search_space = search_space[ind]
            if !CS.isfixed(c_search_space)
                for pv in current_fixed_vals
                    if has(c_search_space, pv)
                        if !rm!(com, c_search_space, pv)
                            logs && @warn "The problem is infeasible"
                            return false
                        end

                        if nvalues(c_search_space) == 1
                            only_value = CS.value(c_search_space)
                            push!(fixed_vals, only_value)
                            push!(new_fixed_vals, only_value)
                            bfixed = true
                            break
                        end
                    end
                end
            end
        end
        current_fixed_vals = new_fixed_vals
    end

    if length(fixed_vals) == nindices
        return true
    end

    min_pvals, max_pvals = extrema(pvals)
    len_range = max_pvals-min_pvals+1

    # find maximum_matching for infeasible check and Berge's lemma
    # building graph
    # ei = indices, ej = possible values
    # i.e ei=[1,2,1] ej = [1,2,2] => 1->[1,2], 2->[2]
    min_pvals_m1 = min_pvals-1

    pval_mapping = zeros(Int, length(pvals))
    vertex_mapping = zeros(Int, len_range)
    vertex_mapping_bw = zeros(Int, nindices+length(pvals))

    vc = 1
    for i in indices
        vertex_mapping_bw[vc] = i
        vc += 1
    end
    pvc = 1
    for pv in pvals
        pval_mapping[pvc] = pv
        vertex_mapping[pv-min_pvals_m1] = vc
        vertex_mapping_bw[vc] = pv
        vc += 1
        pvc += 1
    end
    num_nodes = vc

    # count the number of edges
    num_edges = 0
    @inbounds for i in indices
        num_edges += nvalues(search_space[i])
    end

    di_ei = Vector{Int64}(undef,num_edges)
    di_ej = Vector{Int64}(undef,num_edges)

    # add edge from each index to the possible values
    edge_counter = 0
    vc = 0
    @inbounds for i in indices
        vc += 1
        if isfixed(search_space[i])
            edge_counter += 1
            di_ei[edge_counter] = vc
            di_ej[edge_counter] = vertex_mapping[CS.value(search_space[i])-min_pvals_m1]-nindices
        else
            for pv in view_values(search_space[i])
                edge_counter += 1
                di_ei[edge_counter] = vc
                di_ej[edge_counter] = vertex_mapping[pv-min_pvals_m1]-nindices
            end
        end
    end

    maximum_matching = bipartite_cardinality_matching(di_ei, di_ej, vc, len_range; l_sorted=true)
    if maximum_matching.weight != nindices
        logs && @warn "Infeasible (No maximum matching was found)"
        return false
    end

    # directed edges for strongly connected components
    vc = 0
    edge_counter = 0
    @inbounds for i in indices
        vc += 1
        if isfixed(search_space[i])
            edge_counter += 1
            di_ei[edge_counter] = vc
            di_ej[edge_counter] = vertex_mapping[CS.value(search_space[i])-min_pvals_m1]
        else
            for pv in view_values(search_space[i])
                edge_counter += 1
                if pv == pval_mapping[maximum_matching.match[vc]]
                    di_ei[edge_counter] = vc
                    di_ej[edge_counter] = vertex_mapping[pv-min_pvals_m1]
                else
                    di_ei[edge_counter] = vertex_mapping[pv-min_pvals_m1]
                    di_ej[edge_counter] = vc
                end
            end
        end
    end

    # if we have more values than indices
    if length(pvals) > nindices
        # add extra node which has all values as input which are used in the maximum matching
        # and the ones in maximum matching as inputs see: http://www.minicp.org (Part 6)
        # the direction is opposite to that of minicp
        used_in_maximum_matching = Dict{Int, Bool}()
        for pval in pvals
            used_in_maximum_matching[pval] = false
        end
        vc = 0
        for i in indices
            vc += 1
            for pv in view_values(search_space[i])
                if pv == pval_mapping[maximum_matching.match[vc]]
                    used_in_maximum_matching[pv] = true
                    break
                end
            end
        end
        new_vertex = num_nodes+1
        for kv in used_in_maximum_matching
            # not in maximum matching
            if !kv.second
                push!(di_ei, new_vertex)
                push!(di_ej, vertex_mapping[kv.first-min_pvals_m1])
            else
                push!(di_ei, vertex_mapping[kv.first-min_pvals_m1])
                push!(di_ej, new_vertex)
            end
        end
    end

    sccs_map = strong_components_map(di_ei, di_ej)

    # remove the left over edges from the search space
    vmb = vertex_mapping_bw
    for (src,dst) in zip(di_ei, di_ej)
        # edges to the extra node don't count
        if src > num_nodes || dst > num_nodes
            continue
        end

        # if in same strong component -> part of a cycle -> part of a maximum matching
        if sccs_map[src] == sccs_map[dst]
            continue
        end
        #  remove edges in maximum matching and then edges which are part of a cycle
        if src <= nindices && dst == vertex_mapping[pval_mapping[maximum_matching.match[src]]-min_pvals_m1]
            continue
        end

        cind = vmb[dst]
        if !rm!(com, search_space[cind], vmb[src])
            logs && @warn "The problem is infeasible"
            return false
        end

        # if only one value possible make it fixed
        if nvalues(search_space[cind]) == 1
            only_value = CS.value(search_space[cind])
        end
    end

    return true
end

"""
    all_different(com::CoM, constraint::Constraint, value::Int, index::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `index` to `value`.
"""
function all_different(com::CoM, constraint::Constraint, value::Int, index::Int)
    indices = constraint.indices
    for i=1:length(indices)
        if indices[i] == index
            continue
        end
        if issetto(com.search_space[indices[i]], value)
            return false
        end
    end
    return true
end


text:
"""
    QuantisedDiatomic

This module exports two user facing functions:

- `generate_configurations`
    Creates a set of velocities and positions for diatomic molecule with specified
    vibrational `ν` and rotational `J` quantum numbers.

- `quantise_diatomic`
    Obtains vibrational `ν` and rotational `J` quantum numbers for a diatomic molecule
    with a given set of velocities and positions.

The central concept of this module is the EBK procedure which is nicely detailed here:
[Larkoski2006](@cite)

Inspired by VENUS96: [Hase1996](@cite)
"""
module QuantisedDiatomic

using QuadGK: QuadGK
using Rotations: Rotations
using LsqFit: LsqFit
using LinearAlgebra: norm, normalize!, nullspace, cross, I
using UnitfulAtomic: austrip
using Distributions: Uniform
using UnicodePlots: lineplot, lineplot!, DotCanvas
using Optim: Optim

using NQCDynamics: Simulation, Calculators, DynamicsUtils
using NQCBase: Atoms, PeriodicCell, InfiniteCell
using NQCModels: NQCModels
using NQCModels.AdiabaticModels: AdiabaticModel

export generate_configurations
export quantise_diatomic

struct SurfaceParameters{T}
    reduced_mass::T
    atom_indices::Vector{Int}
    slab::Matrix{T}
    height::T
    surface_normal::Vector{T}
    orthogonal_vectors::Matrix{T}
end

function SurfaceParameters(masses::AbstractVector, atom_indices::AbstractVector, slab::Matrix, height, surface_normal)
    μ = reduced_mass([masses[i] for i in atom_indices])
    orthogonal_vectors = nullspace(reshape(surface_normal, 1, :)) # Plane parallel to surface
    if size(slab, 2) == 0
        surface_top = 0
    else
        surface_top = maximum(slab[3,:])
    end
    return SurfaceParameters{eltype(masses)}(μ, atom_indices, slab, height + surface_top,
        surface_normal, orthogonal_vectors)
end

struct GenerationParameters{T}
    direction::Vector{T}
    translational_energy::T
    samples::Int
end

Base.broadcastable(p::SurfaceParameters) = Ref(p)
Base.broadcastable(p::GenerationParameters) = Ref(p)

"""
    generate_configurations(sim, ν, J; samples=1000, height=10, normal_vector=[0, 0, 1],
        translational_energy=0, direction=[0, 0, -1], position=[0, 0, height])

Generate positions and momenta for given quantum numbers

`translational_energy`, `direction` and `position` specify the kinetic energy in
a specific direction with the molecule placed with centre of mass at `position`.

Keyword arguments `height` and `normal_vector` become relevant if the potential
requires specific placement of the molecule.
These allow the molecule to be placed at a distance `height` in the direction
`normal_vector` when performing potential evaluations.
"""
function generate_configurations(sim, ν, J;
    samples=1000,
    height=10.0,
    surface_normal=[0, 0, 1.0],
    translational_energy=0.0,
    direction=[0, 0, -1.0],
    atom_indices=[1, 2],
    r=zeros(size(sim)))

    _, slab = separate_slab_and_molecule(atom_indices, r)
    surface = SurfaceParameters(sim.atoms.masses, atom_indices, slab, height, surface_normal)
    generation = GenerationParameters(direction, austrip(translational_energy), samples)

    total_energy, bounds = find_total_energy(sim.calculator.model, ν, J, surface)

    radial_momentum = get_radial_momentum_function(total_energy, surface, J, sim.calculator.model)
    bonds, momenta = select_random_bond_lengths(bounds, radial_momentum, samples)

    @info "Generating the requested configurations..."
    configure_diatomic.(sim, bonds, momenta, J, surface, generation)
end

"""
    separate_slab_and_molecule(atom_indices, r)

Get the coordinates of the molecule and slab separately.
"""
function separate_slab_and_molecule(atom_indices, r)
    molecule = r[:,atom_indices]
    slab_indices = [i for i in axes(r, 2) if i ∉ atom_indices]
    slab = r[:,slab_indices]
    return (molecule, slab)
end

"""
    combine_slab_and_molecule(atom_indices, molecule, slab)

Revert the transformation `separate_slab_and_molecule`
"""
function combine_slab_and_molecule(atom_indices, molecule, slab)
    r = zeros(3, size(molecule, 2) + size(slab, 2))
    slab_indices = [i for i in axes(r, 2) if i ∉ atom_indices]
    r[:,atom_indices] .= molecule
    r[:,slab_indices] .= slab
    return r
end

function combine_slab_and_molecule(surface::SurfaceParameters, molecule)
    combine_slab_and_molecule(surface.atom_indices, molecule, surface.slab)
end

"""
    find_total_energy(sim, ν, J, surface)

Returns the energy associated with the specified quantum numbers
"""
function find_total_energy(model::AdiabaticModel, ν, J, surface::SurfaceParameters)

    μ = surface.reduced_mass

    bond_lengths = 0.5:0.01:5.0
    potential = calculate_diatomic_energy.(model, bond_lengths, surface) # Calculate binding curve
    potential_minimum, index = findmin(potential)
    r₀ = bond_lengths[index]
    potential .-= potential_minimum # Shift minimum to zero.

    k = calculate_force_constant(bond_lengths, potential)
    E = guess_initial_energy(k, μ, ν, potential_minimum)

    V(r) = effective_potential(r, J, model, surface)

    "Evaluate vibrational quantum number."
    function nᵣ(E, r₁, r₂)
        kernel(r) = sqrt(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/π * integral - 1/2
    end

    "Derivative of the above"
    function ∂nᵣ∂E(E, r₁, r₂)
        kernel(r) = 1 / sqrt(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/2π * integral
    end

    "Second derivative of the above"
    function ∂²nᵣ∂E²(E, r₁, r₂)
        kernel(r) = 1 / sqrt(E - V(r))^3
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return -sqrt(2μ)/4π * integral
    end

    """
    Evaluates the function, gradient and hessian as described in the Optim documentation.
    """
    function fgh!(F,G,H,x)
        E = x[1]
        r₁, r₂ = find_integral_bounds(model, E, J, surface, r₀)

        νtmp = nᵣ(E, r₁, r₂)
        difference = νtmp - ν

        if G !== nothing
            G[1] = ∂nᵣ∂E(E, r₁, r₂) * sign(difference)
        end
        if H !== nothing
            H[1] = ∂²nᵣ∂E²(E, r₁, r₂) * sign(difference)
        end
        if F !== nothing
            return abs(difference)
        end
    end

    @info "Converging the total energy to match the chosen quantum numbers..."
    optim = Optim.optimize(Optim.only_fgh!(fgh!), [E], Optim.Newton(),
        Optim.Options(show_trace=true, x_tol=1e-3)
    )
    E = Optim.minimizer(optim)[1]
    bounds = find_integral_bounds(model, E, J, surface, r₀)

    E, bounds
end

"""
    quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix;
        height=10, normal_vector=[0, 0, 1])

Quantise the vibrational and rotational degrees of freedom for the specified
positions and velocities

When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`.
If the potential is independent of centre of mass position, this has no effect.
Otherwise, be sure to modify these parameters to give the intended behaviour.
"""
function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix;
    height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2])

    r, slab = separate_slab_and_molecule(atom_indices, r)
    surface = SurfaceParameters(sim.atoms.masses, atom_indices, slab, austrip(height), surface_normal)

    r_com = subtract_centre_of_mass(r, [sim.atoms.masses[i] for i in atom_indices])
    v_com = subtract_centre_of_mass(v, [sim.atoms.masses[i] for i in atom_indices])
    p_com = v_com .* [sim.atoms.masses[i] for i in atom_indices]'

    k = DynamicsUtils.classical_kinetic_energy(sim, v_com)
    p = calculate_diatomic_energy(sim.calculator.model, bond_length(r_com), surface)
    E = k + p

    L = total_angular_momentum(r_com, p_com)
    J = (sqrt(1+4*L^2) - 1) / 2 # L^2 = J(J+1)ħ^2

    r₁, r₂ = find_integral_bounds(sim.calculator.model, E, J, surface, bond_length(r_com))

    μ = surface.reduced_mass
    V(r) = effective_potential(r, J, sim.calculator.model, surface)

    function nᵣ(E, r₁, r₂)
        kernel(r) = sqrt(E - V(r))
        integral, _ = QuadGK.quadgk(kernel, r₁, r₂; maxevals=100)
        return sqrt(2μ)/π * integral - 1/2
    end

    ν = nᵣ(E, r₁, r₂)

    ν, J
end

rotational_energy(r, μ, J) = J*(J+1) / (2μ*r^2)

function effective_potential(r, J, model, surface)
    rotational = rotational_energy(r, surface.reduced_mass, J)
    potential = calculate_diatomic_energy(model, r, surface)
    return rotational + potential
end

function get_radial_momentum_function(total_energy, surface, J, model)
    radial_momentum(r) = sqrt(2surface.reduced_mass * (total_energy - effective_potential(r, J, model, surface)))
end

function find_integral_bounds(model, total_energy, J, surface, r₀, bond_limits=(0.5, 5.0))
    optim_func(r) = abs(total_energy - effective_potential(r, J, model, surface))
    optim = Optim.optimize(optim_func, first(bond_limits), r₀)
    r₁ = Optim.minimizer(optim)
    optim = Optim.optimize(optim_func, r₀, last(bond_limits))
    r₂ = Optim.minimizer(optim)
    return r₁, r₂
end

"""
Pick a random bond length and corresponding radial momentum that matches the
radial probability distribution.

Uses rejection sampling: https://en.wikipedia.org/wiki/Rejection_sampling
"""
function select_random_bond_lengths(bounds, momentum_function, samples)
    @info "Generating distribution of bond lengths and radial momenta..."
    M = 1 / momentum_function(bounds[1]+1e-3)
    radius_proposal = Uniform(bounds...)
    bonds = zeros(samples)
    momenta = zeros(samples)
    for i=1:samples
        keep_going = true
        while keep_going
            r = rand(radius_proposal) # Generate radius
            p = momentum_function(r) # Evaluate probability
            probability = 1 / p # Probability inversely proportional to momentum
            if rand() < probability / M
                bonds[i] = r
                momenta[i] = rand([-1, 1]) * p
                keep_going = false
            end
        end
    end
    bonds, momenta
end

"""
Randomly orient molecule in space for a given bond length and radial momentum
"""
function configure_diatomic(sim, bond, momentum, J, surface::SurfaceParameters, generation::GenerationParameters)
    r = zeros(3,2)
    v = zeros(3,2)
    r[1,1] = bond
    v[1,1] = momentum ./ surface.reduced_mass
    r = subtract_centre_of_mass(r, [sim.atoms.masses[i] for i in surface.atom_indices])
    v = subtract_centre_of_mass(v, [sim.atoms.masses[i] for i in surface.atom_indices])

    L = sqrt(J * (J + 1))

    random = 2π*rand()
    ω = L .* [0, sin(random), cos(random)] ./ (surface.reduced_mass * bond^2)

    v[:,1] .+= cross(ω, r[:,1])
    v[:,2] .+= cross(ω, r[:,2])

    apply_random_rotation!(v, r)
    position_above_surface!(r, surface.height, sim.cell)
    apply_translational_impulse!(v, [sim.atoms.masses[i] for i in surface.atom_indices], generation.translational_energy, generation.direction)

    v, r
end

"""
Randomly rotate each column of two 3*N matrix, same rotation for all columns.
"""
function apply_random_rotation!(x, y)
    rotation = rand(Rotations.RotMatrix{3})
    for (col1, col2) in zip(eachcol(x), eachcol(y))
        col1 .= rotation * col1
        col2 .= rotation * col2
    end
end

function position_above_surface!(r, height, cell::PeriodicCell)
    r[3,:] .+= height
    a1 = cell.vectors[:,1]
    a2 = cell.vectors[:,2]
    displacement = rand()*a1+rand()*a2
    r .+= displacement
    return nothing
end

function position_above_surface!(r, height, ::InfiniteCell)
    r[3,:] .+= height
    return nothing
end

function apply_translational_impulse!(v, masses, translational_energy, direction)
    velocity_impulse = velocity_from_energy(masses, translational_energy)
    v .+= velocity_impulse .* normalize!(direction)
end

function velocity_from_energy(masses, energy)
    m = sum(masses)
    sqrt(2austrip(energy) / m)
end

"""
    calculate_diatomic_energy(model::AdiabaticModel, bond_length::Real;
        height=10, normal_vector=[0, 0, 1])

Returns potential energy of diatomic with `bond_length` at `height` from surface.

Orients molecule parallel to the surface at the specified height, the surface is assumed
to intersect the origin.
This requires that the model implicitly provides the surface, or works fine without one.
"""
function calculate_diatomic_energy(
    model::AdiabaticModel, r₀::Real, surface::SurfaceParameters
)

    (;surface_normal, orthogonal_vectors, height) = surface
    molecule = surface_normal .* height .+ orthogonal_vectors .* r₀ ./ sqrt(2)
    R = combine_slab_and_molecule(surface, molecule)

    return NQCModels.potential(model, R)
end

"""
    calculate_force_constant(sim::Simulation, surface::SurfaceParameters)

Evaluate energy for different bond lengths and identify force constant.

This uses LsqFit.jl to fit a morse potential with optimised minimum.
"""
function calculate_force_constant(bond_lengths, potential)

    morse!(f, x, p) = (@. f = p[1]*(1 - exp(-p[2]*(x-p[3])))^2)
    function jacobian!(J, x, p)
        @. J[:,1] = (1 - exp(-p[2]*(x-p[3])))^2
        @. J[:,2] = 2p[1]*(1 - exp(-p[2]*(x-p[3]))) * exp(-p[2]*(x-p[3])) * (x - p[3])
        @. J[:,3] = -2p[1]*(1 - exp(-p[2]*(x-p[3]))) * exp(-p[2]*(x-p[3])) * p[2]
    end

    σ = potential[end]
    weights = exp.(-potential.^2 ./ σ^2) # Weight the area around the minimum most heavily

    fit = LsqFit.curve_fit(morse!, jacobian!, bond_lengths, potential, weights, [1.0, 1.0, 1.0];
        lower=[0.0, 0.0, 0.0], inplace=true)

    plot_binding_curve(morse!, bond_lengths, potential, fit.param)

    force_constant = fit.param[2]^2 * 2 * fit.param[1]
    return force_constant
end

function plot_binding_curve(morse!, bond_lengths, V, parameters)
    fitted_points = zero(bond_lengths)
    morse!(fitted_points, bond_lengths, parameters)
    plt = lineplot(bond_lengths, V;
        title="Binding curve", xlabel="Bond length / bohr", ylabel="Energy / Hartree",
        name="Actual values", canvas=DotCanvas, border=:ascii
    )
    lineplot!(plt, bond_lengths, fitted_points; name="Fitted curve")
    show(plt)
    println()
    @info "The two lines shown above should closely match. \
        This indicates the evaluation of the potential is working correctly."
end

function guess_initial_energy(k, μ, ν, potential_minimum)
    ω = sqrt(k / μ)
    E_guess = (ν + 1/2) * ω
    return E_guess + potential_minimum
end

subtract_centre_of_mass(x, m) = x .- centre_of_mass(x, m)
@views centre_of_mass(x, m) = (x[:,1].*m[1] .+ x[:,2].*m[2]) ./ (m[1]+m[2])
reduced_mass(atoms::Atoms) = reduced_mass(atoms.masses)
reduced_mass(m::AbstractVector) = m[1]*m[2]/(m[1]+m[2])
@views bond_length(r) = norm(r[:,1] .- r[:,2])
@views total_angular_momentum(r, p) = norm(cross(r[:,1], p[:,1]) + cross(r[:,2], p[:,2]))

end # module


text:
############################################################################################
"""
$(SIGNATURES)
Check if 'θ' is of finite value and return Bool. Not exported.

# Examples
```julia
```

"""
function _checkfinite(θ::T, max_val::R=max_val) where {T<:Real,R<:Real}
    return isfinite(θ) && θ < max_val
end
function _checkfinite(θ::Array{T}, max_val::R=max_val) where {T<:Real,R<:Real}
    @inbounds @simd for iter in eachindex(θ)
        if !_checkfinite(θ[iter], max_val)
            return false
        end
    end
    return true
end
function _checkfinite(θ::AbstractArray, max_val::R=max_val) where {R<:Real}
    @inbounds @simd for iter in eachindex(θ)
        if !_checkfinite(θ[iter], max_val)
            return false
        end
    end
    return true
end
function _checkfinite(θ::Tuple, max_val::R=max_val) where {R<:Real}
    return all(map(_checkfinite, θ))
end
function _checkfinite(θ::N, max_val::R=max_val) where {N<:NamedTuple} where {R<:Real}
    return all(map(_checkfinite, θ))
end

############################################################################################
"""
$(SIGNATURES)
Check if 'prior' is a valid density and return Bool. Not exported.

# Examples
```julia
```

"""
function _checkprior(prior)
    return false
end
function _checkprior(prior::S) where {S<:Distributions.Distribution}
    return true
end
function _checkprior(prior::AbstractVector{T}) where {T}
    @inbounds @simd for iter in eachindex(prior)
        if !_checkprior(prior[iter])
            return false
        end
    end
    return true
end
function _checkprior(prior::AbstractArray)
    @inbounds @simd for iter in eachindex(prior)
        if !_checkprior(prior[iter])
            return false
        end
    end
    return true
end
function _checkprior(prior::N) where {N<:NamedTuple}
    return all(map(_checkprior, prior))
end

############################################################################################
"""
$(SIGNATURES)
Check if argument is not fixed. Returns NamedTuple with true/false. Needed in addition to _checkprior for nested NamedTuples. Not exported.

# Examples
```julia
```

"""
function _checksampleable(constraint)
    return false
end
function _checksampleable(constraint::S) where {S<:Distributions.Distribution}
    return true
end
function _checksampleable(constraint::AbstractVector{T}) where {T}
    @inbounds @simd for iter in eachindex(constraint)
        if !_checksampleable(constraint[iter])
            return false
        end
    end
    return true
end
function _checksampleable(constraint::AbstractArray)
    @inbounds @simd for iter in eachindex(constraint)
        if !_checksampleable(constraint[iter])
            return false
        end
    end
    return true
end
function _checksampleable(constraint::NamedTuple{names}) where {names}
    return NamedTuple{names}(map(_checksampleable, constraint))
end

############################################################################################
#=
"""
$(SIGNATURES)
Check if provided val-constraint combination is valid for Param struct. Not exported.

# Examples
```julia
```

"""
function _checkparam(
    _rng::Random.AbstractRNG,
    constraint::AbstractArray,
    val::AbstractArray,
)
    return all(
        map(constraint, val) do constraintᵢ, valᵢ
            _checkparam(_rng, constraintᵢ, valᵢ)
        end,
    )
end
_checkparam(constraint, val) = _checkparam(Random.GLOBAL_RNG, constraint, val)
=#
############################################################################################
"""
$(SIGNATURES)
Check if all values in (Nested) NamedTuple are a 'Param' struct and return Bool. Not exported.

# Examples
```julia
```

"""
function _checkparams(param)
    return false
end
function _checkparams(param::Param)
    return true
end

function _checkparams(param::N) where {N<:NamedTuple}
    return all(map(_checkparams, param))
end

############################################################################################
# Export


text:
"""
    CR_logreg(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame)

Tests and ranks (F1 score) all possible set of variables fitting a logistic regression.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the comprehensive research.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
"""
function CR_logreg(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame)

    diagnostic = DataFrame(Model = Array{Symbol}[], Deviance = Float64[], F₁ = Float64[], k = Int64[])
    threshold = collect(0.1:.01:.6)
    model_var = collect(combinations(varnames))

    for i in 1:length(model_var)
        vars = model_var[i]
        model = glm(Term(Symbol(topredict)) ~ sum(Term.(vars)), train, Bernoulli(), LogitLink())
        Ŷ = convert(Vector{Float64}, GLM.predict(model, test))
        r = MLBase.roc(test[:, Symbol(topredict)], Ŷ, threshold)
        push!(diagnostic, [model_var[i], -deviance(model), maximum(f1score.(r)), length(vars)])
    end

    sort!(diagnostic, [:F₁, :Deviance], rev=true)
end

"""
    stepF1_logreg(topredict::String, varnames::Vector{Symbol}, initial_list::Union{Vector{Symbol}, Vector{Any}}, train::DataFrame, test::DataFrame, verbose::Bool=true)

Performs a forward feature selection with logistic regression selecting on the F1 score.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the regression.
- `initial_list::Union{Vector{Symbol}, Vector{Any}}`: List of features to start with.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
- `verbose::Bool=true`: whether to print the sequence of inclusions and exclusions. Default = true.
"""

function stepF1_logreg(topredict::String, varnames::Vector{Symbol}, initial_list::Union{Vector{Symbol}, Vector{Any}}, train::DataFrame, test::DataFrame, verbose::Bool=true)
    included = initial_list
    while true
        changed=false

        # forward step
        excluded = setdiff(varnames, included)
        new_F1 = DataFrame(var = [excluded[i] for i in 1:length(excluded)])
        vec_F1 = [] #Will store the F1 scores to then add them to the new_F1 df
        for i=1:length(excluded)
            concat = vcat(included,excluded[i]) #concatenates included and one variable from excluded
            model = glm(Term(Symbol(topredict)) ~ sum(Term.(concat)), train, Bernoulli(), LogitLink())
            predictions = convert(Vector{Float64}, GLM.predict(model, test))
            threshold = collect(0.1:.01:.6)
            r = MLBase.roc(test[:, Symbol(topredict)], predictions, threshold)
            F1 = maximum(f1score.(r))
            push!(vec_F1, F1)
        end
        new_F1.F1score = vec_F1

        best_F1 = maximum(new_F1.F1score)

        if included != []
            model = glm(Term(Symbol(topredict)) ~ sum(Term.(included)), train, Bernoulli(), LogitLink())
            predictions = convert(Vector{Float64},GLM.predict(model, test))
            threshold = collect(0.1:.01:.6)
            r = MLBase.roc(test[:, Symbol(topredict)], predictions, threshold)
            ex_F1 = maximum(f1score.(r))
        else
            ex_F1 = 0
        end
        if  best_F1 > ex_F1
            best_feature = new_F1[new_F1.F1score .== best_F1,:var][1] #[1] is to extract a Symbol and not a Vector
            push!(included, best_feature)
            changed=true
            if verbose
                println("Add $(best_feature) with F1 $(best_F1)")
            end
        end

        if changed == false
            println("F1 score: $(ex_F1)")
            return convert(Vector{Symbol}, included)
            break
        end
    end
    return convert(Vector{Symbol}, included)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
module ReinforcementLearningAnIntroduction

using Reexport

const RLIntro = ReinforcementLearningAnIntroduction
export RLIntro

@reexport using ReinforcementLearningBase
@reexport using ReinforcementLearningCore

# !!! just a quick and dirty fix for
# https://github.com/JuliaReinforcementLearning/ReinforcementLearningAnIntroduction.jl/issues/17
function RLBase.update!(model::AbstractEnvironmentModel, buffer::AbstractTrajectory)
    transitions = extract_experience(buffer, model)
    isnothing(transitions) || update!(model, transitions)
end

include("environments/environments.jl")
include("extensions/extensions.jl")

end # module