text:
module AWS

using Compat: Compat
using Base64
using Dates
using HTTP
using MbedTLS
using Mocking
using OrderedCollections: LittleDict, OrderedDict
using Retry
using Sockets
using UUIDs: UUIDs
using XMLDict

export @service
export _merge, AbstractAWSConfig, AWSConfig, AWSExceptions, AWSServices, Request
export global_aws_config, generate_service_url, set_user_agent, sign!, sign_aws2!, sign_aws4!
export JSONService, RestJSONService, RestXMLService, QueryService

include("utilities.jl")
include("AWSExceptions.jl")
include("AWSCredentials.jl")
include("AWSConfig.jl")
include("AWSMetadataUtilities.jl")
include("AWSMetadata.jl")

using ..AWSExceptions
using ..AWSExceptions: AWSException

const user_agent = Ref("AWS.jl/1.0.0")
const aws_config = Ref{AbstractAWSConfig}()

"""
    global_aws_config()

Retrieve the global AWS configuration.
If one is not set, create one with default configuration options.

# Keywords
- `kwargs...`: AWSConfig kwargs to be passed along if the global configuration is not already set

# Returns
- `AWSConfig`: The global AWS configuration
"""
function global_aws_config(; kwargs...)
    if !isassigned(aws_config) || !isempty(kwargs)
        aws_config[] = AWSConfig(; kwargs...)
    end

    return aws_config[]
end

"""
    global_aws_config(config::AbstractAWSConfig)

Set the global AWSConfig.

# Arguments
- `config::AWSConfig`: The AWSConfig to set in the global state

# Returns
- `AWSConfig`: Global AWSConfig
"""
function global_aws_config(config::AbstractAWSConfig)
    return aws_config[] = config
end

"""
    set_user_agent(new_user_agent::String)

Set the global user agent when making HTTP requests.

# Arguments
- `new_user_agent::String`: User agent to set when making HTTP requests

# Return
- `String`: The global user agent
"""
set_user_agent(new_user_agent::String) = return user_agent[] = new_user_agent

"""
    macro service(module_name::Symbol)

Include a high-level service wrapper based off of the module_name parameter.

When calling the macro you cannot match the predefined constant for the lowl level API.
The low level API constants are named in all lowercase, and spaces replaced with underscores.

Examples:
```julia
using AWS.AWSServices: secrets_manager
using AWS: @service

# This matches the constant and will error!
@service secrets_manager
> ERROR: cannot assign a value to variable AWSServices.secrets_manager from module Main

# This does NOT match the filename structure and will error!
@service secretsmanager
> ERROR: could not open file /.julia/dev/AWS.jl/src/services/secretsmanager.jl

# All of the examples below are valid!
@service Secrets_Manager
@service SECRETS_MANAGER
@service sECRETS_MANAGER
```

# Arguments
- `module_name::Symbol`: Name of the service to include high-level API wrappers in your namespace

# Return
- `Expression`: Base.include() call to introduce the high-level service API wrapper functions in your namespace
"""
macro service(module_name::Symbol)
    service_name = joinpath(@__DIR__, "services", lowercase(string(module_name)) * ".jl")

    return Expr(:toplevel,
    :(module($(esc(module_name)))
        Base.include($(esc(module_name)), $(esc(service_name)))
     end))
end

struct RestXMLService
    name::String
    api_version::String
end

struct QueryService
    name::String
    api_version::String
end

struct JSONService
    name::String
    api_version::String

    json_version::String
    target::String
end

struct RestJSONService
    name::String
    api_version::String

    service_specific_headers::LittleDict{String, String}
end

RestJSONService(name::String, api_version::String) = RestJSONService(name, api_version, LittleDict{String, String}())

Base.@kwdef mutable struct Request
    service::String
    api_version::String
    request_method::String

    headers::AbstractDict{String, String}=LittleDict{String, String}()
    content::Union{String, Vector{UInt8}}=""
    resource::String=""
    url::String=""

    return_stream::Bool=false
    response_stream::Union{Base.BufferStream, Nothing}=nothing
    http_options::Array{String}=[]
    return_raw::Bool=false
    response_dict_type::Type{<:AbstractDict}=LittleDict
end

# Needs to be included after the definition of struct otherwise it cannot find them
include("AWSServices.jl")

_http_status(e::HTTP.StatusError) = e.status
_header(e::HTTP.StatusError, k, d="") = HTTP.header(e.response, k, d)

function sign!(aws::AbstractAWSConfig, request::Request; time::DateTime=now(Dates.UTC))
    if request.service in ("sdb", "importexport")
        sign_aws2!(aws, request, time)
    else
        sign_aws4!(aws, request, time)
    end
end

function sign_aws2!(aws::AbstractAWSConfig, request::Request, time::DateTime)
    # Create AWS Signature Version 2 Authentication query parameters.
    # http://docs.aws.amazon.com/general/latest/gr/signature-version-2.html

    query = Dict{String, String}()
    for elem in split(request.content, '&', keepempty=false)
        (n, v) = split(elem, "=")
        query[n] = HTTP.unescapeuri(v)
    end

    request.headers["Content-Type"] = "application/x-www-form-urlencoded; charset=utf-8"

    creds = check_credentials(credentials(aws))
    query["AWSAccessKeyId"] = creds.access_key_id
    query["Expires"] = Dates.format(time + Dates.Minute(2), dateformat"yyyy-mm-dd\THH:MM:SS\Z")
    query["SignatureVersion"] = "2"
    query["SignatureMethod"] = "HmacSHA256"

    if !isempty(creds.token)
        query["SecurityToken"] = creds.token
    end

    query = [k => query[k] for k in sort!(collect(keys(query)))]
    uri = HTTP.URI(request.url)
    to_sign = "POST\n$(uri.host)\n$(uri.path)\n$(HTTP.escapeuri(query))"
    push!(query, "Signature" => digest(MD_SHA256, to_sign, creds.secret_key) |> base64encode |> strip)

    request.content = HTTP.escapeuri(query)

    return request
end

function sign_aws4!(aws::AbstractAWSConfig, request::Request, time::DateTime)
    # Create AWS Signature Version 4 Authentication Headers.
    # http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html

    date = Dates.format(time, dateformat"yyyymmdd")
    datetime = Dates.format(time, dateformat"yyyymmdd\THHMMSS\Z")

    # Authentication scope...
    authentication_scope = [date, region(aws), request.service, "aws4_request"]

    creds = check_credentials(credentials(aws))
    signing_key = "AWS4$(creds.secret_key)"

    for scope in authentication_scope
        signing_key = digest(MD_SHA256, scope, signing_key)
    end

    # Authentication scope string...
    authentication_scope = join(authentication_scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(digest(MD_SHA256, request.content))

    # HTTP headers...
    delete!(request.headers, "Authorization")

    merge!(request.headers, Dict(
        "x-amz-content-sha256" => content_hash,
        "x-amz-date" => datetime,
        "Content-MD5" => base64encode(digest(MD_MD5, request.content))
    ))

    if !isempty(creds.token)
        request.headers["x-amz-security-token"] = creds.token
    end

    # Sort and lowercase() Headers to produce canonical form...
    canonical_headers = join(sort!(["$(lowercase(k)):$(strip(v))" for (k,v) in request.headers]), "\n")
    signed_headers = join(sort!([lowercase(k) for k in keys(request.headers)]), ";")

    # Sort Query String...
    uri = HTTP.URI(request.url)
    query = HTTP.URIs.queryparams(uri.query)
    query = [k => query[k] for k in sort!(collect(keys(query)))]

    # Create hash of canonical request...
    canonical_form = string(
        request.request_method, "\n",
        request.service == "s3" ? uri.path : HTTP.escapepath(uri.path), "\n",
        HTTP.escapeuri(query), "\n",
        canonical_headers, "\n\n",
        signed_headers, "\n",
        content_hash
    )

    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$authentication_scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    # Append Authorization header...
    request.headers["Authorization"] = join([
            "AWS4-HMAC-SHA256 Credential=$(creds.access_key_id)/$authentication_scope",
            "SignedHeaders=$signed_headers",
            "Signature=$signature",
        ],
        ", ",
    )

    return request
end

function generate_service_url(aws::AbstractAWSConfig, service::String, resource::String)
    SERVICE_HOST = "amazonaws.com"
    reg = region(aws)

    regionless_services = ("iam", "route53")

    if service in regionless_services || (service == "sdb" && reg == "us-east-1")
        reg = ""
    end

    return string("https://", service, ".", isempty(reg) ? "" : "$reg.", SERVICE_HOST, resource)
end

function _http_request(request::Request)
    @repeat 4 try
        http_stack = HTTP.stack(redirect=false, retry=false, aws_authorization=false)

        if request.return_stream && request.response_stream === nothing
            request.response_stream = Base.BufferStream()
        end

        return HTTP.request(
            http_stack,
            request.request_method,
            HTTP.URI(request.url),
            HTTP.mkheaders(request.headers),
            request.content;
            require_ssl_verification=false,
            response_stream=request.response_stream,
            request.http_options...
        )
    catch e
        # Base.IOError is needed because HTTP.jl can often have errors that aren't
        # caught and wrapped in an HTTP.IOError
        # https://github.com/JuliaWeb/HTTP.jl/issues/382
        @delay_retry if isa(e, Sockets.DNSError) ||
                        isa(e, HTTP.ParseError) ||
                        isa(e, HTTP.IOError) ||
                        isa(e, Base.IOError) ||
                        (isa(e, HTTP.StatusError) && _http_status(e) >= 500)
        end
    end
end

"""
    submit_request(aws::AbstractAWSConfig, request::Request; return_headers::Bool=false)

Submit the request to AWS.

# Arguments
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `request::Request`: All the information about making a request to AWS

# Keywords
- `return_headers::Bool=false`: True if you want the headers from the response returned back

# Returns
- `Tuple or Dict`: Tuple if returning_headers, otherwise just return a Dict of the response body
"""
function submit_request(aws::AbstractAWSConfig, request::Request; return_headers::Bool=false)
    response = nothing
    TOO_MANY_REQUESTS = 429
    EXPIRED_ERROR_CODES = ["ExpiredToken", "ExpiredTokenException", "RequestExpired"]
    REDIRECT_ERROR_CODES = [301, 302, 303, 304, 305, 307, 308]
    THROTTLING_ERROR_CODES = [
        "Throttling",
        "ThrottlingException",
        "ThrottledException",
        "RequestThrottledException",
        "TooManyRequestsException",
        "ProvisionedThroughputExceededException",
        "LimitExceededException",
        "RequestThrottled",
        "PriorRequestNotComplete"
    ]

    request.headers["User-Agent"] = user_agent[]
    request.headers["Host"] = HTTP.URI(request.url).host

    @repeat 3 try
        credentials(aws) === nothing || sign!(aws, request)
        response = @mock _http_request(request)

        if response.status in REDIRECT_ERROR_CODES && HTTP.header(response, "Location") != ""
            request.url = HTTP.header(response, "Location")
        end
    catch e
        if e isa HTTP.StatusError
            e = AWSException(e)
        end

        @retry if :message in fieldnames(typeof(e)) && occursin("Signature expired", e.message) end

        # Handle ExpiredToken...
        # https://github.com/aws/aws-sdk-go/blob/v1.31.5/aws/request/retryer.go#L98
        @retry if ecode(e) in EXPIRED_ERROR_CODES
            check_credentials(credentials(aws), force_refresh=true)
        end

        # Throttle handling
        # https://github.com/boto/botocore/blob/1.16.17/botocore/data/_retry.json
        # https://docs.aws.amazon.com/general/latest/gr/api-retries.html
        @delay_retry if e isa AWSException &&
            (_http_status(e.cause) == TOO_MANY_REQUESTS || ecode(e) in THROTTLING_ERROR_CODES)
        end

        # Handle BadDigest error and CRC32 check sum failure
        @retry if e isa AWSException && (
            _header(e.cause, "crc32body") == "x-amz-crc32" ||
            ecode(e) in ("BadDigest", "RequestTimeout", "RequestTimeoutException")
        )
        end

        if e isa AWSException && occursin("Missing Authentication Token", e.message) && aws.credentials === nothing
            return throw(NoCredentials("You're attempting to perform a request without credentials set."))
        end
    end

    response_dict_type = request.response_dict_type

    # For HEAD request, return headers...
    if request.request_method == "HEAD"
        return response_dict_type(response.headers)
    end

    # Return response stream if requested...
    if request.return_stream
        return request.response_stream
    end

    # Return raw data if requested...
    if request.return_raw
        return (return_headers ? (response.body, response.headers) : response.body)
    end

    # Parse response data according to mimetype...
    mime = HTTP.header(response, "Content-Type", "")

    if isempty(mime)
        if length(response.body) > 5 && response.body[1:5] == b"<?xml"
            mime = "text/xml"
        end
    end

    body = String(copy(response.body))

    if occursin(r"/xml", mime)
        xml_dict_type = response_dict_type{Union{Symbol, String}, Any}
        body = parse_xml(body)
        root = XMLDict.root(body.x)

        return (return_headers ? (xml_dict(root, xml_dict_type), response_dict_type(response.headers)) : xml_dict(root, xml_dict_type))
    elseif occursin(r"/x-amz-json-1.[01]$", mime) || endswith(mime, "json")
        info = isempty(response.body) ? nothing : JSON.parse(body, dicttype=response_dict_type)
        return (return_headers ? (info, response_dict_type(response.headers)) : info)
    elseif startswith(mime, "text/")
        return (return_headers ? (body, response_dict_type(response.headers)) : body)
    else
        return (return_headers ? (response.body, response.headers) : response.body)
    end
end

function _generate_rest_resource(request_uri::String, args::AbstractDict{String, <:Any})
    # There maybe a time where both $k and $k+ are in the request_uri, in which case this needs to be updated
    # From looking around, I have not seen an example yet, however that doesn't mean it doesn't exist

    for (k,v) in args
        if occursin("{$k}", request_uri)
            request_uri = replace(request_uri, "{$k}" => v)
        elseif occursin("{$k+}", request_uri)
            request_uri = replace(request_uri, "{$k+}" => HTTP.escapepath(v))
        end
    end

    return request_uri
end

function _return_headers(args::AbstractDict{String, <:Any})
    return_headers = get(args, "return_headers", false)

    if return_headers
        delete!(args, "return_headers")
    end

    return return_headers
end

function _flatten_query(service::String, query::AbstractDict{String, <: Any})
    return _flatten_query!(Pair{String, String}[], service, query)
end

function _flatten_query!(result::Vector{Pair{String, String}}, service::String, query::AbstractDict{String, <:Any}, prefix::String="")
    for (k, v) in query
        if v isa AbstractDict
            _flatten_query!(result, service, v, string(prefix, k, "."))
        elseif v isa AbstractArray
            for (i, j) in enumerate(v)
                suffix = service in ("ec2", "sqs") ? "" : ".member"
                prefix_key = string(prefix, k, suffix, ".", i)

                if j isa AbstractDict
                    _flatten_query!(result, service, j, string(prefix_key, "."))
                else
                    push!(result, Pair(prefix_key, string(j)))
                end
            end
        else
            push!(result, Pair(string(prefix, k), string(v)))
        end
    end

    return result
end

"""
    (service::RestXMLService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestXML request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::RestXMLService)(
    request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    """
        _clean_s3_uri(uri::AbstractString)

    Escape special AWS S3 characters properly.

    AWS S3 allows for various special characters in file names, these characters are not being
    properly escaped before we make the requests.

    We cannot call `HTTP.escapeuri(request.uri)` because this will escape `/` characters which
    are used in the filepathing for sub-directories.

    # Arguments
    - `uri::AbstractString`: URI to to cleaned

    # Returns
    - `String`: URI with characters escaped
    """
    function _clean_s3_uri(uri::AbstractString)
        chars_to_clean = (
            ' ' => "%20",
            '!' => "%21",
            ''' => "%27",
            '(' => "%28",
            ')' => "%29",
            '*' => "%2A",
            '+' => "%2B",
            '=' => "%3D",
        )
        return reduce(replace, chars_to_clean, init=uri)
    end

    request = Request(
        service=service.name,
        api_version=service.api_version,
        request_method=request_method,
        content=get(args, "body", ""),
        headers=LittleDict{String, String}(get(args, "headers", [])),
        return_stream=get(args, "return_stream", false),
        http_options=get(args, "http_options", []),
        response_stream=get(args, "response_stream", nothing),
        return_raw=get(args, "return_raw", false),
    )

    if haskey(args, "response_dict_type")
        request.response_dict_type = pop!(args, "response_dict_type")
    end

    delete!(args, "headers")
    delete!(args, "body")
    return_headers = _return_headers(args)

    if request.service == "s3"
        request_uri = _clean_s3_uri(request_uri)
    end
    request.resource = _generate_rest_resource(request_uri, args)
    query_str = HTTP.escapeuri(args)

    if !isempty(query_str)
        if occursin('?', request.resource)
            request.resource *= "&$query_str"
        else
            request.resource *= "?$query_str"
        end
    end

    request.url = generate_service_url(aws_config, request.service, request.resource)

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::QueryService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a Query request to AWS.

# Arguments
- `operation::String`:
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::QueryService)(
    operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    POST_RESOURCE = "/"
    return_headers = _return_headers(args)

    request = Request(
        service=service.name,
        api_version=service.api_version,
        resource=POST_RESOURCE,
        request_method="POST",
        headers=LittleDict{String, String}(get(args, "headers", [])),
        url=generate_service_url(aws_config, service.name, POST_RESOURCE),
        return_stream=get(args, "return_stream", false),
        http_options=get(args, "http_options", []),
        response_stream=get(args, "response_stream", nothing),
        return_raw=get(args, "return_raw", false),
    )

    if haskey(args, "response_dict_type")
        request.response_dict_type = pop!(args, "response_dict_type")
    end

    request.headers["Content-Type"] = "application/x-www-form-urlencoded; charset=utf-8"

    args["Action"] = operation
    args["Version"] = service.api_version
    request.content = HTTP.escapeuri(_flatten_query(service.name, args))

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::JSONService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a JSON request to AWS.

# Arguments
- `operation::String`: Name of the operation to perform
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::JSONService)(
    operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    POST_RESOURCE = "/"
    return_headers = _return_headers(args)

    request = Request(
        service=service.name,
        api_version=service.api_version,
        resource=POST_RESOURCE,
        request_method="POST",
        headers=LittleDict{String, String}(get(args, "headers", [])),
        content=json(args),
        url=generate_service_url(aws_config, service.name, POST_RESOURCE),
        return_stream=get(args, "return_stream", false),
        http_options=get(args, "http_options", []),
        response_stream=get(args, "response_stream", nothing),
        return_raw=get(args, "return_raw", false),
    )

    if haskey(args, "response_dict_type")
        request.response_dict_type = pop!(args, "response_dict_type")
    end

    request.headers["Content-Type"] = "application/x-amz-json-$(service.json_version)"
    request.headers["X-Amz-Target"] = "$(service.target).$(operation)"

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::RestJSONService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestJSON request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::RestJSONService)(
    request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return_headers = _return_headers(args)

    request = Request(
        service=service.name,
        api_version=service.api_version,
        request_method=request_method,
        headers=LittleDict{String, String}(get(args, "headers", [])),
        resource=_generate_rest_resource(request_uri, args),
        return_stream=get(args, "return_stream", false),
        http_options=get(args, "http_options", []),
        response_stream=get(args, "response_stream", nothing),
        return_raw=get(args, "return_raw", false),
    )

    if haskey(args, "response_dict_type")
        request.response_dict_type = pop!(args, "response_dict_type")
    end

    request.url = generate_service_url(aws_config, request.service, request.resource)

    if !isempty(service.service_specific_headers)
        merge!(request.headers, service.service_specific_headers)
    end

    request.headers["Content-Type"] = "application/json"
    delete!(args, "headers")
    request.content = json(args)

    return submit_request(aws_config, request; return_headers=return_headers)
end

end  # module AWS


text:
using CourseMatch
using Base.Test

@testset "demand function tests" begin
    @testset "unimodal preferences 1 student" begin
	    p = rand(3)
	    pref = [sparse(1:3,1:3,[100;0;0])] #one class prefered
	    budg = [100] #arbitrary budget
	    req = [1] #1 class is required
	    solution = CourseMatch.demand(p, pref, budg, req)
	    @test solution[:ind_demands] == [1 0 0]
	    @test solution[:course_demand] == [1 0 0]
	    @test solution[:total] == 1
    end

    @testset "unimodal preferences 2 students" begin
	    p = rand(3)
	    pref = []
	    push!(pref,sparse(1:3,1:3,[100;0;0]))
	    push!(pref,sparse(1:3,1:3,[0;100;0]))
	    budg = [100;100] #arbitrary budget
	    req = [1;1] #1 class is required
	    solution = CourseMatch.demand(p, pref, budg, req)
	    @test solution[:ind_demands] == [1 0 0;0 1 0]
	    @test solution[:course_demand] == [1 1 0]
	    @test solution[:total] == 2
    end

    @testset "unaffordable most preferred course" begin
	    p = rand(3)
	    p[1] = 101  # high price on course 1
	    pref = [sparse(1:3,1:3,[100;50;0])] #one class prefered
	    budg = [100] #arbitrary budget
	    req = [1] #1 class is required
	    solution = CourseMatch.demand(p, pref, budg, req)
	    @test solution[:ind_demands] == [0 1 0]
	    @test solution[:course_demand] == [0 1 0]
	    @test solution[:total] == 1
 	end

    @testset "unaffordable most preferred course 2 students" begin
	    p = rand(3)
	    p[1] = 101  # high price on course 1
	    pref = []
	    push!(pref,sparse(1:3,1:3,[100;50;0]))
	    push!(pref,sparse(1:3,1:3,[100;0;50]))
	    budg = [100;100] #arbitrary budget
	    req = [1;1] #1 class is required
	    solution = CourseMatch.demand(p, pref, budg, req)
	    @test solution[:ind_demands] == [0 1 0;0 0 1]
	    @test solution[:course_demand] == [0 1 1]
	    @test solution[:total] == 2
 	end

    @testset "substitutable classes" begin
    	#Here I have one student, 3 available classes, and it is required to enroll in two classes in total. The student as a stric ordering on the three classes
    	#but they have a strong preference over not taking their two most preferred ones together such that they prefer o take their preferred one
    	#and their least preferred one. The test checks if the demand functions does indeed select
    	#the most preferred combination of courses in this context.
    	# p = rand(3)
    	# pref = [[100 -100 0 ; -100 50 0; 0 0 0]]
    	# budg = [100] #arbitrary budget
    	# req = [2] #2 classes are required
    	# solution = CourseMatch.demand(p, pref, budg, req)
    	# @test solution[:ind_demands][1] == [1, 0, 1]
    	# @test solution[:total] == [1, 0, 1]
 	end
    @testset "identical students" begin
	    #Here I have two strictly identical students, 3 available classes, and it is required to enroll in one classes n total. The test checks if the demand
	    #functions does indeed give the same demand for the two students.
	    # p = rand(3)
	    # pref = Matrix(diagm(rand(0:100, 3)))
	    # pref = [pref, pref]
	    # budg = rand(150:200, 1)[1]
	    # budg = [budg, budg]
	    # req = [1, 1] #1 class is required
	    # solution = CourseMatch.demand(p, pref, budg, req)
	    # @test solution[:ind_demands][1] == solution[:ind_demands][2]
	    # @test solution[:total] == 2*solution[:ind_demands][2]
 	end
end

ext: jl

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
function ShotProfileLSEWEM(m::Array{AbstractString,1},m0::Array{AbstractString,1},d::Array{AbstractString,1},param=Dict())

	# least squares shot profile wave equation migration of isotropic 3C data.
	cost = get(param,"cost","cost.txt")   # cost output text file
	precon = get(param,"precon",false)   # flag for preconditioning by smoothing the image
	wd = join(["tmp_LSM_wd_",string(int(rand()*100000))])
	CalculateSampling(d[1],wd,param)
	param["wd"] = wd
	param["tmute"] = 0.;
	param["vmute"] = 999999.;
	if (precon == true)
		param["operators"] = [ApplyDataWeights SeisMute ShotProfileEWEM SeisSmoothGathers]
	else
		param["operators"] = [ApplyDataWeights SeisMute ShotProfileEWEM]
	end
	ConjugateGradients(m,m0,d,cost,param)
	SeisRemove(wd)

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# Passing Julia callback functions to Python routines.
#
# Note that this will typically involve two functions: the
# desired Julia function/closure, and a top-level C-callable
# wrapper function used with PyCFunction_NewEx -- the latter
# is called from Python and calls the former as needed.

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

# Define a Python method/function object from f(PyPtr,PyPtr)::PyPtr.
# Requires f to be a top-level function.
function pymethod(f::Function, name::AbstractString, flags::Integer)
    # Python expects the PyMethodDef structure to be a *constant*,
    # so we define an anonymous global to hold it.
    def = gensym("PyMethodDef")
    @eval const $def = PyMethodDef[PyMethodDef($name, $f, $flags)]
    PyObject(@pycheckn ccall((@pysym :PyCFunction_NewEx), PyPtr,
                             (Ptr{PyMethodDef}, Ptr{Void}, Ptr{Void}),
                             eval(def), C_NULL, C_NULL))
end

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

# To pass an arbitrary Julia object to Python, we wrap it
# in a jlwrap Python class, where jlwrap.__call__
# executes the pyjlwrap_call function in Julia, which
# in turn fetches the actual Julia value from the python wrapper
# and calls it.

function pyjlwrap_call(self_::PyPtr, args_::PyPtr, kw_::PyPtr)
    args = PyObject(args_) # don't need pyincref because of finally clause below
    try
        f = unsafe_pyjlwrap_to_objref(self_)

        # on 0.6 we need to use invokelatest to get execution in newest world
        @static if isdefined(Base, :invokelatest)
            if kw_ == C_NULL
                ret = PyObject(Base.invokelatest(f, convert(PyAny, args)...))
            else
                kw = PyDict{Symbol,PyAny}(pyincref(kw_))
                kwargs = [ (k,v) for (k,v) in kw ]

                # 0.6 `invokelatest` doesn't support kwargs, instead
                # use a closure over kwargs. see:
                #   https://github.com/JuliaLang/julia/pull/22646
                f_kw_closure() = f(convert(PyAny, args)...; kwargs...)
                ret = PyObject(Core._apply_latest(f_kw_closure))
            end
        else # 0.5 support
            if kw_ == C_NULL
                ret = PyObject(f(convert(PyAny, args)...))
            else
                kw = PyDict{Symbol,PyAny}(pyincref(kw_))
                kwargs = [ (k,v) for (k,v) in kw ]
                ret = PyObject(f(convert(PyAny, args)...; kwargs...))
            end
        end

        return pystealref!(ret)
    catch e
        pyraise(e)
    finally
        args.o = PyPtr_NULL # don't decref
    end
    return PyPtr_NULL
end


text:
# Lee, E. C., Kelly, M. R., Ochocki, B. M., Akinwumi, S. M., Hamre, K. E., Tien, J. H., Eisenberg, M. C.,
# Model distinguishability and inference robustness in mechanisms of cholera transmission and loss of immunity
# https://doi.org/10.1016/j.jtbi.2017.01.032
# Eq. (3)
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
    I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
    W'(t) = xi * (I(t) - W(t)),
    R'(t) = gam * I(t) - (mu + a) * R(t),
    y(t) = k * I(t)
)

@time println(assess_identifiability(ode))


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")
fun_number = 2
@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data/meshs_ks_func_count_$(fun_number).csv",DataFrame)
ks_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_ks_func_count_$(fun_number).csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error between true CDF and reconstruction")
# plot!(legend=(0.7,0.8))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data/meshs_l2_pdf_func_count_$(fun_number).csv",DataFrame)
l2_pdf_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_l2_pdf_func_count_$(fun_number).csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross,:diamond]
for (c,col) in enumerate(names(l2_pdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
# plot!(legend=:outerbottomright)
plot!()
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l2_pdf_error_formatted.pdf")
end
p()

function p()
l1_cdf_data = CSV.read((@__DIR__)*"/../data/meshs_l1_cdf_func_count_$(fun_number).csv",DataFrame)
l1_cdf_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_l1_cdf_func_count_$(fun_number).csv",DataFrame).DG_limiter

@halfwidth_plot_defaults()
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross,:diamond]
for (c,col) in enumerate(names(l1_cdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l1_cdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.225,0.625))
plot!()
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l1_cdf_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()
dir = @__DIR__()
os = 9:2:15
plot(layout=(5,length(os)))
linestyles_vec = [:solid,:solid,:solid,:solid,:solid]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(2+eps(),7-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c+4,1]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,cell_probs(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
            linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[1],
            linewidth=2,xticks=false)
    end
    plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "DG" : false), grid=false,
        title = string("Dim. ",o))

    mesh = DGMesh(0.0:(1.0./o):10.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c+4,2]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,cell_probs(dt).(xs,1),subplot=c+4,label=(c==1 ? "Unif." : false),
            linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[2],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "Unif." : false), grid=false)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c+4,3]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,cell_probs(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG lim" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG lim" : false), grid=false)

    mesh = FRAPMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c+4,4]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,cell_probs(dt).(xs,1),subplot=c+12,label=(c==1 ? "Q-R" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false, grid=false)
    end
    (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "Q-R" : false), grid=false)
    # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
    #     linestyle=:dot,
    #     # marker=markerstyles_vec[3],
    #     color=:grey,
    #     linewidth=2)

    # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
    #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
    #     ylabel=(c==1 ? "probability" : false), grid=false,
    #     title = string("Dim. ",o),
    #     legend=:outertopright)
    plot!(legend=false)

    model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = interior_point_mass(eps(),1,dq)
    dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)

    plot!(xs,pdf(dt).(xs,1),subplot=c+16,label=(c==1 ? "DG lin lim" : false),
        linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
        # marker=markerstyles_vec[3],
        linewidth=2,xticks=0:2:10, grid=false)
    (c==1)&&plot!(subplot=c+16, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG lin lim" : false), grid=false)
end
plot!()
savefig((@__DIR__)*"/cell_probs_formatted.pdf")
end
p()

function p()
@halfwidth_plot_defaults()
L1_cell_probs_errors = DataFrame(DG = Float64[], Unif = Float64[], DG_limiter = Float64[], QBDRAP = Float64[], DG_lin_lim = Float64[])
truth = zeros(10)
truth[5] = 1.0
truth_pm = zeros(3)
for (c_o,o) in enumerate(1:2:21)
    row = zeros(5)
    for (c_m,m) in enumerate(names(L1_cell_probs_errors))
        (c_m<5)&&jldopen((@__DIR__)*"/../../../coeffs_matrix.jld2") do f
            model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

            mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
            dq = DiscretisedFluidQueue(model,mesh)

            coeffs = f["coeff_matrix"][fun_number,c_o,c_m]
            dt = SFMDistribution(coeffs,dq)

            x_vals = range(0.5,9.5;length=10)

            approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
        (c_m==5)&&begin
            model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

            mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
            dq = DiscretisedFluidQueue(model,mesh)
            d0 = interior_point_mass(0.5,1,dq)
            global dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)
            coeffs = dt.coeffs

            x_vals = (mesh.nodes[1:end-1]+mesh.nodes[2:end])./2

            approx = sum(reshape(cell_probs(dt).(x_vals,1),Int((o+1)/2),10),dims=1)[:]
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
    end
    push!(L1_cell_probs_errors,row)
end

file = (@__DIR__)*"/../data/meshs_l1_cell_"*"func_count_$(fun_number)"
CSV.write(file*".csv",L1_cell_probs_errors)
q = plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
        linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension", ylabel="Error",
    title="Error between cell probabilities",)
    # legend=(0.225,0.625))
error_ticks!(plot!())
lines = @add_lines!(L1_cell_probs_errors,("Unif","QBDRAP"),@__DIR__)

plot!(
    lines["DG"].x,lines["DG"].y,linestyle=:dot,label=false,
    color=:black,linewidth=2,
)
annotate!((lines["DG"].x[end]+0.09,lines["DG"].y[end],string(round(lines["DG"].slope,digits=2))))
savefig((@__DIR__)*"/L1_cell_probs.pdf")
end
p()

# os = 1:2:7
# plot(layout=(1,length(os)))
# linestyles_vec = [:solid,:dash,:dashdot]
# markerstyles_vec = [:cross,:diamond,:circle]
# xs = range(eps(),1-eps(),length=200)
# for (c,o) in enumerate(os)
#     model = BoundedFluidQueue(fill(0.0,1,1),[1.0],1.0)

#     mesh = DGMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
#         linestyle=linestyles_vec[1],
#         # marker=markerstyles_vec[1],
#         linewidth=2)

#     mesh = DGMesh(0.0:(1.0./o):1.0,1)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "Unif." : false),
#         linestyle=linestyles_vec[2],
#         # marker=markerstyles_vec[2],
#         linewidth=2)

#     mesh = FRAPMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "QBD-RAP" : false),
#         linestyle=linestyles_vec[3],
#         # marker=markerstyles_vec[3],
#         linewidth=2)

#     plot!(xs,x->2*Float64(x<=0.5),subplot=c,label=(c==1 ? "truth" : false),
#         linestyle=:dot,
#         # marker=markerstyles_vec[3],
#         color=:grey,
#         linewidth=2)

#     plot!(subplot=c, yticks=(c==1 ? (-2:2:6) : false),
#         ylims = (-3,7.2), xticks = 0.25:0.5:0.75,
#         ylabel=(c==1 ? "probability" : false), grid=false,
#         title = string("Dimension ",o),
#         legend=(0.1,0.8))
# end
# plot!()
# savefig((@__DIR__)*"/pdfs_formatted.pdf")


text:
using Immerse, Colors

hfig = figure()
x = linspace(0,4pi,101)
display(plot(x=x, y=sin(x), Geom.line(tag=:line)))
val = setproperty!((hfig,:line), rand(1:5), :linewidth)
@test getproperty((hfig,:line), :linewidth) == [val*mm]
val = setproperty!((hfig,:line), RGB(rand(),rand(),rand()), :stroke)
@test getproperty((hfig,:line), :stroke) == [coloralpha(val)]
setproperty!((hfig,:line), false, :visible)
@test getproperty((hfig,:line), :visible) == [false]
setproperty!((hfig,:line), true, :visible)
@test getproperty((hfig,:line), :visible) == [true]


text:
module LimitOrderBook
    import Base.@kwdef
    include("orderqueue.jl")
    include("sidebook.jl")
    include("book.jl")
    export OrderBook, Order, submit_limit_order!, cancel_limit_order!, submit_market_order!,
        book_depth_info, volume_bid_ask, best_bid_ask, n_orders_bid_ask,
        get_acct, write_csv
end


text:
function free(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintFree), Void, (Ptr{Void},), c.ptr)
end

function get_space(c::Constraint)
    Space(ccall(dlsym(libchipmunk, :cpConstraintGetSpace), Ptr{Void}, (Ptr{Void},), c.ptr))
end

function get_bodyA(c::Constraint)
    Body(ccall(dlsym(libchipmunk, :cpConstraintGetBodyA), Ptr{Void}, (Ptr{Void},), c.ptr))
end

function get_bodyB(c::Constraint)
    Body(ccall(dlsym(libchipmunk, :cpConstraintGetBodyB), Ptr{Void}, (Ptr{Void},), c.ptr))
end

function get_max_force(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetMaxForce), Cdouble, (Ptr{Void},), c.ptr)
end

function set_max_force(c::Constraint, max_force::Real)
    ccall(dlsym(libchipmunk, :cpConstraintSetMaxForce), Void, (Ptr{Void}, Cdouble,), c.ptr, max_force)
end

function get_error_bias(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetErrorBias), Cdouble, (Ptr{Void},), c.ptr)
end

function set_error_bias(c::Constraint, error_bias::Real)
    ccall(dlsym(libchipmunk, :cpConstraintSetErrorBias), Void, (Ptr{Void}, Cdouble,), c.ptr, error_bias)
end

function get_max_bias(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetMaxBias), Cdouble, (Ptr{Void},), c.ptr)
end

function set_max_bias(c::Constraint, max_bias::Real)
    ccall(dlsym(libchipmunk, :cpConstraintSetMaxBias), Void, (Ptr{Void}, Cdouble,), c.ptr, max_bias)
end

function get_collide_bodies(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetCollideBodies), Bool, (Ptr{Void},), c.ptr)
end

function set_collide_bodies(c::Constraint, collide_bodies::Bool)
    ccall(dlsym(libchipmunk, :cpConstraintSetCollideBodies), Void, (Ptr{Void}, Int32,), c.ptr, collide_bodies)
end

function get_userdata(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetUserData), Ptr{Void}, (Ptr{Void},), c.ptr)
end

function set_userdata(c::Constraint, userdata::Ptr{Void})
    ccall(dlsym(libchipmunk, :cpConstraintSetUserData), Void, (Ptr{Void}, Ptr{Void},), c.ptr, userdata)
end

function get_impulse(c::Constraint)
    ccall(dlsym(libchipmunk, :cpConstraintGetImpulse), Cdouble, (Ptr{Void},), c.ptr)
end

export free, get_space, get_bodyA, get_bodyB, get_max_force, set_max_force, get_error_bias,
set_error_bias, get_max_bias, set_max_bias, get_collide_bodies, set_collide_bodies, get_userdata,
set_userdata, get_impulse


text:
using Base64

# Compatibility for data backends which implement DataSets 0.2.0 storage
# interface.

#-------------------------------------------------------------------------------
struct OldBackendAPI
    data
end

function _lookup_path(storage::OldBackendAPI, path)
    # For backends which try to avoid depending on DataSets types
    path = string(path)
    isempty(path) ? storage.data : get(storage.data, path, nothing)
end

function _get_data(storage::OldBackendAPI, path)
    base64decode(_lookup_path(storage, path))
end

function Base.open(f::Function, ::Type{IO}, storage::OldBackendAPI, path; kws...) where {T}
    f(IOBuffer(_get_data(storage, path)))
end

function Base.read(storage::OldBackendAPI, path, ::Type{T}) where {T}
    read(IOBuffer(_get_data(storage, path)), T)
end

function Base.read(storage::OldBackendAPI, path)
    _get_data(storage, path)
end

function Base.readdir(storage::OldBackendAPI, path)
    if isempty(path)
        sort!(collect(keys(storage.data)))
    else
        []
    end
end

function Base.isdir(storage::OldBackendAPI, path)
    path = string(path)
    @assert storage.data isa Dict
    isempty(path)
end

function Base.isfile(storage::OldBackendAPI, path)
    _lookup_path(storage, path) isa String
end

function Base.ispath(storage::OldBackendAPI, path)
    !isnothing(_lookup_path(storage, path))
end

function connect_old_backend(f, config, ds)
    storage = OldBackendAPI(config["data"])
    if config["type"] == "Blob"
        f(Blob(storage))
    else
        f(BlobTree(storage))
    end
end

DataSets.add_storage_driver("OldBackendAPI"=>connect_old_backend)

#-------------------------------------------------------------------------------
@testset "OldBackendAPI" begin
    proj = DataSets.load_project("Data.toml")

    @test open(IO, dataset(proj, "old_backend_blob")) do io
           read(io, String)
    end == "x"
    @test String(open(read, IO, dataset(proj, "old_backend_blob"))) == "x"
    @test open(Vector{UInt8}, dataset(proj, "old_backend_blob")) == UInt8['x']
    @test read(open(dataset(proj, "old_backend_blob")), String) == "x"
    @test read(open(dataset(proj, "old_backend_blob"))) == UInt8['x']

    @test readdir(open(dataset(proj, "old_backend_tree"))) == ["a.txt", "b.txt"]
    @test open(dataset(proj, "old_backend_tree"))[path"a.txt"] isa Blob
    @test read(open(dataset(proj, "old_backend_tree"))[path"a.txt"], String) == "a"
    @test read(open(dataset(proj, "old_backend_tree"))[path"b.txt"], String) == "b"
end


text:
# Implement DistanceBetweenPatternAndStrings

calc_hamming_distance(s1, s2) = count(((c1, c2),) -> c1 != c2, zip(s1, s2))

function distance_between_pattern_and_text(text, pattern, k)
    mindist = k
    for ix in 1:(length(text)-k+1)
        kmer = text[ix:ix+k-1]
        distance = calc_hamming_distance(pattern, kmer)
        if distance < mindist
            mindist = distance
        end
    end
    mindist
end

distance_between_pattern_and_collection_of_texts(texts, pattern, k) =
    foldl((distance, text) -> distance + distance_between_pattern_and_text(text, pattern, k), texts; init=0)

function main()
    pattern = readline()
    k = length(pattern)
    texts = split(readline())
    result = distance_between_pattern_and_collection_of_texts(texts, pattern, k)
    println(result)
end

main()


text:
using Random; Random.seed!(1234);

include("src/data.jl")
using .DataUtil
include("src/stepwise/SelectiveInference.jl")
using .SelectiveInference
include("src/selective_p.jl")
using .MultiTruncatedDistributions
using LinearAlgebra:I
using Distributions
using IntervalSets
using Parameters
using ProgressMeter
using Plots

function selective_p_rand(test_prob::SelectiveHypothesisTest; isOC=false)
    @unpack alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺 = test_prob
    M_obs = collect(alg_result.model.M)
    isempty(M_obs) && return NaN
    test_stat = test_stats[rand(1:length(M_obs))]
    z_obs, σ²_z = test_stat.z_obs, test_stat.σ²_z

    # compute Z = {z | selection_event}
    z_min = -abs(z_obs) - 10 * √(σ²_z); z_max = abs(z_obs) + 10 * √(σ²_z)
    Z = isOC ? SelectiveInference.region(𝐗, 𝐲_obs, 𝚺, test_stat, alg)[2] ∩ [z_min..z_max] : SelectiveInference.truncated_interval(alg_result.model, 𝐗, 𝚺, test_stat, alg; z_min=z_min, z_max=z_max)

    𝜋 = cdf(TruncatedDistribution(Normal(0, √(σ²_z)), Z), z_obs)
    return 2 * min(𝜋, 1 - 𝜋)
end

# FPR experiment
function FPR_experiment(param::DataParam, alg::StepwiseFeatureSelection;
    isOC=false,
    α=0.05,
    denominator_upper=100,
    iter_size=5)
    FPR = fill(NaN, iter_size)
    @showprogress for l = 1:iter_size
        numerator = 0
        denominator = 0
        for _ = 1:denominator_upper
            try
                𝐗, 𝐲_obs = make_dataset(param)
                alg_result = stepwise(𝐗, 𝐲_obs, I, alg)
                p_value = selective_p_rand(SelectiveHypothesisTest(alg, alg_result, 𝐗, 𝐲_obs, I); isOC=isOC)
                isnan(p_value) && continue
                numerator += p_value < α ? 1 : 0
                denominator += 1
            catch e
println(e)
                continue
            end
        end
        FPR[l] = numerator / denominator
    end
    return FPR
end

using JLD
function main(alg::StepwiseFeatureSelection)
    # parameters
    β₀ = 0.0
    Ns = (50, 100, 150)
    ps = (10, 20, 50)

    iter_size = 20
        FPR = Array{Float64}(undef, iter_size, length(Ns), length(ps))
    FPR_OC = similar(FPR)
    α = 0.05
    for (i, N) in enumerate(Ns)
        for (j, p) in enumerate(ps)
            FPR[:,i,j] = FPR_experiment(DataParam(N=N, β₀=β₀, 𝛃=zeros(p)), alg; isOC=false, α=α, iter_size=iter_size)
            FPR_OC[:,i,j] = FPR_experiment(DataParam(N=N, β₀=β₀, 𝛃=zeros(p)), alg; isOC=true, α=α, iter_size=iter_size)
        end
    end

    FPR_mean = reshape(mean(FPR, dims=1), size(FPR, 2), size(FPR, 3))
    FPR_OC_mean = reshape(mean(FPR_OC, dims=1), size(FPR_OC, 2), size(FPR_OC, 3))

    h_cs = (:red2, :violetred2, :darkred)
    o_cs = (:blue2, :skyblue, :navyblue)
    p = plot()
    for j = 1:size(FPR, 3) # σ
        plot!(p, 1:length(Ns), FPR_mean[:,j], label="Homotopy, p=$(ps[j])", color=h_cs[j], markershape=:circle)
        plot!(p, 1:length(Ns), FPR_OC_mean[:,j], label="Quadratic, p=$(ps[j])", color=o_cs[j], markershape=:rect)
    end
    # plot setting
    plot!(p; ylims=(0, 0.3))
    plot!(p; xlabel="Sample size", ylabel="False Positive Rate (FPR)")
    plot!(p; legend=:topright)
    plot!(p; guidefontsize=16, legendfontsize=12, tickfontsize=10)
    plot!(p; xticks=(1:length(Ns), Ns))
    savefig(p, "img/FPR_forward_backward.pdf")
end


text:
using ImGuiMakieBackend
using Test

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


text:
using Unitful: @u_str
using UnitfulAtomic: austrip

"""
    ASEFrictionProvider{A} <: ElectronicFrictionProvider

Obtain the electronic friction from an ASE calculator that implements `get_friction_tensor`.
Assumes that the units of friction are "eV/Å/Å".
Construct by passing the ase atoms object with the calculator already attached.
"""
struct ASEFrictionProvider{A} <: ElectronicFrictionProvider
    atoms::A
end

NQCModels.ndofs(::ASEFrictionProvider) = 3

function friction!(model::ASEFrictionProvider, F::AbstractMatrix, R::AbstractMatrix)
    set_coordinates!(model, R)
    F .= model.atoms.get_friction_tensor()
    @. F = austrip(F * u"eV/Å/Å")
end

function set_coordinates!(model::ASEFrictionProvider, R)
    model.atoms.set_positions(ustrip.(auconvert.(u"Å", R')))
end


text:
module Reader

using ..Client: Connection, Core.read_one

import ..Wrapper

include("decoder.jl")

function read_msg(socket)

  msg = read_one(socket)

  # Assert and chop last null char
  @assert length(msg) > 1 && msg[1] != pop!(msg) == 0x00

  split(String(msg), '\0')
end

"""
    check_msg(ib, wrap)

Process one message and dispatch the appropriate callback. **Blocking**.
"""
function check_msg(ib::Connection, w::Wrapper)

  msg = read_msg(ib.socket)

  Decoder.decode(msg, w, ib.version)
end

"""
    check_all(ib, wrap, flush=false)

Process all messages waiting in the queue. **Not blocking**.
If `flush=true`, messages are read but callbacks are not dispatched.

Return number of messages processed.
"""
function check_all(ib::Connection, w::Wrapper, flush::Bool=false)

  count = 0
  while bytesavailable(ib.socket) > 0 || ib.socket.status == Base.StatusOpen # =3

    msg = read_msg(ib.socket)

    flush || Decoder.decode(msg, w, ib.version)

    count += 1
  end

  count
end

"""
    start_reader(ib, wrap)

Start a new [`Task`](@ref) to process messages asynchronously.
"""
function start_reader(ib::Connection, w::Wrapper)

  @async  begin
            try
              while true
                check_msg(ib, w)
              end

            catch e

              e isa EOFError && @warn "Connection Terminated"
              println(e)
            end

            @info "Reader Exiting"
          end
end

end


text:
export ConvDims

"""
    ConvDims

Type system-level information about convolution dimensions. Critical for things like
`im2col!()` to generate efficient code, and helpful to reduce the number of kwargs
getting passed around.

We don't want to specialize on things like image size/channel count, so we generally
store those as fields, just for convenience, and to allow for non-breaking changes when
we decide we _do_ want to specialize on those values.  We always want to specialize on
things like stride, padding, dilation, and kernel flipping though.
"""
abstract type ConvDims{N, S, P, D, F} end

# Hack to get rid of type parameters
function basetype(::Type{C}) where {C <: ConvDims}
    if C <: DepthwiseConvDims
        return DepthwiseConvDims
    elseif C <: DenseConvDims
        return DenseConvDims
    elseif C <: PoolDims
        return PoolDims
    else
        return nothing
    end
end

# Obvious getter definitions for the type system-level definitions
spatial_dims(c::ConvDims{N,S,P,D,F}) where {N, S, P, D, F} = N
stride(c::ConvDims{N,S,P,D,F}) where {N, S, P, D, F} = S
padding(c::ConvDims{N,S,P,D,F}) where {N, S, P, D, F} = P
dilation(c::ConvDims{N,S,P,D,F}) where {N, S, P, D, F} = D
flipkernel(c::ConvDims{N,S,P,D,F}) where {N, S, P, D, F} = F

"""
    im2col_dims(c::ConvDims)

im2col calculates, for each output pixel, the "convolution" of N kernels where N is the
number of output channels, by doing a matrix multiply.  The dimensions of that matrix
are given by this function.

Note that because im2col is multithreaded, we need to allocate a separate workspace of
memory per-thread; hence the dimensions returned by this will depend on the number of
threads Julia is currently running with.
"""
function im2col_dims(c::ConvDims)
    return (
        # Output size
        prod(output_size(c)),
        # Size of single dotproduct within convolution
        prod(kernel_size(c))*channels_in(c),
        # One workspace per thread
        Threads.nthreads(),
    )
end

# Protect your skin, kids.  Also do common validation of stride, padding, etc...
function check_spdf(x_size::NTuple{N}, w_size::NTuple{N}, stride, padding, dilation) where {N}
    # Number of spatial dimensions in `x` and `w`.
    nd = N - 2

    # Given a number, duplicate it out to have `nd` length.  If it's already a collection,
    # just splat it out into a tuple so it's always a tuple.  We'll lint length later.
    expand_size(p::Number) = ntuple(_ -> Int(p), nd)
    expand_size(p) = tuple(p...)

    # Convert stride, padding, dilation, etc.. to fully-specified tuples
    pstride = expand_size(stride)
    pdilation = expand_size(dilation)
    ppadding = expand_size(padding)

    if length(pstride) != nd
        throw(DimensionMismatch("Stride $(length(stride))d, should be $(nd)d!"))
    end
    if length(pdilation) != nd
        throw(DimensionMismatch("Dilation $(length(pdilation))d, should be $(nd)d!"))
    end

    # padding is kind of a special case; we allow it to be either 2-length or 4-length,
    # since we support asymmetrical padding
    if length(ppadding) != 2*nd
        if length(ppadding) == nd
            # Do this repeat dance so that we get lo/hi symmetrical padding
            ppadding = tuple(repeat(collect(ppadding), inner=2)...)
        else
            throw(DimensionMismatch("Padding $(length(ppadding))d, should be either $(nd)d or $(2*nd)d!"))
        end
    end

    # Assert that kernel size * dilation is <= padded input size
    for idx in 1:nd
        Is = x_size[idx]
        Pl = ppadding[(idx - 1)*2 + 1]
        Ph = ppadding[(idx - 1)*2 + 2]
        Ks = w_size[idx]
        Ds = pdilation[idx]
        if Is + Pl + Ph < (Ks - 1)*Ds + 1
            throw(DimensionMismatch("Kernel * dilation (($Ks - 1) * $Ds + 1) cannot be larger than input + padding ($Is + $Pl + $Ph)!"))
        end
    end

    return pstride, ppadding, pdilation
end

"""
    output_size(c::ConvDims)

Calculate the output (spatial) dimensions of the convolution.  Get channel count via
`channels_out(c)`, and batch count is unknowable.
"""
function output_size(c::ConvDims)
    I = input_size(c)
    K = kernel_size(c)
    S = stride(c)
    P = padding(c)
    D = dilation(c)

    return ntuple(spatial_dims(c)) do i
        return div(I[i] + P[(i-1)*2 + 1] + P[(i-1)*2 + 2] - (K[i] - 1) * D[i] - 1, S[i]) + 1
    end
end

# Override show() for these beauties
function Base.show(io::IO, cdims::C) where {C <: ConvDims}
    I = (input_size(cdims)..., channels_in(cdims))
    O = (output_size(cdims)..., channels_out(cdims))
    K = kernel_size(cdims)
    S = stride(cdims)
    P = padding(cdims)
    D = dilation(cdims)
    F = flipkernel(cdims)
    print(io, "$(basetype(C)): $I * $K -> $O, stride: $S, pad: $P, dil: $D, flip: $F")
end


text:
timestwo(a) = (a ≈ 0) ? 0 : a+timestwo(a/2)

timestwo(1234.12341)

using Plots

gr()
plot(2.^(1./(1:200)),[((i*2 - timestwo(i))^2)/i for i in 2.^(1./(1:200))],
     label = "", size = (1000, 500))
ylabel!("Percentage error")
xlabel!("x")
title!(string("Error of recursive timestwo as fractions become small"))
xflip!()
savefig("timestwo-recursive")

max_stars_repo_name: vankesteren/julia-tryout

max_issues_repo_name: vankesteren/julia-tryout

max_forks_repo_name: vankesteren/julia-tryout


text:
using PascalsTriangle
const PT = PascalsTriangle
using Test

@testset "PascalsTriangle.jl" begin

    include("entries.jl")
    include("zerorange.jl")
    include("rows.jl")
    include("columns.jl")
    include("centres.jl")

end


text:
using StatsBase

function randhands!(a::Matrix{UInt8})
    for i in eachindex(a)
        a[i] = rand(0x0:(ranks-0x1))
    end
    sort!(a, dims = 2)
    return a
end

function atleast(cards, ranks, k, n)
    s = 0
    a = Vector{Int}(undef, cards)
    for j in 1:n
        for i in eachindex(a)
            a[i] = rand(1:ranks)
        end
        cs = counts(a)
        if maximum(cs) >= k
            s += 1
        end
    end
    s / n
end

function exactly(cards, ranks, k, n)
    s = 0
    a = Vector{Int}(undef, cards)
    for j in 1:n
        for i in eachindex(a)
            a[i] = rand(1:ranks)
        end
        cs = counts(a)
        if maximum(cs) == k
            s += 1
        end
    end
    s / n
end

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
"""
plot_rtm(image, spacing)

Plots an rtm image. Image needs to be (x, z).

# Arguments

- `image`: RTM image shape (nx, nz)
- `spacing`: Grid spacing (dx, dz)
- `perc`: perc for vmin/vmax, default at 98
- `cmap`: Colormap, default is "PuOr"
- `o`: Origin of image, default is (0, 0)
- `interp`: matplotlib interplolation, default is "hanning"
- `aspect`: Aspect ratio of plot, defaults to "auto"
- `d_scale`: Depth scaling power, depth^d_scale, default is 1.5
- `name`: Plot title, default is "RTM"
- `units`: Unit of spacing, default is "m"
- `new_fig`: Whether to create a new figure, default is true
- `save`: Whether to save the figure, default name is "name" if just true, or save if save is a string
"""

function plot_rtm(image, spacing; perc=98, cmap="PuOr", o=(0, 0),
                interp="hanning", aspect=nothing, d_scale=1.5,
                name="RTM", units="m", new_fig=true, save=nothing)
nx, nz = size(image)
dx, dz = spacing
ox, oz = o
depth = range(oz, oz + (nz - 1)*spacing[2], length=nz).^d_scale

scaled = image .* depth'

a = quantile(abs.(vec(scaled)), perc/100)
extent = [ox, ox+ (nx-1)*dx, oz+(nz-1)*dz, oz]
isnothing(aspect) && (aspect = .5 * nx/nz)

if new_fig
    figure()
end
imshow(scaled', vmin=-a, vmax=a, cmap=cmap, aspect=aspect,
       interpolation=interp, extent=extent)
xlabel("X ($units)")
ylabel("Depth ($units)")
title("$name")

if ~isnothing(save)
    save == True ? filename=name : filename=save
    savefig(filename, bbox_inches="tight", dpi=150)
end
end


text:
using stoosc
rc("font", family="arial")

fig, ax_arr = subplots(2, 2, figsize=(7, 6))

I_max = 0.3

filename = joinpath("arnold_tongues", "model01", "sde_n=100_Ω=1000.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
plot_arnold_tongue(arnold_tongue, I_max, ax=ax_arr[1, 1])
ax_arr[1, 1].set_title("SDE model; n = 100; Ω = 1000")

filename = joinpath("arnold_tongues", "model01", "sde_n=100_Ω=1000_jump.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
plot_arnold_tongue(arnold_tongue, I_max, ax=ax_arr[1, 2])
ax_arr[1, 2].set_title("Gillespie method; n = 100; Ω = 1000")

filename = joinpath("arnold_tongues", "model01", "sde_n=1000_Ω=1000.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
plot_arnold_tongue(arnold_tongue, I_max, ax=ax_arr[2, 1])
ax_arr[2, 1].set_title("SDE model; n = 1000; Ω = 1000")

filename = joinpath("arnold_tongues", "model01", "sde_n=1000_Ω=1000_jump.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
plot_arnold_tongue(arnold_tongue, I_max, ax=ax_arr[2, 2])
ax_arr[2, 2].set_title("Gillespie method; n = 1000; Ω = 1000")

fig.tight_layout()
save_figure(fig, "arnold_tonuges_sde_jump.svg")
# ax_arr[1, 1].tight_layout()
# filename = joinpath("arnold_tongue", "arnold_deterministic.svg")
#
# close("all")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
function local_search!(ip::Vector{Int64}, population::Matrix{Float64},
                        n_population::Int64, n_children::Int64, n_gene::Int64)::Matrix{Float64}
    idx::BitArray{1} = trues(n_population)
    idx[ip[1]] = false

    children::Matrix{Float64} = zeros(n_children, n_gene+1)

    for i in 1:n_children
        ip[2:end] = sample(
            collect(1:n_population)[idx], n_gene+1, replace=false
        )
        children[i,:] = mutation(
            population[ip,:], n_gene
        )
    end

    family::Matrix{Float64} = zeros(n_children+1, n_gene+1)
    @inbounds for i in 1:n_gene+1
        @simd for j in 1:n_children
            family[j,i] = children[j,i]
        end
        family[n_children+1,i] = population[ip[1],i]
    end

    family = sortslices(family, dims=1, by=x->x[end])

    for i in 1:n_gene+1
        @inbounds population[ip[1],i] = family[1,i]  # Best
    end

    population = sortslices(population, dims=1, by=x->x[end])

    return population
end

function mutation(parents::Matrix{Float64}, n_gene::Int64)::Vector{Float64}
    child::Vector{Float64} = NDM(parents, n_gene)
    for i in 1:n_gene
        @inbounds child[i] = clamp(child[i], 0.0, 1.0)
    end

    child[end] = objective(
        child[1:n_gene]
    )

    return child
end

# Normal Distribution Mutation
function NDM(parents::Matrix{Float64},n_gene::Int64)::Vector{Float64}
    child::Vector{Float64} = zeros(n_gene+1)

    GAMMA::Float64 = 0.35/sqrt(n_gene)

    t2::Vector{Float64} = zeros(n_gene)
    centroid::Vector{Float64} = reshape(
        mean(parents[2:end,1:n_gene], dims=1), n_gene
    )

    for i in 1:n_gene+1
        t2 += randn()*GAMMA*(parents[i+1,1:n_gene] - centroid)
    end

    for i in 1:n_gene
        @inbounds child[i] = parents[1,i] + t2[i]
    end

    return child
end


text:
using GBGraph
using Test

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


text:
using PyPlot
include("structures.jl")
include("utilitaires.jl")
include("instance_sur_mesure.jl")

function main()
    instanceSurMesure(false)
end

main()

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
# https://en.wikipedia.org/wiki/MU_puzzle#Solution
#

miu = @theory begin
    # Composition of the string monoid is associative
    x ⋅ (y ⋅ z) => (x ⋅ y) ⋅ z
    # Add a U to the end of any string ending in I
    x ⋅ :I =>  x ⋅ :I ⋅ :U
    # Double the string after the M
    :M ⋅ x  => :M ⋅ x ⋅ x
    # Replace any III with a U
    x ⋅ :I ⋅ :I ⋅ :I ⋅ y => x ⋅ :U ⋅ y
    # Remove any UU
    x ⋅ :U ⋅ :U ⋅ y => x ⋅ y
end

@testset "MU puzzle" begin
    # no matter the timeout we set here,
    # MU is not a theorem of the MIU system
    @test false == areequal(miu, :(M⋅U), :(M⋅I), timeout=12, sizeout=2^13)
end


text:
@testset "REST Client" begin
    url = "http://btc.brane.cc:8332"
    @testset "GET Tx" begin
        key = "f4184fc596403b9d638783cf57adfe4c75c605f6356fbc91338530e9831e9e16"
        @test typeof(get_tx(key)) == Tx
    end
    @testset "GET Headers" begin
        key = "000000000019d6689c085ae165831e934ff763ae46a2a6c172b3f1b60a8ce26f"
        headers = get_headers(key, amount=2)
        @test length(headers) == 2
        @test typeof(headers[1]) == BlockHeader
    end
    @testset "GET BlockHashByHeight" begin
        want = "000000000019d6689c085ae165831e934ff763ae46a2a6c172b3f1b60a8ce26f"
        blockhash = get_blockhashbyheight(0)
        @test blockhash == want
    end
end


text:
module BenchComponent

using BenchmarkTools
import Random
using Overseer
using Parameters

@component @with_kw struct Spatial
    position::NTuple{3, Float64} = (1.0,1.0,1.0)
    velocity::NTuple{3, Float64} = (1.0,1.0,1.0)
end

@component struct Spring
    center::NTuple{3, Float64}
    spring_constant::Float64
end

@component mutable struct Rotation
    omega::Float64
    center::NTuple{3, Float64}
    axis::NTuple{3, Float64}
    function Rotation()
        new(1.0, (2.0,2.0,2.0), (2.0,2.0,2.0))
    end
end

struct Oscillator <: System end

Overseer.requested_components(::Oscillator) = (Spatial, Spring)

function Overseer.update(::Oscillator, m::AbstractLedger)
    spatial = m[Spatial]
    spring = m[Spring]
    # g = group(m, Spatial, Spring)
    @inbounds for e in @entities_in(spatial && spring)
        e_spat  = spatial[e]
        spr     = spring[e]
        v_prev  = e_spat.velocity
        new_v   = v_prev .- (e_spat.position .- spr.center) .* spr.spring_constant
        spatial[e] = Spatial(e_spat.position, new_v)
    end
end

struct Rotator <: System  end
Overseer.requested_components(::Rotator) = (Spatial, Rotation)

function Overseer.update(::Rotator, dio::AbstractLedger)
    rotation  = dio[Rotation]
    spatial   = dio[Spatial]
    dt = 0.01
    @inbounds for e in @entities_in(rotation && spatial)
        e_rotation = rotation[e]
        e_spatial  = spatial[e]
        n          = e_rotation.axis
        r          = e_spatial.position .- e_rotation.center
        theta      = e_rotation.omega * dt
        nnd        = n .* sum(n .* r)
        t = (r[2] * n[3] - r[3] * n[2], r[3]*n[1] - r[1] * n[3] , r[1] * n[2] - r[2] * n[1])
        spatial[e] = Spatial(e_rotation.center .+ nnd .+ (r .- nnd) .* cos(theta) .+ t .* sin(theta), e_spatial.velocity)
    end
end

struct Mover <: System end

Overseer.requested_components(::Mover) = (Spatial, )

function Overseer.update(::Mover, m::AbstractLedger)
    dt = 0.01
    spat = m[Spatial]
    @inbounds for e in @entities_in(spat)
        e_spat = spat[e]
        spat[e] = Spatial(e_spat.position .+ e_spat.velocity.*dt, e_spat.velocity)
    end
end

# Benchmark defs
const suite = BenchmarkGroup()

Random.seed!(1234)

function bench_insertion(c::Component, ids::Vector{Entity}, v)
    for e in ids
        c[e] = v
    end
    return c
end
function bench_access(c::Component, ids::Vector{Entity}, v)
    for e in ids
        v += c[e].position[1]
    end
    return v
end
function bench_access_inbounds(c::Component, ids::Vector{Entity}, v)
    @inbounds for e in ids
        v += c[e].position[1]
    end
    return v
end
function bench_delete(c::Component, ids::Vector{Entity})
    delete!(c, ids)
    return c
end

suite["basic"] = BenchmarkGroup()

const ids = Entity.(unique(rand(1:1000, 1000)))

suite["basic"]["insertion"] = @benchmarkable bench_insertion(c, $ids, $(Spatial((123.0, 1.2, 1.0), (0.4, 12.0, 234.9)))) setup=(c=Component{Spatial}()) evals=1

suite["basic"]["access"] = @benchmarkable bench_access(c, $ids, 0.0) setup=(c = Component{Spatial}(); bench_insertion(c, $ids, Spatial((123.0, 1.2, 1.0), (0.4, 12.0, 234.9))))

const c_full = Component{Spatial}()
bench_insertion(c_full, ids, Spatial((123.0, 1.2, 1.0), (0.4, 12.0, 234.9)))

suite["basic"]["access inbounds"] = @benchmarkable bench_access_inbounds(c, $ids, 0.0) setup=(c = deepcopy(c_full))

suite["basic"]["deletion"] = @benchmarkable bench_delete(c, $ids) setup=(c = deepcopy(c_full)) evals=1

suite["real life"] = BenchmarkGroup()
suite["real life"]["creation"] = BenchmarkGroup()

st = Stage(:simulation, [Oscillator(), Rotator(), Mover()])

suite["real life"]["creation"]["ledger"] = @benchmarkable m = Ledger($st)
function e_fill(m)
    for i=1:1000
        e1 = Entity(m,
                    Spatial((i, 1.0, 1.0), (0.0, 0.0, 0.0)),
                    Spring((1.0, 0.0, 0.0), 0.01))

        e1 = Entity(m,
                    Spatial((i, 1.0, 1.0), (0.0, 0.0, 0.0)),
                    Spring((1.0, 0.0, 0.0), 0.01),
                    Rotation())

    end
end

suite["real life"]["creation"]["filling entities"] = @benchmarkable e_fill(m) setup=(m = Ledger(st))

suite["real life"]["update"] = BenchmarkGroup()

suite["real life"]["update"]["old school empty"] = @benchmarkable update(m) setup=(m = Ledger(st))
suite["real life"]["update"]["old school full"] = @benchmarkable update(m) setup=(m = Ledger(st); e_fill(m))

function update_new(::Oscillator, m::AbstractLedger)
    @inbounds for e in @entities_in(m, Spatial && Spring)
        new_v   = e.velocity .- (e.position .- e[Spring].center) .* e.spring_constant
        e[Spatial] = Spatial(e.position, new_v)
    end
end

function update_new(::Rotator, m::AbstractLedger)
    dt = 0.01
    @inbounds for e in @entities_in(m, Rotation && Spatial)
        n          = e.axis
        r          = e.position .- e.center
        theta      = e.omega * dt
        nnd        = n .* sum(n .* r)
        t = (r[2] * n[3] - r[3] * n[2], r[3]*n[1] - r[1] * n[3] , r[1] * n[2] - r[2] * n[1])
        e[Spatial] = Spatial(e.center .+ nnd .+ (r .- nnd) .* cos(theta) .+ t .* sin(theta), e.velocity)
    end
end

function update_new(::Mover, m::AbstractLedger)
    dt = 0.01
    @inbounds for e in @entities_in(m, Spatial)
        e[Spatial] = Spatial(e.position .+ e.velocity.*dt, e.velocity)
    end
end

suite["real life"]["update"]["new school empty"] = @benchmarkable (update_new(Oscillator(), m);
                                                      update_new(Rotator(), m);
                                                      update_new(Mover(), m)) setup=(m = Ledger(st))
suite["real life"]["update"]["new school full"] = @benchmarkable (update_new(Oscillator(), m);
                                                      update_new(Rotator(), m);
                                                      update_new(Mover(), m)) setup=(m = Ledger(st); e_fill(m))

end
BenchComponent.suite


text:
# -*- coding: utf-8 -*-
module GeneralizedPhaseReduction
    export eye, vec,
           update_z!, find_stable_periodic_solution, phase_sensitivity_func, approx_z,
           generalized_phase_sensitivity_func,
           get_ode_integrator, get_ode_solution,
           compute_QΘ, compute_IsΘ, conventinal_coupled_updateΘ,
           generalized_coupled_updateΘ_I, generalized_coupled_updateΘ_PQ,
           phase2cum_phase, coupled_original_system, coupled_conventinal_phase_model,
           coupled_generalized_phase_model_I, coupled_generalized_phase_model_PQ

    using ForwardDiff
    using LinearAlgebra
    using Random
    using Statistics
    using ProgressMeter
    using Interpolations
    using DifferentialEquations
    using Suppressor
    using Printf

    include("utils.jl")
    include("phase_sensitivity.jl")
    include("generalized_phase_sensitivity.jl")
    include("differential_equations_extensions.jl")
    include("coupled_phase_equation.jl")
end


text:
using Test

include("loocv.jl")
#test problem
A = rand(5, 3)
b = rand(5, 1)
xstar = A\b
r = b-A*xstar
aqr = qr(A)
remainders, minimizers = lsq_loocv(A, aqr, r, xstar)
rr, mm = lsq_loocv(A, aqr, r, xstar, 3)

remainders_actual = similar(remainders)
minimizers_actual = similar(minimizers)
for i = 1:size(A, 1)
    A1 = A[[1:i-1;i+1:end], :]
    b1 = b[[1:i-1;i+1:end]]
    x1 = A1\b1
    minimizers_actual[:, i] = x1
    remainders_actual[:, i] = b-A*x1
end
@testset "lsq_loocv, compute all rem and min" begin
    @test norm(remainders_actual - remainders)<1e-10
    @test norm(minimizers_actual - minimizers)<1e-10
end

#println("remainder error: ", norm(remainders_actual - remainders))
#println("minimizer error: ", norm(minimizers_actual - minimizers))

###
### test lsq_loocv for Delete-one-point cross validation
###
xstar = A\b
r = b-A*xstar
aqr = qr(A)
remainders1 = zeros(5, 5)
minimizers1 = zeros(3, 5)

for i = 1:size(A, 1)
    rr, mm = lsq_loocv(A, aqr, r, xstar, i)
    remainders1[:, i] = rr
    minimizers1[:, i] = mm
end

remainders_actual = similar(remainders1)
minimizers_actual = similar(minimizers1)
for i = 1:size(A, 1)
    A1 = A[[1:i-1;i+1:end], :]
    b1 = b[[1:i-1;i+1:end]]
    x1 = A1\b1
    minimizers_actual[:, i] = x1
    remainders_actual[:, i] = b-A*x1
end
#println("remainder error: ", norm(remainders_actual - remainders1))
#println("minimizer error: ", norm(minimizers_actual - minimizers1))

@testset "lsq_loocv single deletion" begin
    @test norm(remainders_actual - remainders1)<1e-10
    @test norm(minimizers_actual - minimizers1)<1e-10
end

###
### lin_sys_loocv test
###

K = [1.0         0.00369786  0.0608101;
0.00369786  1.0         0.0608101;
0.0608101   0.0608101   1.0]
println("K: ", K)
y = rand(3, 10)
#cholK = incremental_cholesky!(K, 3)

U = cholesky(K).U
c = K\y
#solve subsystem of Kc=y
c1 = lin_sys_loocv(c, U, 1)
c2 = lin_sys_loocv(c, U, 2)
c3 = lin_sys_loocv(c, U, 3)
Cs = zeros(6, 10)
for i = 1:3
    Ki = K[[1:i-1;i+1:end], [1:i-1;i+1:end]]
    Cs[2*i-1:2*i, :] = Ki\y[[1:i-1;i+1:end], :]
end
#println("error of indirect error computation: ", norm(Cs - vcat(vcat(c1, c2), c3)))

cholK = incremental_cholesky!(K, 3)
c1IC = lin_sys_loocv_IC(c, cholK, 1)
c2IC = lin_sys_loocv_IC(c, cholK, 2)
c3IC = lin_sys_loocv_IC(c, cholK, 3)
#println("error of indirect error computation using Incremental_Cholesky (IC): ", norm(Cs - vcat(vcat(c1IC, c2IC), c3IC)))
@testset "lin_sys_loocv_IC" begin
    @test norm(Cs - vcat(vcat(c1IC, c2IC), c3IC))<1e-10
end

###
### Larger test problem. Not activated for now.
###
if true
    n = 50
    K = rand(n, n)
    K = K'*K+UniformScaling(10.0)
    y = rand(n)
    U = cholesky(K).U
    c = K\y
    #solve subsystem of Kc=y
    @time begin
    indirectCs = zeros(n-1, n)
    for i = 1:n
        indirectCs[:, i] = lin_sys_loocv(c, U, i)
    end
    end
    @time begin
    Cs = zeros(n-1, n)
    for i = 1:n
        Ki = K[[1:i-1;i+1:end], [1:i-1;i+1:end]]
        Cs[:, i] = Ki\y[[1:i-1;i+1:end]]
    end
    end
    #println("error of indirect error computation for size $n problem: ", norm(Cs - indirectCs))
    @testset "lin_sys_loocv_IC" begin
    @test norm(Cs - indirectCs) <1e-10
end

    cholK = incremental_cholesky!(K, n)
    indirectCs2 = zeros(n-1, n)
    for i = 1:n
        indirectCs2[:, i] = lin_sys_loocv_IC(c, cholK, i)
    end
    #println("error of indirect error computation for size $n problem using Incremental Cholesky: ", norm(Cs - indirectCs2))

    @testset "lin_sys_loocv_IC" begin
    @test norm(Cs - indirectCs2)<1e-10
end

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    SVD(; kwargs...)

Imputes the missing values in a matrix using an expectation maximization (EM) algorithm
over low-rank SVD approximations.

# Keyword Arguments
* `init::Imputor`: initialization method for missing values (default: Substitute())
* `rank::Union{Int, Nothing}`: rank of the SVD approximation (default: nothing meaning start and 0 and increase)
* `tol::Float64`: convergence tolerance (default: 1e-10)
* `maxiter::Int`: Maximum number of iterations if convergence is not achieved (default: 100)
* `limits::Unoin{Tuple{Float64, Float64}, Nothing}`: Bound the possible approximation values (default: nothing)
* `verbose::Bool`: Whether to display convergence progress (default: true)

# References
* Troyanskaya, Olga, et al. "Missing value estimation methods for DNA microarrays." Bioinformatics 17.6 (2001): 520-525.
"""
struct SVD <: Imputor
    init::Imputor
    rank::Union{Int, Nothing}
    tol::Float64
    maxiter::Int
    limits::Union{Tuple{Float64, Float64}, Nothing}
    verbose::Bool
end

function SVD(;
    init=Substitute(), rank=nothing, tol=1e-10, maxiter=100, limits=nothing, verbose=true
)
    SVD(init, rank, tol, maxiter, limits, verbose)
end

function impute!(data::AbstractMatrix{Union{T, Missing}}, imp::SVD; dims=nothing) where T<:Real
    n, p = size(data)
    k = imp.rank === nothing ? 0 : min(imp.rank, p-1)
    S = zeros(T, min(n, p))
    X = zeros(T, n, p)

    # Get our before and after views of our missing and non-missing data
    mmask = ismissing.(data)
    omask = .!mmask

    mdata = data[mmask]
    mX = X[mmask]
    odata = data[omask]
    oX = X[omask]

    # Fill in the original data
    impute!(data, imp.init; dims=dims)

    C = sum(abs2, mdata - mX) / sum(abs2, mdata)
    err = mean(abs.(odata - oX))
    @debug("Before", Diff=sum(mdata - mX), MAE=err, convergence=C, normsq=sum(abs2, mdata), mX[1])

    for i in 1:imp.maxiter
        if imp.rank === nothing
            k = min(k + 1, p - 1, n - 1)
        end

        # Compute the SVD and produce a low-rank approximation of the data
        F = LinearAlgebra.svd(data)
        # println(join([size(S), size(F.S), size(F.U), size(F.Vt)], ", "))

        S[1:k] .= F.S[1:k]
        X = F.U * Diagonal(S) * F.Vt

        # Clamp the values if necessary
        imp.limits !== nothing && clamp!(X, imp.limits...)

        # Test for convergence
        mdata = data[mmask]
        mX = X[mmask]
        odata = data[omask]
        oX = X[omask]

        # println(join([size(mdata), size(mX)], ", "))
        C = sum(abs2, mdata - mX) / sum(abs2, mdata)

        # Print the error between reconstruction and observed inputs
        if imp.verbose
            err = mean(abs.(odata - oX))
            @debug("Iteration", i, Diff=sum(mdata - mX), MAE=err, MSS=sum(abs2, mdata), convergence=C)
        end

        # Update missing values
        data[mmask] .= X[mmask]

        isfinite(C) && C < imp.tol && break
    end

    return data
end

impute!(data::AbstractMatrix{Missing}, imp::SVD; kwargs...) = data

function impute(data::AbstractMatrix{Union{T, Missing}}, imp::SVD; kwargs...) where T<:Real
    return impute!(trycopy(data), imp; kwargs...)
end

impute(data::AbstractMatrix{Missing}, imp::SVD; kwargs...) = trycopy(data)


text:
include("./src/trieMusic.jl")
include("./src/midi.jl")

############## midiデータからtrie木を作る ##############
pitchSequence, positionSequence, velocitySequence, = midiModule.getMidiDatas("./midi/bach_2.mid")
#階差数列に変換
differenceSequence = trieMusic.seqenceToDifferenceSequence(pitchSequence)
#2値の数列に変換
binarySequence = trieMusic.sequenceToBinarySequence(differenceSequence)
#LOUDSのTrie木に変換し、トライのビット配列、ラベル（[そのノードのラベル、そのラベルが出現したシーケンス番号１つ目、２つ目、...]という形式）、
#トライ木の各階層のノードの数、ビット配列の階層を示す区切り位置を取得する
bitAry, label, eachLevelNodesAccumulatedNumber, levelBoundaryBaIdx = @time trieMusic.sequenceToLouds(binarySequence)
#全部分列の情報を取得する
subSequences, startSeqIdxes, endSeqIdxes = trieMusic.getAllSubSequences(label, pitchSequence, eachLevelNodesAccumulatedNumber)
#全部分列のDTW距離を計算する
println(trieMusic.calcDtwDistances(subSequences))
#ルートノードの直下から指定したノードまでの部分列の形状を持つ、実際の部分列群を取得
# subSequences, startSeqIdxes, endSeqIdxes = trieMusic.getMatchSubsequences(pitchSequence, label, eachLevelNodesAccumulatedNumber, 421)
# #取得した部分列群を表示
# trieMusic.plotSubsequences(pitchSequence, subSequences, startSeqIdxes, endSeqIdxes)
#PlantUMLファイルを生成　長さによっては出力ファイルが大きくなるので注意
trieMusic.trieToPuml(bitAry, label)

include("./src/musicXML.jl")
musicXML.generateMusicXML()


text:
import DielectricUpdateTechique as DUTCH
using Test

x = collect(0:1f-3:2f-3)
y = collect(0:1f-3:2f-3)
z = collect(0:1f-3:2f-3)

vσ = rand(Float32,3, 3, 3);
vϵ = rand(Float32,3, 3, 3);

bg = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

vσ[1:2,2,2] = 2 .* vσ[1:2,2,2];

imp = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

Δ, loc = DUTCH.get_update_maps(bg,imp)

@test round(Δ[:σˣ][2,2,2], digits=3) == round(imp[:σˣ][2,2,2] - bg[:σˣ][2,2,2], digits=3)
@test round(Δ[:σˣ][1,2,2], digits=3) == round(imp[:σˣ][1,2,2] - bg[:σˣ][1,2,2], digits=3)

@test round(Δ[:σʸ][2,2,2], digits=3) == round(imp[:σʸ][2,2,2] - bg[:σʸ][2,2,2], digits=3)
@test round(Δ[:σʸ][1,2,2], digits=3) == round(imp[:σʸ][1,2,2] - bg[:σʸ][1,2,2], digits=3)

@test round(Δ[:σᶻ][2,2,2], digits=3) == round(imp[:σᶻ][2,2,2] - bg[:σᶻ][2,2,2], digits=3)
@test round(Δ[:σᶻ][1,2,2], digits=3) == round(imp[:σᶻ][1,2,2] - bg[:σᶻ][1,2,2], digits=3)

# there are 20 edges to be updated for 2 neighbouring cells
@test sum(vec(loc[:updateLocˣ])) + sum(vec(loc[:updateLocʸ])) + sum(vec(loc[:updateLocᶻ])) == 20

vσ = rand(Float32,3, 3, 3);
vϵ = rand(Float32,3, 3, 3);

bg = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

vσ[1,2,2] = 2 .* vσ[1,2,2];
vσ[3,2,2] = 2 .* vσ[3,2,2];

imp = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

Δ, loc = DUTCH.get_update_maps(bg,imp)
# there are 24 edges to be updated for 2 non neighbouring cells
@test sum(vec(loc[:updateLocˣ])) + sum(vec(loc[:updateLocʸ])) + sum(vec(loc[:updateLocᶻ])) == 24

# test for get_S_matrix(_gpu) and get_C_matrix

x = collect(0:1f-3:5f-3)
y = collect(0:1f-3:5f-3)
z = collect(0:1f-3:5f-3)

vσ = rand(Float32,5, 5, 5);
vϵ = rand(Float32,5, 5, 5);

bg = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

vσ[2,2,2] = 2 .* vσ[2,2,2];
vσ[4,2,2] = 2 .* vσ[4,2,2];

imp = DUTCH.cell_to_yee_dielectric(vσ, vϵ, x, y, z)

Δ, loc = DUTCH.get_update_maps(bg,imp)

S = DUTCH.get_S_matrix(loc)

@test count(S .== 1) == 24
@test size(S) == ((length(x)-1)*(length(y)-2)*(length(z)-2) + (length(x)-2)*(length(y)-1)*(length(z)-2) + (length(x)-2)*(length(y)-2)*(length(z)-1), 24)

#=
S = DUTCH.get_S_matrix_gpu(loc) |> collect

@test count(S .== 1) == 24
@test size(S) == ((length(x)-1)*(length(y)-2)*(length(z)-2) + (length(x)-2)*(length(y)-1)*(length(z)-2) + (length(x)-2)*(length(y)-2)*(length(z)-1), 24)
=#

@test sum(DUTCH._get_number_of_updated_edges(loc)) == 24

f = 300e6; #300MHz
m = DUTCH.get_constants(f);
C = DUTCH.get_C_matrix(Δ,m)

@test length(C) == 24
@test sum(imag.(C)) == 0
@test round(sum(real.(C)),digits=1)  == round(sum(real.(Δ.σˣ)) + sum(real.(Δ.σʸ)) + sum(real.(Δ.σᶻ)),digits = 1)
@test typeof(C) == Vector{ComplexF32}

χ = DUTCH._set_dielectric(bg,m,[3,3,3]);
@test abs.(χ.vac) == 1
@test size(χ.patient) == (16,16,16,3)
@test abs.(χ.patient[8,7,7,1]) == 0
@test abs.(χ.patient[3,3,3,1]) != 0
@test abs.(χ.patient[7,8,7,2]) == 0
@test abs.(χ.patient[3,3,3,2]) != 0
@test abs.(χ.patient[7,7,8,3]) == 0
@test abs.(χ.patient[3,3,3,3]) != 0


text:
typealias Steps Vector{Int}

immutable Chord <: Sequencable
  steps :: Steps
  scale_size :: Int
  name :: Nullable{Symbol}
  seq :: Sequencable
end

Chord(s) = Chord(s, 7)
Chord(s::Steps, size::Int) = Chord(s, size, nothing)
Chord(steps, size, name) = Chord(steps, size, name, Event())

function Chord(s::Symbol, size::Int)
  if s == :triad && size == 7
    Chord([0, 2, 4], 7, s)
  end

  # TODO this is really slow to compile...
  # @match (s, size) begin
  #   (:triad, 7) => Chord([0, 2, 4], 7, s)
  #   (:sixth, 7) => Chord([0, 2, 4, 5], 7, s)
  #   (:seventh, 7) => Chord([0, 2, 4, 6], 7, s)
  #   (:ninth, 7) => Chord([0, 2, 4, 8], 7, s)
  #   (:eleventh, 7) => Chord([0, 2, 4, 8, 10], 7, s)
  # end
end

@delegate Chord.steps [
  Base.get, Base.getindex, Base.haskey,
  Base.getkey, Base.start, Base.done, Base.next,
  Base.isempty, Base.length, Base.hash
]

@delegate Chord.seq [
  Base.rot180, Base.rotr90, Base.rotl90
]

Base.convert(::Type{Vector{Sequencable}}, c::Chord) =
  Vector{Sequencable}(map(step -> merge(c.seq, :degree => step // c.scale_size), c.steps))

Base.convert(::Type{HSeq}, c::Chord) = HSeq(Vector{Sequencable}(c))
Base.convert(::Type{VSeq}, c::Chord) = VSeq(Vector{Sequencable}(c))
Base.merge(c::Chord, s::Sequencable) = Chord(c.steps, c.scale_size, c.name, merge(c.seq, s))
Base.merge(s::Sequencable, c::Chord) = Chord(c.steps, c.scale_size, c.name, merge(s, c.seq))

function circshift(c::Chord, count)
  # TODO wrap around the scale_size
  c
end


text:
NaiveNASlib.nin(l) = nin(layertype(l), l)
NaiveNASlib.nout(l) = nout(layertype(l), l)

NaiveNASlib.nin(t::FluxLayer, l) = throw(ArgumentError("Not implemented for $t"))
NaiveNASlib.nout(t::FluxLayer, l) = throw(ArgumentError("Not implemented for $t"))

NaiveNASlib.nin(::FluxParLayer, l) = size(weights(l), indim(l))
NaiveNASlib.nout(::FluxParLayer, l) = size(weights(l), outdim(l))
NaiveNASlib.nout(::FluxDepthwiseConv, l) = size(weights(l), outdim(l)) * nin(l)

NaiveNASlib.nin(::FluxParInvLayer, l) = nout(l)

NaiveNASlib.nout(::FluxDiagonal, l) = length(weights(l))
NaiveNASlib.nout(::FluxParInvLayer, l::LayerNorm) = nout(l.diag)
NaiveNASlib.nout(::FluxParNorm, l) = length(l.β)

NaiveNASlib.nout(::FluxRecurrent, l) = div(size(weights(l), outdim(l)), outscale(l))

NaiveNASlib.minΔninfactor(::FluxLayer, l) = 1
NaiveNASlib.minΔnoutfactor(::FluxLayer, l) = 1

outscale(l) = outscale(layertype(l))
outscale(::FluxRnn) = 1
outscale(::FluxLstm) = 4
outscale(::FluxGru) = 3

indim(l) = indim(layertype(l))
outdim(l) = outdim(layertype(l))
actdim(l) = actdim(layertype(l))
actrank(l) = actrank(layertype(l))

indim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
outdim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
actdim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
actrank(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))

indim(::FluxDense) = 2
outdim(::FluxDense) = 1
actdim(::FluxDense) = 1
actrank(::FluxDense) = 1

indim(::FluxDiagonal) = 1
outdim(::FluxDiagonal) = 1
actdim(::FluxDiagonal) = 1
actrank(::FluxDiagonal) = 1

indim(::FluxRecurrent) = 2
outdim(::FluxRecurrent) = 1
actdim(::FluxRecurrent) = 1
actrank(::FluxRecurrent) = 2

indim(::FluxConvolutional{N}) where N = 2+N
outdim(::FluxConvolutional{N}) where N = 1+N
actdim(::FluxConvolutional{N}) where N = 1+N
actrank(::FluxConvolutional{N}) where N = 1+N
indim(::Union{FluxConv{N}, FluxCrossCor{N}}) where N = 1+N
outdim(::Union{FluxConv{N}, FluxCrossCor{N}}) where N = 2+N
# Note: Contrary to other ML frameworks, bias seems to always be present in Flux
weights(l) = weights(layertype(l), l)
bias(l) = bias(layertype(l), l)

weights(::FluxDense, l) = l.W.data
bias(::FluxDense, l) = l.b.data

weights(::FluxConvolutional, l) = l.weight.data
bias(::FluxConvolutional, l) = l.bias.data

weights(::FluxDiagonal, l) = l.α.data
bias(::FluxDiagonal, l) = l.β.data

weights(::FluxRecurrent, l) = l.cell.Wi.data
bias(::FluxRecurrent, l) = l.cell.b.data

hiddenweights(l) = hiddenweights(layertype(l), l)
hiddenweights(::FluxRecurrent, l) = l.cell.Wh.data
hiddenstate(l) = hiddenstate(layertype(l), l)
hiddenstate(::FluxRecurrent, l) = Flux.hidden(l.cell).data
hiddenstate(::FluxLstm, l) = [h.data for h in Flux.hidden(l.cell)]
state(l) = state(layertype(l), l)
state(::FluxRecurrent, l) = l.state.data
state(::FluxLstm, l) = [h.data for h in l.state]

function NaiveNASlib.compconstraint!(s, ::FluxLayer, data) end

function NaiveNASlib.compconstraint!(s::NaiveNASlib.AbstractJuMPΔSizeStrategy, ::FluxDepthwiseConv, data)
  # Add constraint that nout(l) == n * nin(l) where n is integer
  fv_out = @variable(data.model, integer=true)
  ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))

  # Would have preferred to have "data.noutdict[data.vertex] == data.noutdict[ins[i]] * fv_out", but it is not linear
  @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] == data.noutdict[ins[i]] + nin(data.vertex)[] * fv_out)
  @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] >= data.noutdict[ins[i]])

  # Inputs which does not have a variable, possibly because all_in_Δsize_graph did not consider it to be part of the set of vertices which may change
  fixedins = filter(vin -> vin ∉ ins, inputs(data.vertex))
  @constraint(data.model, [i=1:length(fixedins)], data.noutdict[data.vertex] == nout(fixedins[i]) + nin(data.vertex)[] * fv_out)
end
# compconstraint! for AbstractJuMPSelectionStrategy not needed as there currently is no strategy which allows size changes