text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: groundstation

using AWS.Compat
using AWS.UUIDs
"""
    CancelContact()

Cancels a contact with a specified contact ID.

# Required Parameters
- `contactId`: UUID of a contact.

"""

cancel_contact(contactId; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/contact/$(contactId)"; aws_config=aws_config)
cancel_contact(contactId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/contact/$(contactId)", args; aws_config=aws_config)

"""
    CreateConfig()

Creates a Config with the specified configData parameters. Only one type of configData can be specified.

# Required Parameters
- `configData`: Parameters of a Config.
- `name`: Name of a Config.

# Optional Parameters
- `tags`: Tags assigned to a Config.
"""

create_config(configData, name; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/config", Dict{String, Any}("configData"=>configData, "name"=>name); aws_config=aws_config)
create_config(configData, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/config", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("configData"=>configData, "name"=>name), args)); aws_config=aws_config)

"""
    CreateDataflowEndpointGroup()

Creates a DataflowEndpoint group containing the specified list of DataflowEndpoint objects. The name field in each endpoint is used in your mission profile DataflowEndpointConfig to specify which endpoints to use during a contact. When a contact uses multiple DataflowEndpointConfig objects, each Config must match a DataflowEndpoint in the same group.

# Required Parameters
- `endpointDetails`: Endpoint details of each endpoint in the dataflow endpoint group.

# Optional Parameters
- `tags`: Tags of a dataflow endpoint group.
"""

create_dataflow_endpoint_group(endpointDetails; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/dataflowEndpointGroup", Dict{String, Any}("endpointDetails"=>endpointDetails); aws_config=aws_config)
create_dataflow_endpoint_group(endpointDetails, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/dataflowEndpointGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endpointDetails"=>endpointDetails), args)); aws_config=aws_config)

"""
    CreateMissionProfile()

Creates a mission profile.  dataflowEdges is a list of lists of strings. Each lower level list of strings has two elements: a from ARN and a to ARN.

# Required Parameters
- `dataflowEdges`: A list of lists of ARNs. Each list of ARNs is an edge, with a from Config and a to Config.
- `minimumViableContactDurationSeconds`: Smallest amount of time in seconds that you’d like to see for an available contact. AWS Ground Station will not present you with contacts shorter than this duration.
- `name`: Name of a mission profile.
- `trackingConfigArn`: ARN of a tracking Config.

# Optional Parameters
- `contactPostPassDurationSeconds`: Amount of time after a contact ends that you’d like to receive a CloudWatch event indicating the pass has finished.
- `contactPrePassDurationSeconds`: Amount of time prior to contact start you’d like to receive a CloudWatch event indicating an upcoming pass.
- `tags`: Tags assigned to a mission profile.
"""

create_mission_profile(dataflowEdges, minimumViableContactDurationSeconds, name, trackingConfigArn; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/missionprofile", Dict{String, Any}("dataflowEdges"=>dataflowEdges, "minimumViableContactDurationSeconds"=>minimumViableContactDurationSeconds, "name"=>name, "trackingConfigArn"=>trackingConfigArn); aws_config=aws_config)
create_mission_profile(dataflowEdges, minimumViableContactDurationSeconds, name, trackingConfigArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/missionprofile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("dataflowEdges"=>dataflowEdges, "minimumViableContactDurationSeconds"=>minimumViableContactDurationSeconds, "name"=>name, "trackingConfigArn"=>trackingConfigArn), args)); aws_config=aws_config)

"""
    DeleteConfig()

Deletes a Config.

# Required Parameters
- `configId`: UUID of a Config.
- `configType`: Type of a Config.

"""

delete_config(configId, configType; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/config/$(configType)/$(configId)"; aws_config=aws_config)
delete_config(configId, configType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/config/$(configType)/$(configId)", args; aws_config=aws_config)

"""
    DeleteDataflowEndpointGroup()

Deletes a dataflow endpoint group.

# Required Parameters
- `dataflowEndpointGroupId`: UUID of a dataflow endpoint group.

"""

delete_dataflow_endpoint_group(dataflowEndpointGroupId; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)"; aws_config=aws_config)
delete_dataflow_endpoint_group(dataflowEndpointGroupId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)", args; aws_config=aws_config)

"""
    DeleteMissionProfile()

Deletes a mission profile.

# Required Parameters
- `missionProfileId`: UUID of a mission profile.

"""

delete_mission_profile(missionProfileId; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
delete_mission_profile(missionProfileId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/missionprofile/$(missionProfileId)", args; aws_config=aws_config)

"""
    DescribeContact()

Describes an existing contact.

# Required Parameters
- `contactId`: UUID of a contact.

"""

describe_contact(contactId; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/contact/$(contactId)"; aws_config=aws_config)
describe_contact(contactId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/contact/$(contactId)", args; aws_config=aws_config)

"""
    GetConfig()

Returns Config information. Only one Config response can be returned.

# Required Parameters
- `configId`: UUID of a Config.
- `configType`: Type of a Config.

"""

get_config(configId, configType; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/config/$(configType)/$(configId)"; aws_config=aws_config)
get_config(configId, configType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/config/$(configType)/$(configId)", args; aws_config=aws_config)

"""
    GetDataflowEndpointGroup()

Returns the dataflow endpoint group.

# Required Parameters
- `dataflowEndpointGroupId`: UUID of a dataflow endpoint group.

"""

get_dataflow_endpoint_group(dataflowEndpointGroupId; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)"; aws_config=aws_config)
get_dataflow_endpoint_group(dataflowEndpointGroupId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)", args; aws_config=aws_config)

"""
    GetMinuteUsage()

Returns the number of minutes used by account.

# Required Parameters
- `month`: The month being requested, with a value of 1-12.
- `year`: The year being requested, in the format of YYYY.

"""

get_minute_usage(month, year; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/minute-usage", Dict{String, Any}("month"=>month, "year"=>year); aws_config=aws_config)
get_minute_usage(month, year, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/minute-usage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("month"=>month, "year"=>year), args)); aws_config=aws_config)

"""
    GetMissionProfile()

Returns a mission profile.

# Required Parameters
- `missionProfileId`: UUID of a mission profile.

"""

get_mission_profile(missionProfileId; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
get_mission_profile(missionProfileId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile/$(missionProfileId)", args; aws_config=aws_config)

"""
    GetSatellite()

Returns a satellite.

# Required Parameters
- `satelliteId`: UUID of a satellite.

"""

get_satellite(satelliteId; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/satellite/$(satelliteId)"; aws_config=aws_config)
get_satellite(satelliteId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/satellite/$(satelliteId)", args; aws_config=aws_config)

"""
    ListConfigs()

Returns a list of Config objects.

# Optional Parameters
- `maxResults`: Maximum number of Configs returned.
- `nextToken`: Next token returned in the request of a previous ListConfigs call. Used to get the next page of results.
"""

list_configs(; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/config"; aws_config=aws_config)
list_configs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/config", args; aws_config=aws_config)

"""
    ListContacts()

Returns a list of contacts. If statusList contains AVAILABLE, the request must include groundStation, missionprofileArn, and satelliteArn.

# Required Parameters
- `endTime`: End time of a contact.
- `startTime`: Start time of a contact.
- `statusList`: Status of a contact reservation.

# Optional Parameters
- `groundStation`: Name of a ground station.
- `maxResults`: Maximum number of contacts returned.
- `missionProfileArn`: ARN of a mission profile.
- `nextToken`: Next token returned in the request of a previous ListContacts call. Used to get the next page of results.
- `satelliteArn`: ARN of a satellite.
"""

list_contacts(endTime, startTime, statusList; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/contacts", Dict{String, Any}("endTime"=>endTime, "startTime"=>startTime, "statusList"=>statusList); aws_config=aws_config)
list_contacts(endTime, startTime, statusList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/contacts", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "startTime"=>startTime, "statusList"=>statusList), args)); aws_config=aws_config)

"""
    ListDataflowEndpointGroups()

Returns a list of DataflowEndpoint groups.

# Optional Parameters
- `maxResults`: Maximum number of dataflow endpoint groups returned.
- `nextToken`: Next token returned in the request of a previous ListDataflowEndpointGroups call. Used to get the next page of results.
"""

list_dataflow_endpoint_groups(; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup"; aws_config=aws_config)
list_dataflow_endpoint_groups(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup", args; aws_config=aws_config)

"""
    ListGroundStations()

Returns a list of ground stations.

# Optional Parameters
- `maxResults`: Maximum number of ground stations returned.
- `nextToken`: Next token that can be supplied in the next call to get the next page of ground stations.
- `satelliteId`: Satellite ID to retrieve on-boarded ground stations.
"""

list_ground_stations(; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/groundstation"; aws_config=aws_config)
list_ground_stations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/groundstation", args; aws_config=aws_config)

"""
    ListMissionProfiles()

Returns a list of mission profiles.

# Optional Parameters
- `maxResults`: Maximum number of mission profiles returned.
- `nextToken`: Next token returned in the request of a previous ListMissionProfiles call. Used to get the next page of results.
"""

list_mission_profiles(; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile"; aws_config=aws_config)
list_mission_profiles(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile", args; aws_config=aws_config)

"""
    ListSatellites()

Returns a list of satellites.

# Optional Parameters
- `maxResults`: Maximum number of satellites returned.
- `nextToken`: Next token that can be supplied in the next call to get the next page of satellites.
"""

list_satellites(; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/satellite"; aws_config=aws_config)
list_satellites(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/satellite", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of tags for a specified resource.

# Required Parameters
- `resourceArn`: ARN of a resource.

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    ReserveContact()

Reserves a contact using specified parameters.

# Required Parameters
- `endTime`: End time of a contact.
- `groundStation`: Name of a ground station.
- `missionProfileArn`: ARN of a mission profile.
- `satelliteArn`: ARN of a satellite
- `startTime`: Start time of a contact.

# Optional Parameters
- `tags`: Tags assigned to a contact.
"""

reserve_contact(endTime, groundStation, missionProfileArn, satelliteArn, startTime; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/contact", Dict{String, Any}("endTime"=>endTime, "groundStation"=>groundStation, "missionProfileArn"=>missionProfileArn, "satelliteArn"=>satelliteArn, "startTime"=>startTime); aws_config=aws_config)
reserve_contact(endTime, groundStation, missionProfileArn, satelliteArn, startTime, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/contact", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "groundStation"=>groundStation, "missionProfileArn"=>missionProfileArn, "satelliteArn"=>satelliteArn, "startTime"=>startTime), args)); aws_config=aws_config)

"""
    TagResource()

Assigns a tag to a resource.

# Required Parameters
- `resourceArn`: ARN of a resource tag.
- `tags`: Tags assigned to a resource.

"""

tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Deassigns a resource tag.

# Required Parameters
- `resourceArn`: ARN of a resource.
- `tagKeys`: Keys of a resource tag.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateConfig()

Updates the Config used when scheduling contacts. Updating a Config will not update the execution parameters for existing future contacts scheduled with this Config.

# Required Parameters
- `configData`: Parameters of a Config.
- `configId`: UUID of a Config.
- `configType`: Type of a Config.
- `name`: Name of a Config.

"""

update_config(configData, configId, configType, name; aws_config::AWSConfig=global_aws_config()) = groundstation("PUT", "/config/$(configType)/$(configId)", Dict{String, Any}("configData"=>configData, "name"=>name); aws_config=aws_config)
update_config(configData, configId, configType, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("PUT", "/config/$(configType)/$(configId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("configData"=>configData, "name"=>name), args)); aws_config=aws_config)

"""
    UpdateMissionProfile()

Updates a mission profile. Updating a mission profile will not update the execution parameters for existing future contacts.

# Required Parameters
- `missionProfileId`: UUID of a mission profile.

# Optional Parameters
- `contactPostPassDurationSeconds`: Amount of time after a contact ends that you’d like to receive a CloudWatch event indicating the pass has finished.
- `contactPrePassDurationSeconds`: Amount of time after a contact ends that you’d like to receive a CloudWatch event indicating the pass has finished.
- `dataflowEdges`: A list of lists of ARNs. Each list of ARNs is an edge, with a from Config and a to Config.
- `minimumViableContactDurationSeconds`: Smallest amount of time in seconds that you’d like to see for an available contact. AWS Ground Station will not present you with contacts shorter than this duration.
- `name`: Name of a mission profile.
- `trackingConfigArn`: ARN of a tracking Config.
"""

update_mission_profile(missionProfileId; aws_config::AWSConfig=global_aws_config()) = groundstation("PUT", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
update_mission_profile(missionProfileId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = groundstation("PUT", "/missionprofile/$(missionProfileId)", args; aws_config=aws_config)


text:
#
# User facing struct that unifies a ModelDef and a ModelInstance and delegates
# function calls to one or the other as appropriate.
#
using MacroTools

"""
    modeldef(m)

Return the `ModelDef` contained by Model `m`.
"""
modeldef(m::Model) = m.md

modelinstance(m::Model) = m.mi
modelinstance_def(m::Model) = modeldef(modelinstance(m))

is_built(m::Model) = !(dirty(m.md) || modelinstance(m) === nothing)

is_built(mm::MarginalModel) = (is_built(mm.base) && is_built(mm.marginal))

@delegate compinstance(m::Model, name::Symbol) => mi
@delegate has_comp(m::Model, name::Symbol) => md

@delegate number_type(m::Model) => md

@delegate internal_param_conns(m::Model) => md
@delegate external_param_conns(m::Model) => md

@delegate external_params(m::Model) => md
@delegate external_param(m::Model, name::Symbol; missing_ok=false) => md

@delegate connected_params(m::Model) => md
@delegate unconnected_params(m::Model) => md

@delegate add_connector_comps!(m::Model) => md

"""
    connect_param!(m::Model, dst_comp_path::ComponentPath, dst_par_name::Symbol, src_comp_path::ComponentPath,
        src_var_name::Symbol, backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false, offset::Int=0)

Bind the parameter `dst_par_name` of one component `dst_comp_path` of model `md`
to a variable `src_var_name` in another component `src_comp_path` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `offset` argument indicates the offset
between the destination and the source ie. the value would be `1` if the destination
component parameter should only be calculated for the second timestep and beyond.
"""
@delegate connect_param!(m::Model,
                         dst_comp_path::ComponentPath, dst_par_name::Symbol,
                         src_comp_path::ComponentPath, src_var_name::Symbol,
                         backup::Union{Nothing, Array}=nothing;
                         ignoreunits::Bool=false, offset::Int=0) => md

@delegate connect_param!(m::Model,
                         dst_comp_name::Symbol, dst_par_name::Symbol,
                         src_comp_name::Symbol, src_var_name::Symbol,
                         backup::Union{Nothing, Array}=nothing;
                         ignoreunits::Bool=false, offset::Int=0) => md

@delegate connect_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol) => md

@delegate connect_param!(m::Model, dst::AbstractString, src::AbstractString, backup::Union{Nothing, Array}=nothing;
                         ignoreunits::Bool=false, offset::Int=0) => md

"""
    connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol}, backup::Array; ignoreunits::Bool=false)

Bind the parameter `dst[2]` of one component `dst[1]` of model `md`
to a variable `src[2]` in another component `src[1]` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `offset` argument indicates the offset
between the destination and the source ie. the value would be `1` if the destination
component parameter should only be calculated for the second timestep and beyond.

"""
function connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
                           backup::Union{Nothing, Array}=nothing;
                           ignoreunits::Bool=false, offset::Int=0)
    connect_param!(m.md, dst[1], dst[2], src[1], src[2], backup; ignoreunits=ignoreunits, offset=offset)
end

@delegate disconnect_param!(m::Model, comp_path::ComponentPath, param_name::Symbol) => md
@delegate disconnect_param!(m::Model, comp_name::Symbol, param_name::Symbol) => md

@delegate set_external_param!(m::Model, name::Symbol, value::ModelParameter) => md

@delegate set_external_param!(m::Model, name::Symbol, value::Number;
                              param_dims::Union{Nothing,Array{Symbol}} = nothing) => md

@delegate set_external_param!(m::Model, name::Symbol, value::Union{AbstractArray, AbstractRange, Tuple};
                              param_dims::Union{Nothing,Array{Symbol}} = nothing) => md

@delegate add_internal_param_conn!(m::Model, conn::InternalParameterConnection) => md

# @delegate doesn't handle the 'where T' currently. This is the only instance of it for now...
function set_leftover_params!(m::Model, parameters::Dict{T, Any}) where T
    set_leftover_params!(m.md, parameters)
end

"""
    update_param!(m::Model, name::Symbol, value; update_timesteps = false)

Update the `value` of an external model parameter in model `m`, referenced by
`name`. Optional boolean argument `update_timesteps` with default value `false`
indicates whether to update the time keys associated with the parameter values
to match the model's time index.
"""
@delegate update_param!(m::Model, name::Symbol, value; update_timesteps = false) => md

"""
    update_params!(m::Model, parameters::Dict{T, Any}; update_timesteps = false) where T

For each (k, v) in the provided `parameters` dictionary, `update_param!``
is called to update the external parameter by name k to value v, with optional
Boolean argument update_timesteps. Each key k must be a symbol or convert to a
symbol matching the name of an external parameter that already exists in the
model definition.
"""
@delegate update_params!(m::Model, parameters::Dict; update_timesteps = false) => md

"""
    add_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component indicated by `comp_id` to the model indicated by `m`. The component is added
at the end of the list unless one of the keywords `before` or `after` is specified. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`.
"""
function add_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    comp_def = add_comp!(m.md, comp_id, comp_name; kwargs...)
    return ComponentReference(m.md, comp_name)
end

"""
    add_comp!(
        m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component `comp_def` to the model indicated by `m`. The component is added at
the end of the list unless one of the keywords, `first`, `last`, `before`, `after`. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`.
"""
function add_comp!(m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    return add_comp!(m, comp_def.comp_id, comp_name; kwargs...)
end

"""
    replace_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

Replace the component with name `comp_name` in model `m` with the component
`comp_id` using the same name.  The component is added in the same position as
the old component, unless one of the keywords `before` or `after` is specified.
The component is added with the same first and last values, unless the keywords
`first` or `last` are specified. Optional boolean argument `reconnect` with
default value `true` indicates whether the existing parameter connections
should be maintained in the new component.
"""
function replace_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    comp_def = replace_comp!(m.md, comp_id, comp_name; kwargs...)
    return ComponentReference(m.md, comp_name)
end

function replace_comp!(m::Model, comp_def::ComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    return replace_comp!(m, comp_def.comp_id, comp_name; kwargs...)
end

@delegate ComponentReference(m::Model, name::Symbol) => md

"""
    components(m::Model)

Return an iterator on the components in a model's model instance.
"""
@delegate components(m::Model) => mi

@delegate compdefs(m::Model) => md

@delegate compdef(m::Model, comp_name::Symbol) => md

@delegate Base.length(m::Model) => md

@delegate first_and_step(m::Model) => md

@delegate time_labels(m::Model) => md

# Return the number of timesteps a given component in a model will run for.
@delegate getspan(m::Model, comp_name::Symbol) => md

"""
    datumdef(comp_def::ComponentDef, item::Symbol)

Return a DatumDef for `item` in the given component `comp_def`.
"""
function datumdef(comp_def::ComponentDef, item::Symbol)
    if has_variable(comp_def, item)
        return variable(comp_def, item)

    elseif has_parameter(comp_def, item)
        return parameter(comp_def, item)
    else
        error("Cannot access data item; :$item is not a variable or a parameter in component $(comp_def.comp_id).")
    end
end

datumdef(m::Model, comp_name::Symbol, item::Symbol) = datumdef(compdef(m.md, comp_name), item)

"""
    dim_names(m::Model, comp_name::Symbol, datum_name::Symbol)

Return the dimension names for the variable or parameter `datum_name`
in the given component `comp_name` in model `m`.
"""
dim_names(m::Model, comp_name::Symbol, datum_name::Symbol) = dim_names(compdef(m, comp_name), datum_name)
dim_names(mm::MarginalModel, comp_name::Symbol, datum_name::Symbol) = dim_names(mm.base, comp_name, datum_name)

@delegate dimension(m::Model, dim_name::Symbol) => md

# Allow access of the form my_model[:grosseconomy, :tfp]
@delegate Base.getindex(m::Model, comp_name::Symbol, datum_name::Symbol) => mi

"""
    dim_count(m::Model, dim_name::Symbol)

Return the size of index `dim_name` in model `m`.
"""
@delegate dim_count(m::Model, dim_name::Symbol) => md
@delegate dim_counts(m::Model, dims::Vector{Symbol}) => md
@delegate dim_count_dict(m::Model) => md

"""
    dim_keys(m::Model, dim_name::Symbol)

Return keys for dimension `dim-name` in model `m`.
"""
@delegate dim_keys(m::Model, dim_name::Symbol) => md
"""
    dim_key_dict(m::Model)

Return a dict of dimension keys for all dimensions in model `m`.
"""
@delegate dim_key_dict(m::Model) => md
"""
    dim_values(m::Model, name::Symbol)

Return values for dimension `name` in Model `m`.
"""
@delegate dim_values(m::Model, name::Symbol) => md
"""
    dim_value_dict(m::Model)

Return a dictionary of the values of all dimensions in Model `m`.
"""
@delegate dim_value_dict(m::Model) => md
"""
    set_dimension!(m::Model, name::Symbol, keys::Union{Vector, Tuple, AbstractRange})

Set the values of `m` dimension `name` to integers 1 through `count`, if `keys`` is
an integer; or to the values in the vector or range if `keys`` is either of those types.
"""
@delegate set_dimension!(m::Model, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange}) => md

@delegate check_parameter_dimensions(m::Model, value::AbstractArray, dims::Vector, name::Symbol) => md

@delegate parameter_names(m::Model, comp_name::Symbol) => md

@delegate parameter_dimensions(m::Model, comp_name::Symbol, param_name::Symbol) => md

@delegate parameter_unit(m::Model, comp_name::Symbol, param_name::Symbol) => md

parameter(m::Model, comp_def::ComponentDef, param_name::Symbol) = parameter(comp_def, param_name)

parameter(m::Model, comp_name::Symbol, param_name::Symbol) = parameter(m, compdef(m, comp_name), param_name)

"""
    parameters(m::Model, comp_name::Symbol)

Return a list of the parameter definitions for `comp_name` in model `m`.
"""
parameters(m::Model, comp_name::Symbol) = parameters(compdef(m, comp_name))

variable(m::Model, comp_name::Symbol, var_name::Symbol) = variable(compdef(m, comp_name), var_name)

@delegate variable_unit(m::Model, comp_path::ComponentPath, var_name::Symbol) => md

@delegate variable_dimensions(m::Model, comp_path::ComponentPath, var_name::Symbol) => md

"""
    variables(m::Model, comp_name::Symbol)

Return an iterator on the variable definitions for `comp_name` in model `m`.
"""
variables(m::Model, comp_name::Symbol) = variables(compdef(m, comp_name))

@delegate variable_names(m::Model, comp_name::Symbol) => md

"""
    set_external_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims)

Add a one or two dimensional (optionally, time-indexed) array parameter `name`
with value `value` to the model `m`.
"""
@delegate set_external_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims) => md

"""
    set_external_scalar_param!(m::Model, name::Symbol, value::Any)

Add a scalar type parameter `name` with value `value` to the model `m`.
"""
@delegate set_external_scalar_param!(m::Model, name::Symbol, value::Any) => md

"""
    delete!(m::Model, component::Symbol

Delete a `component`` by name from a model `m`'s ModelDef, and nullify the ModelInstance.
"""
@delegate Base.delete!(m::Model, comp_name::Symbol) => md

"""
    set_param!(m::Model, comp_name::Symbol, name::Symbol, value, dims=nothing)

Set the parameter of a component `comp_name` in a model `m` to a given `value`.
The `value` can by a scalar, an array, or a NamedAray. Optional argument 'dims'
is a list of the dimension names of the provided data, and will be used to check
that they match the model's index labels.
"""
@delegate set_param!(m::Model, comp_name::Symbol, param_name::Symbol, value, dims=nothing) => md

"""
    set_param!(m::Model, path::AbstractString, param_name::Symbol, value, dims=nothing)

Set a parameter for a component with the given relative path (as a string), in which "/x" means the
component with name `:x` beneath `m.md`. If the path does not begin with "/", it is treated as
relative to `m.md`, which at the top of the hierarchy, produces the same result as starting with "/".
"""
@delegate set_param!(m::Model, path::AbstractString, param_name::Symbol, value, dims=nothing) => md

"""
    set_param!(m::Model, path::AbstractString, value, dims=nothing)

Similar to above but param_name appears in `path` after a colon delimiter.
"""
@delegate set_param!(m::Model, path::AbstractString, value, dims=nothing) => md

"""
    set_param!(m::Model, param_name::Symbol, value, dims=nothing)

Set the value of a parameter exposed in the ModelDef (m.md).
"""
@delegate set_param!(m::Model, param_name::Symbol, value, dims=nothing) => md

@delegate set_param!(m::Model, comp_path::ComponentPath, param_name::Symbol, value, dims=nothing) => md

"""
    run(m::Model)

Run model `m` once.
"""
function Base.run(m::Model; ntimesteps::Int=typemax(Int), rebuild::Bool=false,
                  dim_keys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)
    if length(m) == 0
        error("Cannot run a model with no components.")
    end

    if (rebuild || ! is_built(m))
        build(m)
    end

    # println("Running model...")
    mi = modelinstance(m)
    run(mi, ntimesteps, dim_keys)
    nothing
end

ext: jl


text:
#=
    The File contains `Sparsik<T>` and related funcs implementation

    `Sparsik` implements `AbstractSparseObject` interface
=#

#------------------------------------------------------------------------------

include("../typics.jl")
include("../gizmos.jl")

#------------------------------------------------------------------------------

import Base: ==, !=, +, -, *
import Nemo: GaloisField, GaloisFmpzField, fmpq
import Distributions: Bernoulli
import AbstractAlgebra: elem_type, Field, FieldElem,
                        characteristic, FracField

#------------------------------------------------------------------------------

"""
Sparsik
"""
mutable struct Sparsik{T<:Field} <: AbstractSparseVector{T}
    dim::Int
    field::T
    nonzero::Vector{Int}
    data::Dict{Int, <:FieldElem}
end

#------------------------------------------------------------------------------

# returns the density of `v`,
# i.e the percentage of nonzero-valued elements in `v`
density(v::Sparsik) = length(v.nonzero) / v.dim

# returns the number of nonzero-valued elements in `v`
nnz(v::Sparsik) = length(v.nonzero)

# base_ring field!
base_ring(v::Sparsik) = v.field

#------------------------------------------------------------------------------

# hitherto `Sparsik(other::Sparik)` has stood for the `deepcopy`,
# from now on we override Base.deepcopy
function Base.deepcopy_internal(x::Sparsik, stackdict::IdDict)
    y = Sparsik(x.dim, x.field,
            Base.deepcopy_internal(x.nonzero, stackdict),
            Base.deepcopy_internal(x.data, stackdict))
    stackdict[x] = y
    return y
end

#------------------------------------------------------------------------------

# returns the index of first nonzero or -1 if not present
# O(1)
function first_nonzero(v::Sparsik)
    if iszero(v)
        return -1
    end
    return first(v.nonzero)
end

#------------------------------------------------------------------------------

# returns v * c
# does modify `v`
# if `v` has k nonzeroes
# O(k)
function scale!(v::Sparsik, c)
    if iszero(base_ring(v)(c))
        return empty!(v)
    end

    for idx in v.nonzero
        v.data[idx] *= c
    end
    return v
end

# returns v * c
function scale(v::Sparsik, c)
    return scale!(deepcopy(v), c)
end

#------------------------------------------------------------------------------

# returns v + c * u
# mutates v
# if the vectors `v` and `u` have k and r nonzeroes respectively
# O(k + r)
function reduce!(v::Sparsik{T}, u::Sparsik{T}, c) where {T}
    new_nonzero = Int[]
    i, j = 1, 1

    v_nnz = deepcopy(v.nonzero)
    u_nnz = u.nonzero
    empty!(v.nonzero)

    while i <= length(v_nnz) || j <= length(u_nnz)
        new_idx = 0

        if i > length(v_nnz)
            new_idx = u_nnz[j]
            j += 1
        elseif j > length(u_nnz)
            new_idx = v_nnz[i]
            i += 1
        else
            new_idx = min(v_nnz[i], u_nnz[j])
            if v_nnz[i] < u_nnz[j]
                i += 1
            elseif v_nnz[i] > u_nnz[j]
                j += 1
            else
                i += 1
                j += 1
            end
        end
        v.data[new_idx] = v[new_idx] + c * u[new_idx]
        if !iszero(v[new_idx])
            push!(v.nonzero, new_idx)
        else
            delete!(v.data, new_idx)
        end
    end

    return v
end

# returns v + c * u
# O(k + r)
function Base.reduce(v::Sparsik{T}, u::Sparsik{T}, c) where {T}
    return reduce!(deepcopy(v), u, c)
end

#------------------------------------------------------------------------------

# returns standard inner product of v and u
# if the vectors `v` and `u` have k and r nonzeroes respectively
# O(min(k, r))
function inner(v::Sparsik{T}, u::Sparsik{T}) where {T}
    ans = zero(v.field)

    if length(v) > length(u)
        inner(u, v)
    end
    for i in v.nonzero
        # zero values can be skipped instantly
        # TODO: fix
        if haskey(u.data, i)
            ans += v.data[i] * u.data[i]
        end
    end

    return ans
end

function inner(A::Sparsik{T}, c::FieldElem) where {T}
    if iszero(c)
        return zero(c)
    end

    error("not implemented!")
end

#------------------------------------------------------------------------------

# constructs a `Sparsik` instance from a dense Vector
# and the provided field
# each `a` element is explicitly converted to its field representation
# O(n) where n = length(a)
function from_dense(a::AbstractVector, field)
    new_nonzero = Int[]
    new_data = Dict{Int, elem_type(field)}()

    for idx in findall(!(iszero ∘ field), a)
        push!(new_nonzero, idx[1])
        new_data[idx[1]] = field(a[idx])
    end

    return Sparsik(length(a), field, new_nonzero, new_data)
end

#------------------------------------------------------------------------------

# empties the `v`
# if the vector `v` has k nonzeroes
# O(k)
function Base.empty!(v::Sparsik)
    empty!(v.nonzero)
    empty!(v.data)
    return v
end

#------------------------------------------------------------------------------

# function for Sparsiks compatibility
# O(1)
function zero_sparsik(dim, field)
    return Sparsik(dim, field, Int[], Dict{Int, elem_type(field)}())
    # return zero(field)
end

#------------------------------------------------------------------------------

Base.zero(v::Sparsik) = zero_sparsik(v.dim, v.field)
Base.iszero(v::Sparsik{T}) where {T} = length(v) == 0

Base.get(v::Sparsik, i::Int) = get(v.data, i, zero(v.field))

function Base.haskey(A::Sparsik, i::Int)
    return haskey(A.data, i)
end

Base.getindex(v::Sparsik, i::Int) = get(v, i)

Base.size(v::Sparsik) = (v.dim, )
Base.size(v::Sparsik, i::Int) = size(v)[i]
Base.length(v::Sparsik) = length(v.nonzero)
dim(v::Sparsik) = v.dim

Base.isempty(v::Sparsik) = length(v.nonzero) == 0

# -----------------------------------------------------------------------------

# Strategy of types promotion:
# explicit restriction to the operands types
+(v::Sparsik{T}, u::Sparsik{T}) where {T} = reduce(v, u, one(v.field))
-(v::Sparsik{T}, u::Sparsik{T}) where {T} = reduce(v, u, - one(v.field))
-(v::Sparsik) = scale(v, - one(v.field))

*(v::Sparsik, c) = scale(v, c)
*(c, v::Sparsik) = v * c

==(v::Sparsik{T}, u::Sparsik{T}) where {T} = (v.dim == u.dim &&
            v.field == u.field &&
            v.data == u.data &&
            v.nonzero == u.nonzero)
!=(v::Sparsik{T}, u::Sparsik{T}) where {T} = !(v == u)

# -----------------------------------------------------------------------------

Base.repr(::MIME"text/plain", A::Sparsik) = "($(join(map(idx -> A[idx], 1 : A.dim), ", ")))"
Base.show(io::IO, A::Sparsik) = print(io, repr(MIME("text/plain"), A))

# -----------------------------------------------------------------------------

# iterates over all nonzeroes in `v`
# each returned element is a pair (idx, value)
# representing a nonzero element from `v`
# linear in number of nonzeroes
function Base.iterate(v::Sparsik, state=1)
    if state <= length(v)
        idx = v.nonzero[state]
        state += 1
        return ((idx, v[idx]), state)
    else
        return nothing
    end
end

# -----------------------------------------------------------------------------

Base.valtype(v::Sparsik) = elem_type(v.field)
Base.eltype(v::Sparsik) = (Int, valtype(v))

#-----------------------------------------------------------------------------

# Returns a new Sparsik object consisting of elements
# from the given `v` each converted to the `field`
# Throws if any of the denominators reduces to zero
#
# O(k) if k is the number of nonzero in `v`
function modular_reduction(v::Sparsik, field)
    new_nonzero = Int[]
    new_data = Dict{Int, elem_type(field)}()

    for (i, x) in v
        # wrapped into a function
        y = modular_reduction(x, field)
        if !iszero(y)
            push!(new_nonzero, i)
            new_data[i] = y
        end
    end

    return Sparsik(dim(v), field, new_nonzero, new_data)
end

#-----------------------------------------------------------------------------

# returns a new Sparsik object consisting of elements
# from the `v` each reconstructed from v.field to QQ
#
# O†(k) if k is the number of nonzero in `v`
# † reconstruction cost not included
function rational_reconstruction(v::Sparsik)
    new_nonzero = Int[]
    new_data = Dict{Int, fmpq}()

    tmp_type = BigInt
    if isa(base_ring(v), GaloisField)
        tmp_type = Int
    end
    ch = convert(tmp_type, characteristic(v.field))

    for (i, x) in v
        y = rational_reconstruction(convert(tmp_type, x), ch)
        if !iszero(y)
            push!(new_nonzero, i)
            new_data[i] = y
        end
    end

    return Sparsik(dim(v), QQ, new_nonzero, new_data)
end

#-----------------------------------------------------------------------------

function unit_sparsik(dim, i, field)
    # Int unfolding results into Int
    # Tuple{Int} unfolding results into Int
    # We can generalize it!
    return Sparsik(dim..., field, [i], Dict(i => one(field)))
end

#-----------------------------------------------------------------------------

#
function random_sparsik(sz::Int, field::T; density=0.1) where {T<:Field}
    if T <: FracField
        # todo!  ??????
        return unit_sparsik(sz, rand(1:sz), field)
    elseif T <: Union{GaloisField, GaloisFmpzField}
        dense_repr = map(field, rand(Bernoulli(density), sz))

        for idx in findall(!iszero, dense_repr)
            x = zero(field)
            while iszero(x)
                x = rand(field)
            end
            dense_repr[idx] = x
        end

        return from_dense(dense_repr, field)
    end

    error("random sparsik over $field is not implemented!")
end

function random_sparsik(sz::Tuple{Int}, field::T; density=0.1) where {T<:Field}
    random_sparsik(sz..., field, density=density)
end

#-----------------------------------------------------------------------------

# Inserts `value` to the `idx`-th position in `v`
# If idx is already present, rewrites it
# if k = length(v)
# O(klogk)
#
# it is not used
function Base.insert!(v::Sparsik, idx::Int, value)
    if iszero(value)
        return v
    end
    if iszero(v[idx])
        pos = searchsortedfirst(v.nonzero, idx)
        insert!(v.nonzero, pos, idx)
    end
    v.data[idx] = value
    return v
end

#-----------------------------------------------------------------------------

# Appends `value` to the vector treated as an array
# and increments the dimension of the sparsik by one
# O(1)
function Base.append!(v::Sparsik, value)
    v.dim += 1
    if iszero(value)
        return v
    end

    v.data[v.dim] = value
    push!(v.nonzero, v.dim)
    return v
end

# Erases the value at the last vector index
# and decrements the dimension of the sparsik by one
# Throws if v is 0-dim
# O(1)
function Base.pop!(v::Sparsik)
    if v.dim == 0
        error("sparsik $v is very very empty")
    end

    if haskey(v.data, v.dim)
        pop!(v.nonzero)
    end
    delete!(v.data, v.dim)
    v.dim -= 1
    return v
end

#-----------------------------------------------------------------------------

function restrict(v::Sparsik, coords)
    field = base_ring(v)
    nonzero = Int[]
    data = Dict{Int, elem_type(field)}()

    for (i, col) in enumerate(coords)
        if ! iszero(v[col])
            push!(nonzero, i)
            data[i] = v[col]
        end
    end

    return Sparsik(length(coords), field, nonzero, data)
end

#-----------------------------------------------------------------------------

function to_dense(A::Sparsik)
    ans = fill(zero(base_ring(A)), size(A))
    for (i, x) in A
        ans[i] = x
    end
    return ans
end

#-----------------------------------------------------------------------------

# v[idx] += c
# mutates v
# O(k)
# ( but O(1) generically )
function reduce!(v::Sparsik{T}, idx::Int, c) where {T}
    if iszero(c)
        return v
    end

    if haskey(v.data, idx)
        v.data[idx] += c
        if !iszero(v.data[idx])

        else
            popat!(v.nonzero, idx)
            delete!(v.data, idx)
        end
    else
        v.data[idx] = c
        pos = searchsortedfirst(v.nonzero, idx)
        insert!(v.nonzero, pos, idx)
    end
    return v
end


text:
###############################################################################
# Damper Force
###############################################################################
function damper_force(relative::Symbol, joint::Rotational{T}, qa::Quaternion, ωa::AbstractVector,
        qb::Quaternion, ωb::AbstractVector, timestep;
        rotate::Bool=true,
        unitary::Bool=false) where T

    damper = unitary ? 1.0 : joint.damper
    A = nullspace_mask(joint)
    Aᵀ = zerodimstaticadjoint(A)
    orientation_offset = joint.orientation_offset

    velocity = minimal_velocities(joint, szeros(3), szeros(3), qa, ωa, szeros(3), szeros(3), qb, ωb, timestep)
    if relative == :parent
        force = 1.0000 * damper * Aᵀ * velocity # currently assumes same damper constant in all directions
        rotate && (force = vector_rotate(force, orientation_offset)) # rotate back to frame a
        return [szeros(T, 3); force]
    elseif relative == :child
        force = - 1.0000 * damper * Aᵀ * velocity # currently assumes same damper constant in all directions
        rotate && (force = vector_rotate(force, inv(qb) * qa * orientation_offset)) # rotate back to frame b
        return [szeros(T, 3); force]
    end
end

damper_impulses(relative::Symbol, joint::Rotational, pbody::Node, cbody::Node, timestep; unitary::Bool=false) =
    timestep * damper_force(relative, joint, current_configuration(pbody.state)[2], pbody.state.ωsol[2],
    current_configuration(cbody.state)[2], cbody.state.ωsol[2], timestep, unitary=unitary)

damper_impulses(relative::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep; unitary::Bool=false) where T = szeros(T, 6)

###############################################################################
# Damper Jacobians
################################################################################

function damper_jacobian_configuration(relative::Symbol, jacobian::Symbol,
        joint::Rotational, pbody::Node, cbody::Node,
        timestep::T;
        ) where T

    Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
    xa, va, qa, ωa = current_configuration_velocity(pbody.state)
    xb, vb, qb, ωb = current_configuration_velocity(cbody.state)
    orientation_offset = joint.orientation_offset
    X = szeros(T, 3, 3)
    Z = szeros(T, 3, 6)

    force = damper_force(relative, joint, qa, ωa, qb, ωb, timestep; rotate = false)[SUnitRange(4,6)]
    ∂vel = minimal_velocities_jacobian_configuration(jacobian, joint, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)[:, SUnitRange(4,6)]

    if relative == :parent
        Q = rotation_matrix(orientation_offset) * 1.0000 * joint.damper * Aᵀ * ∂vel
    elseif relative == :child
        Q = rotation_matrix(inv(qb) * qa * orientation_offset) * -1.0000 * joint.damper * Aᵀ * ∂vel
        if jacobian == :parent
            # Q += rotation_matrix(inv(qb)) * ∂rotation_matrix∂q(qa, vector_rotate(force, orientation_offset), attjac=true)
            Q += rotation_matrix(inv(qb)) * ∂rotation_matrix∂q(qa, vector_rotate(force, orientation_offset)) * LVᵀmat(qa) # ATTJAC
        elseif jacobian == :child
            # Q += ∂rotation_matrix_inv∂q(qb, vector_rotate(force, qa * orientation_offset), attjac=true)
            Q += ∂rotation_matrix_inv∂q(qb, vector_rotate(force, qa * orientation_offset)) * LVᵀmat(qb) # ATTJAC
        end
    end
    return timestep * [Z; X Q]
end

function damper_jacobian_velocity(relative::Symbol, jacobian::Symbol,
    joint::Rotational, pbody::Node, cbody::Node,
    timestep::T) where T

    Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
    xa, va, qa, ωa = current_configuration_velocity(pbody.state)
    xb, vb, qb, ωb = current_configuration_velocity(cbody.state)
    orientation_offset = joint.orientation_offset

    force = damper_force(relative, joint, qa, ωa, qb, ωb, timestep; rotate = false)[SUnitRange(4,6)]
    ∂vel = minimal_velocities_jacobian_velocity(jacobian, joint, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)

    if relative == :parent
        VΩ = rotation_matrix(orientation_offset) * 1.0000 * joint.damper * Aᵀ * ∂vel
    elseif relative == :child
        VΩ = rotation_matrix(inv(qb) * qa * orientation_offset) * -1.0000 * joint.damper * Aᵀ * ∂vel
    end
    return timestep * [szeros(T, 3, 6); VΩ]
end

damper_jacobian_configuration(relative::Symbol, jacobian::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep::T; attjac::Bool = true) where T = attjac ? szeros(T, 6, 6) : szeros(T, 6, 7)
damper_jacobian_velocity(relative::Symbol, jacobian::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep::T) where T = szeros(T, 6, 6)


text:
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test

#this function was removed in Julia 0.7
"""    linreg(x,y)

perform a simple linear regression for scalar case. Coefficients are returned
as (intercept, scalar)
"""
linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y

"""    peakgap(s1, s2)

identify the time delay between the peak in two time series. Time series should
be stored as (x(t), t) for consistency with the tuples(sol) from OrdinaryDiffEq.
"""
function peakgap(s1, s2)
    @show p1, i1 = findmax(first.(s1))
    @show p2, i2 = findmax(first.(s2))
    return s2[i2][end] - s1[i1][end]
end

"""    peakgap(sol, i::Int, j::Int)

identify the time delay between the peak in two dimensions of an ODE solution structure.
"""
function peakgap(sol, i::Int, j::Int)
    s1 = [(u[i],t) for (u,t) in tuples(sol)]
    s2 = [(u[j],t) for (u,t) in tuples(sol)]
    return peakgap(s1,s2)
end

"""    paramsweep(f::Function, m::Petri.Model, initials, tspan, params)

solve an ODEProblem for each parameter setting in params and apply the function f to the solution.

Usage:
    paramsweep(m,initials,tspan,[p1,p2,p3]) do sol
        return f(sol)
    end
"""
function paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
    map(params) do p
        prob = ODEProblem(fluxes(m), initials, tspan, p)
        sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
        return f(sol)
    end
end

"""    fluxes(m::Petri.Model)

a PetriNet interpreter that computes the mass action kinetics of a petri net.

Usage: pass this to ODEProblem to set up an ODE for a given model.
"""
function fluxes(m::Petri.Model)
    S = m.S
    T = m.Δ
    nS = length(S)
    nT = length(T)
    ϕ = zeros(Float64, nT)
    f(du, u, p, t) = begin
        for (i, t) in enumerate(T)
            ins = t[1]
            # TODO: accomodate multiplicites here
            ϕ[i] = p[i]*prod(u[ins])
        end
        for i in S
            du[i] = 0
        end
        for (i, t) in enumerate(T)
            ins = t[1]
            out = t[2]
            for s in ins
                # TODO: accomodate multiplicites here
                du[s] -= ϕ[i]
            end
            for s in out
                # TODO: accomodate multiplicites here
                du[s] += ϕ[i]
            end
        end
        return du
    end
    return f
end

u₀(m::Petri.Model) = begin
    zeros(Float64, length(m.S))
end

u₀(m::Petri.Model, initialS, initialI=1) = begin
    u0=zeros(Float64, length(m.S))
    u0[1] = initialS
    u0[2] = initialI
    return u0
end

function makeplots_seir(sol, prefix)
    mkpath(dirname(prefix))
    p1 = plot(sol,vars=[1,2,3,4], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
    p2 = plot(sol,vars=[5,6,7,8], xlabel="", ylabel="people", linewidth=3, legend=false)
    p3 = plot(sol,vars=[9,10,11,12], xlabel="time", ylabel="people", linewidth=3, legend=false)
    p4 = plot(sol,vars=[2,6,10], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
    p5 = plot(sol,vars=[3,7,11], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
    p6 = plot(sol,vars=[4,8,12], xlabel="time", linewidth=3, labels=["r1" "r2" "r3"], legend=true)
    p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
    savefig(p, "$(prefix)combined.pdf")
    p
end

function makeplots_seird(sol, prefix)
    mkpath(dirname(prefix))
    p1 = plot(sol,vars=[1,2,3,4,5], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
    p2 = plot(sol,vars=[6,7,8,9,10], xlabel="", ylabel="people", linewidth=3, legend=false)
    p3 = plot(sol,vars=[11,12,13,14,15], xlabel="time", ylabel="people", linewidth=3, legend=false)
    p4 = plot(sol,vars=[2,7,12], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
    p5 = plot(sol,vars=[3,8,13], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
    p6 = plot(sol,vars=[5,10,15], xlabel="time", linewidth=3, labels=["d1" "d2" "d3"], legend=true)
    p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
    savefig(p, "$(prefix)combined.pdf")
    p
end

T = [
    # City 1 SEIR
    ([1, 2], [3, 2]),
    ([3], [2]), # E→I
    ([2], [4]), # I→R
    # outflow 1→2
    ([1], [5]), # S→S′
    ([2], [6]), # I→I′
    ([3], [7]), # E→E′
    # City 2 SEIR
    ([5, 6], [7, 6]),
    ([7], [6]), # E→I
    ([6], [8]), # I→R
    # outflow 2→3
    ([5], [9]), # S→S′
    ([6], [10]),# I→I′
    ([7], [11]),# E→E′
    # City 3 SEIR
    ([9, 10], [11, 10]),
    ([11], [10]), # E→I
    ([10], [12]) # I→R
]

m = Petri.Model(1:12, T, missing, missing)
nS = length(m.S)
β = ones(Float64, length(T))
@test fluxes(m)(zeros(Float64, nS), ones(Float64, nS), β, 1) |> length == length(m.S)

tspan = (0.0,60.0)
prob = ODEProblem(fluxes(m), u₀(m, 1,0), tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@test sol.u[end][end-3] > .90
@test sol.u[end][end-2] == 0
@test sol.u[end][end-1] == 0
@test sol.u[end][end] == 0

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[5]  = 10000
u0[9]  = 10000
u0[2]  = 1
u0
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/1000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]

using Plots
makeplots_seir(sol, "img/seir/baseline")

βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/100
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/travel")

βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/200, 1/2]/100
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/screening")

βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/200, 1/200, 1/200]/100000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seir(sol, "img/seir/shutdown")

println("Exploring the SEIRD^3 example")

S = 1:15
T = [([1, 2], [3, 2]),
     ([3], [2]),
     ([2], [4]),
     ([2], [5]),
     ([1], [6]),
     ([2], [7]),
     ([3], [8]),
     ([6, 7], [8, 7]),
     ([8], [7]),
     ([7], [9]),
     ([7], [10]),
     ([6], [11]),
     ([7], [12]),
     ([8], [13]),
     ([11, 12], [13, 12]),
     ([13], [12]),
     ([12], [14]),
     ([12], [15])]
m = Petri.Model(S, T, missing, missing)
u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[6]  = 10000
u0[11]  = 10000
u0[2]  = 1
u0
βseir = [10/sum(u0), 1/2, 1/5, 1/10]
βtravel = [1/2, 1/2, 1/2]/1000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]
makeplots_seird(sol, "img/seird/baseline")

function paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
    map(params) do p
        prob = ODEProblem(fluxes(m), initials, tspan, p)
        sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
        return f(sol)
    end
end

packparams(βlocal, βflow) = vcat(βlocal, βflow, βlocal, βflow, βlocal)

# The βseir parameters are the mass action kinetics rates of
# 1. S+I --> E+I
# 2. E --> I
# 3. I --> R
# The βtravel parameters are the mass action kinetics rates of outflow to the next
# city in the transit network.
# 1. S --> S⁺
# 2. I --> I⁺
# 3. E --> E⁺

# Baseline scenario

βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/100
β1 = packparams(βseir, βtravel)

# Travel scenario
βtravel = [1/2, 1/2, 1/2]/10
β2 = packparams(βseir, βtravel)

# Screening scenario
βtravel = [1/2, 1/20, 1/2]/100
β3 = packparams(βseir, βtravel)

# Shutdown scenario
βtravel = [1/2, 1/2, 1/2]/1000
β4 = packparams(βseir, βtravel)

# Total and Complete Shutdown scenario
βtravel = [1/2, 1/2, 1/2]/10000

β5 = packparams(βseir, βtravel)
βs = [β1, β2, β3, β4, β5]

T = [
    # City 1 SEIR
    ([1, 2], [3, 2]),
    ([3], [2]), # E→I
    ([2], [4]), # I→R
    # outflow 1→2
    ([1], [5]), # S→S′
    ([2], [6]), # I→I′
    ([3], [7]), # E→E′
    # City 2 SEIR
    ([5, 6], [7, 6]),
    ([7], [6]), # E→I
    ([6], [8]), # I→R
    # outflow 2→3
    ([5], [9]), # S→S′
    ([6], [10]),# I→I′
    ([7], [11]),# E→E′
    # City 3 SEIR
    ([9, 10], [11, 10]),
    ([11], [10]), # E→I
    ([10], [12]) # I→R
]

# u0 is the initial configuration of the system
# u0[[1,5,9]] = number of initially Susceptible people in each population
# u0[2] = number of Infected individuals
u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[5]  = 10000
u0[9]  = 10000
u0[2]  = 1
u0

m = Petri.Model(1:12, T, missing, missing)
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, βs)

# prob = ODEProblem(fluxes(m), u0, tspan, βs[1])
# sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
# makeplots_seir(sol, "img/debug_1")

travelparams(k::Number) = begin
    βseir = [10/sum(u0), 1/2, 1/5]
    βtravel = [1/2, 1/2, 1/2]/(100k)
    β = packparams(βseir, βtravel)
    return β
end

travelsweep(r::AbstractVector) = begin
    return travelparams.(2 .^ r)
end

βrange = 1:8
gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, travelsweep(βrange))

coeffs = linreg(βrange, gaps)
@test sum(abs.((coeffs[1] .+ (coeffs[2]*(βrange))) .- gaps))/length(βrange)<= 1e-1

# struct GapResult{P,G}
#     param::P
#     gap::G
# end
#
# GapResult(t::Tuple) = GapResult(t...)
#
# struct GapStudy
#     result::Vector
#     coeffs::Vector
# end
#
# slope(gs::GapStudy) = gs.coeffs[end]

function estimatedelay(r::AbstractVector)
    gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, travelsweep(r))
    coeffs = linreg(βrange, gaps)
    return (results=zip(r, gaps), coeffs=coeffs)
end

function estimatedelay(m::Petri.Model, u0, tspan, parameters)
    gaps = paramsweep(sol->peakgap(sol, 2, 10), m, u0, tspan, parameters)
    coeffs = linreg(βrange, gaps)
    return (results=zip(parameters, gaps), coeffs=coeffs)
end

slope(gs) = gs.coeffs[end]

gs = estimatedelay(βrange)
@test slope(gs) > 1
collect(gs.results)
gs.coeffs

using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
using SemanticModels.CategoryTheory
import SemanticModels.CategoryTheory: undecorate, ⊔
using SemanticModels.PetriModels
using SemanticModels.PetriCospans
import SemanticModels: model

model(c::PetriCospan) = left(c.f).d[1].model

X = FinSet(1)
Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
seir_petri = left(Fseir.f).d[1]
f = FinSetMorph(1:4, [1, 2, 3])
g = FinSetMorph(1:4, [1, 2, 3])
Fseir′ = PetriCospan(Cospan(Decorated(f,seir_petri), Decorated(g, seir_petri)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
@test length(model(Fseir′).S) == 4
@test length(model(Fseir′).Δ) == 3

m = model(Fseir′ ⋅ Fflow ⋅ Fseir′ ⋅ Fflow ⋅ Fseir′)

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[5]  = 10000
u0[9]  = 10000
u0[2]  = 1
res = estimatedelay(m::Petri.Model, u0, tspan, travelsweep(1:8))
seirgaps = collect(map(last, res.results))

@test 1.35 <= slope(res) <= 1.45

Pseird = PetriModel(
          Petri.Model(1:5,[
            ([1,2],[3,2]), # exposure
            ([3],[2]),     # onset
            ([2],[4]),     # recovery
            ([2],[5]),     # death
            ], missing, missing))
inputs = FinSetMorph(1:5, [1,2,3])
outputs = FinSetMorph(1:5, [1,2,3])
Fcityd = PetriCospan(Cospan(Decorated(inputs, Pseird),
                            Decorated(outputs, Pseird)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
Fcity₀ = Fcityd ⋅ Fflow
Fcity₁ = Fcityd ⋅ Fflow
Fcityₑ = Fcityd
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ

m = model(Fcity³)

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[6]  = 10000
u0[11] = 10000
u0[2]  = 1
seirdparams(k) = begin
    βseird = [10/sum(u0), 1/2, 1/5, 1/16]
    βtravel = [1/2, 1/2, 1/2]/(100k)
    β = vcat(βseird, βtravel)
    return vcat(β, β, β)
end

gaps = paramsweep(sol->peakgap(sol, 2, 12), m,
                  u0, tspan, seirdparams.(2 .^ (1:8)))

prob = ODEProblem(fluxes(model(Fcity³)), u0, tspan, seirdparams(2))
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())

seirdgaps = gaps

@test model(Fcity³).S == 1:15
@test length(model(Fcity³).Δ) == 18

makeplots_seird(sol, "img/tmp")

gp = plot([seirgaps seirdgaps], marker=:square, labels=[:seir :seird], xlabel="2^x fold reduction in travel", ylabel="Gap between first and last epidemic", legend=:bottomright, title="Effect of model assumptions on travel restrictions")
savefig(gp, "img/gapplot.pdf")


text:
function add_contact_mark_callbacks(b,handles)

    signal_connect(contact_onset_cb,b["contact_block_button"],"toggled",Void,(),false,(handles,))
    signal_connect(no_contact_cb,b["no_contact_block_button"],"toggled",Void,(),false,(handles,))

    signal_connect(contact_select_cb,b["contact_spin"],"value-changed",Void,(),false,(handles,))

    signal_connect(pro_re_cb,b["protraction_button"],"clicked",Void,(),false,(handles,1))
    signal_connect(pro_re_cb,b["retraction_button"],"clicked",Void,(),false,(handles,2))

    signal_connect(exclude_cb,b["exclude_button"],"toggled",Void,(),false,(handles,))

    signal_connect(advance_slider_cb,b["contact_win"],"key-press-event",Void,(Ptr{Gtk.GdkEventKey},),false,(handles,))

    nothing
end

function pro_re_cb(w::Ptr,user_data::Tuple{Tracker_Handles,Int})

    han, con_type = user_data

    han.man.pro_re[han.displayed_frame]=con_type

    update_pro_re_block(han.man)

    nothing
end

function update_pro_re_block(man::Manual_Class)

    #https://stackoverflow.com/questions/29848734/is-it-possible-to-sort-a-dictionary-in-julia
    sorted_pro_re=sort(collect(man.pro_re), by=x->x[1])

    for i=1:(length(sorted_pro_re)-1)
        man.pro_re_block[sorted_pro_re[i][1]] = sorted_pro_re[i][2]

        if sorted_pro_re[i+1][2] == sorted_pro_re[i][2]
            man.pro_re_block[sorted_pro_re[i][1]:sorted_pro_re[i+1][1]] .= sorted_pro_re[i][2]
        end
    end

    nothing
end

function contact_onset_cb(w::Ptr,user_data::Tuple{Tracker_Handles})
    han, = user_data

    onset = get_gtk_property(han.b["contact_block_button"],:active,Bool)

    if onset

        han.man.partial_contact = han.displayed_frame

        set_gtk_property!(han.b["contact_block_button"],:label,"Mark Contact Offset")
    else

        if han.displayed_frame > han.man.partial_contact
            push!(han.man.contact_block,(han.man.partial_contact,han.displayed_frame))
            update_contact_block(han.man)
            calc_contact_block(han)
        else
            println("error: start with beginning of contact block")
        end
        set_gtk_property!(han.b["contact_block_button"],:label,"Mark Contact Onset")
    end

    nothing
end

function update_contact_block(man::Manual_Class)

    myrange = make_range(man.contact_block[end])

    man.contact[myrange] .= 2 #Contact

    man.contact_block = sort(man.contact_block, by=x->x[1])

    nothing
end

function no_contact_cb(w::Ptr,user_data::Tuple{Tracker_Handles})
    han, = user_data

    onset = get_gtk_property(han.b["no_contact_block_button"],:active,Bool)

    if onset

        han.man.partial_contact = han.displayed_frame

        set_gtk_property!(han.b["no_contact_block_button"],:label,"Mark End of Block")
    else

        if han.displayed_frame > han.man.partial_contact
            push!(han.man.no_contact_block,(han.man.partial_contact,han.displayed_frame))
            update_no_contact_block(han.man)
            calc_contact_block(han)
        else
            println("error: start with beginning of contact block")
        end
        set_gtk_property!(han.b["no_contact_block_button"],:label,"Mark No Contact")
    end
    nothing
end

function update_no_contact_block(man::Manual_Class)

    myrange = make_range(man.no_contact_block[end])

    man.contact[myrange] .= 1 #No Contact

    man.no_contact_block = sort(man.no_contact_block, by=x->x[1])
end

function exclude_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    onset = get_gtk_property(han.b["exclude_button"],:active,Bool)

    if onset

        han.man.partial_contact = han.displayed_frame

        set_gtk_property!(han.b["exclude_button"],:label,"Mark End of Block")
    else

        if han.displayed_frame > han.man.partial_contact
            push!(han.man.exclude,(han.man.partial_contact,han.displayed_frame))
            update_exclude(han.man)
        else
            println("error: start with beginning of contact block")
        end
        set_gtk_property!(han.b["exclude_button"],:label,"Start Exclude")
    end

    nothing
end

function update_exclude(man::Manual_Class)

    myrange = make_range(man.exclude[end])

    man.exclude_block[myrange] .= true #Exclude Frames

    man.exclude = sort(man.exclude, by=x->x[1])

    nothing
end

function save_manual_class(man::Manual_Class,path)

    file = jldopen(path,"w")
    write(file,"Max_Frames",man.max_frames)
    write(file,"Contact_Block",man.contact_block)
    write(file,"No_Contact_Block",man.no_contact_block)
    write(file,"Contact",man.contact)

    write(file,"Pro_Re",man.pro_re)
    write(file,"Pro_Re_Block",man.pro_re_block)

    write(file, "Exclude",man.exclude)
    write(file, "Exclude_Block",man.exclude_block)

    close(file)
    nothing
end

function load_manual_class(man::Manual_Class,path)

    file = jldopen(path,"r")
    man.max_frames = read(file,"Max_Frames")

    man.contact_block = read(file,"Contact_Block")
    man.no_contact_block = read(file,"No_Contact_Block")
    man.contact = read(file,"Contact")

    man.pro_re = read(file,"Pro_Re")
    man.pro_re_block = read(file,"Pro_Re_Block")

    man.exclude = read(file, "Exclude")
    man.exclude_block = read(file, "Exclude_Block")

    close(file)
    nothing
end

function make_range(x::Tuple)
    x[1]:x[2]
end

function draw_touch2(han::Tracker_Handles)

    if length(han.tracked_contact) == han.max_frames

        ctx=Gtk.getgc(han.c)

        if han.man.contact[han.displayed_frame] != 0 #prefer manual
            if han.man.contact[han.displayed_frame] == 2
                set_source_rgb(ctx,1,0,0)
            else
                set_source_rgb(ctx,1,1,1)
            end
        else
            if han.tracked_contact[han.displayed_frame]
                set_source_rgb(ctx,1,0,0)
            else
                set_source_rgb(ctx,1,1,1)
            end
        end

        rectangle(ctx,620,0,20,20)
        fill(ctx)

        for i=1:length(han.man.calc_contact_block)
            if (han.displayed_frame >= han.man.calc_contact_block[i][1]) & (han.displayed_frame <= han.man.calc_contact_block[i][2])
                set_source_rgb(ctx,0,0,0)
                move_to(ctx,620,15)
                show_text(ctx,string(i))
            end
        end

        reveal(han.c)
    end

    nothing
end

function contact_select_cb(w::Ptr,user_data::Tuple{Tracker_Handles})
    han, = user_data

    c_id = get_gtk_property(han.b["contact_spin"],:value,Int64)

    if (c_id > 0)&(c_id <= length(han.man.calc_contact_block))
        frame = han.man.calc_contact_block[c_id][1]
        set_gtk_property!(han.b["adj_frame"],:value,frame)
    end

    nothing
end

function calc_contact_block(han::Tracker_Handles)

    if length(han.tracked_contact) == han.max_frames
        contact = falses(han.max_frames)
        for i=1:length(han.tracked_contact)
            if han.man.contact[i] == 2
                contact[i] = true
            elseif han.man.contact[i] == 1
                contact[i] = false
            else
                contact[i] = han.tracked_contact[i]
            end
        end

        (c_on,c_off) = get_contact_indexes(contact,han.man.exclude_block)
        han.man.calc_contact_block = [(c_on[i],c_off[i]) for i=1:length(c_on)]

        set_gtk_property!(han.b["contact_number_label"],:label,string(length(c_on)))
        set_gtk_property!(han.b["contact_spin_adj"],:upper,length(c_on))
    end
end

function get_contact_indexes(c,exclude,c_dur_min = 4)

    c_ind=Array{Int64,1}()
    c_off=Array{Int64,1}()

    i=1
    while (i < length(c) - c_dur_min)

        if exclude[i] != true
            if ((c[i] - c[i+1])==-1)
                push!(c_ind,i)
                myoff=findnext(c.==false,i+1)
                if typeof(myoff) == Nothing #single frame contact
                    myoff = i+1
                end

                push!(c_off,myoff)

                i=myoff + 1
            else
                i+=1
            end
        else
            i+= 1
        end
    end

    keep=trues(length(c_ind))
    for i=1:length(c_ind)
        if (c_off[i] - c_ind[i]) < c_dur_min
            keep[i] = false
        end
    end

    (c_ind[keep], c_off[keep])
end


text:
# This file was generated, do not modify it. # hide
mean_age = mean(skipmissing(data_nmiss.plant_age))
median_age = median(skipmissing(data_nmiss.plant_age))


text:
module ValueToState

function getState(value_str, downregulatedCutoff, upregulatedCutoff)
    value = typeof(value_str) == Float64 ? value_str : parse(Float64, value_str)
    if upregulatedCutoff < value
        return "Up Regulated"
    elseif downregulatedCutoff > value
        return "Down Regulated"
    else
        return "Normal"
    end
end

function getStateNumber(value_str, downregulatedCutoff, upregulatedCutoff)
    value = parse(Float64, value_str)
    if upregulatedCutoff < value
        return "2"
    elseif downregulatedCutoff > value
        return "0"
    else
        return "1"
    end
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
# Figure 5.3b in "Infinite GMRES for parameterized linear systems"
using SparseArrays
using LinearAlgebra
using PyPlot
using MAT
pygui(true)
include("iar.jl");
include("infgmres.jl");

# Problem represented with explicitly available derivatives
struct GeneralProblem
    sigma
    Aders
    Afun
    b
end

function compute_Mlincomb(prob::GeneralProblem,sigma,Y)
    if (prob.sigma != sigma)
        @show sigma
        @show prob.sigma
        error("Wrong sigma");
    end
    z=zeros(eltype(Y),size(Y,1));
    for i=1:size(Y,2)
        z += prob.Aders[i]*Y[:,i];
    end
    return z;
end
function compute_Mder(prob::GeneralProblem,mu)
    return prob.Afun(mu);
end

## Define problem
data=matread("delay_problem_matrices.mat")
A0=data["A0"];
A1=data["A1"];
b=vec(data["b"]);
n=length(b);

# Identity matrix
II=Matrix{Float64}(I,n,n);

Afun= s-> -s*II+A0+A1*exp(-s);
sigma=0;
N=100;
ders=Vector{Matrix{Float64}}(undef,N);

for i=1:N
    ders[i]=((-1)^(i+1))*A1;
end
ders[1] += A0;
ders[2] += -II;

myproblem=GeneralProblem(sigma,ders,Afun,b);
(f,hist1)=infgmres(myproblem,30,method=:iar);
max_it = 20;
v = range(-5,stop=-2,length=50);
mu_list = [10^i for i in v];
tol1 = 1.e-10;
obs = zeros(length(mu_list),1);
pred1 = zeros(length(mu_list),1);
r1 = 9.304668750449350; #This is taken from TDS Arnoldi

for j = 1:length(mu_list)
    global mu
    res = zeros(max_it,1);
    mu = mu_list[j];
    A_of_mu = myproblem.Afun(mu);

    for m=1:max_it
        global f
        (xx)=f(mu,m)
        res[m] = norm(A_of_mu*xx - myproblem.b)/norm(myproblem.b);
    end

    #Take the worst reducation in the residual
    ratio_list=zeros(max_it,1);
    end_m=0;
    for i=1:length(res)
        if res[i] < tol1
            end_m = i;
            break;
        end
    end
    for i=2:end_m
        ratio_list[i-1]=res[i]/res[i-1];
    end
    obs[j]=maximum(ratio_list);
    pred1[j]=mu*r1;

end

figure(1)
ax = PyPlot.gca()
ax.set_yscale("log")
ax.set_xscale("log")
xlabel("mu")
ylabel("rho")
plot(mu_list,obs,"o")
plot(mu_list,pred1)
legend(["Observed","Predicted"])

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
using Clang.Generators

cd(@__DIR__)

const HEADER_BASE = joinpath(ENV["KNITRODIR"], "include")
const KNITRO_H = joinpath(HEADER_BASE, "knitro.h")

headers = [KNITRO_H]
options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$HEADER_BASE")

ctx = create_context(headers, args, options)
build!(ctx)


text:
using Libxc
include("../xc/xc_evaluate.jl")

"""
Exchange-correlation term, defined by a list of functionals and usually evaluated through libxc.
"""
struct Xc
    functionals::Vector{Symbol}  # Symbols of the functionals (Libxc.jl / libxc convention)
    scaling_factor::Real         # Scales by an arbitrary factor (useful for exploration)

    # Density cutoff for XC computation: Below this value a gridpoint counts as zero
    # `nothing` implies that libxc defaults are used (for each functional a different
    # small positive value like 1e-24)
    density_threshold::Union{Nothing,Float64}
end
Xc(symbols::Symbol...; kwargs...) = Xc([symbols...]; kwargs...)
function Xc(symbols::Vector; scaling_factor=1, density_threshold=nothing)
    Xc(convert.(Symbol, symbols), scaling_factor, density_threshold)
end

function (xc::Xc)(basis)
    functionals = Functional.(xc.functionals; n_spin=basis.model.n_spin_components)
    if !isnothing(xc.density_threshold)
        for func in functionals
            func.density_threshold = xc.density_threshold
        end
    end
    TermXc(basis, functionals, xc.scaling_factor)
end

struct TermXc <: Term
    basis::PlaneWaveBasis
    functionals::Vector{Functional}
    scaling_factor::Real
end

@timing "ene_ops: xc" function ene_ops(term::TermXc, ψ, occ; ρ, ρspin=nothing, kwargs...)
    basis = term.basis
    T     = eltype(basis)
    model = basis.model
    n_spin = model.n_spin_components
    @assert all(xc.family in (:lda, :gga) for xc in term.functionals)

    if isempty(term.functionals)
        ops = [NoopOperator(term.basis, kpoint) for kpoint in term.basis.kpoints]
        return (E=0, ops=ops)
    end

    # Take derivatives of the density if needed.
    max_ρ_derivs = maximum(max_required_derivative, term.functionals)
    density = LibxcDensity(basis, max_ρ_derivs, ρ, ρspin)

    potential = [zeros(T, basis.fft_size) for _ in 1:n_spin]  # TODO CPU arrays
    zk = zeros(T, basis.fft_size)  # Energy per unit particle
    E = zero(T)
    for xc in term.functionals
        # Evaluate the functional and its first derivative (potential)
        terms = evaluate(xc, density; zk=zk)

        # Add energy contribution
        dVol = basis.model.unit_cell_volume / prod(basis.fft_size)
        E += sum(terms.zk .* ρ.real) * dVol

        # Add potential contributions Vρ -2 ∇⋅(Vσ ∇ρ)
        for σ in 1:n_spin
            potential[σ] .+= @view terms.vrho[σ, :, :, :]
        end

        if haskey(terms, :vsigma)
            # TODO Potential to save some FFTs for spin-polarised calculations:
            #      For n_spin == 2 this calls divergence_real 4 times, where only 2 are
            #      needed (one for potential[1] and one for potential[2])
            @views for σ in 1:n_spin, τ in 1:n_spin
                στ = libxc_symmetric_index(σ, τ)
                # Extra factor (1/2) for σ != τ is needed because libxc only keeps σ_{αβ}
                # in the energy expression. See comment block below on spin-polarised XC.
                spinfac = (σ == τ ? 2 : 1)
                potential[σ] .+=
                     -spinfac .* divergence_real(
                        α -> terms.vsigma[στ, :, :, :] .* density.∇ρ_real[τ, :, :, :, α],
                        density.basis,
                    )
            end
        end
    end
    if term.scaling_factor != 1
        E *= term.scaling_factor
        potential = [pot .*= term.scaling_factor for pot in potential]
    end

    ops = [RealSpaceMultiplication(basis, kpoint, potential[kpoint.spin])
           for kpoint in basis.kpoints]
    (E=E, ops=ops)
end

function compute_kernel(term::TermXc; ρ::RealFourierArray, ρspin=nothing, kwargs...)
    @assert term.basis.model.spin_polarization in (:none, :spinless, :collinear)
    density = LibxcDensity(term.basis, 0, ρ, ρspin)
    n_spin = term.basis.model.n_spin_components
    n_fxc  = n_spin == 1 ? 1 : 3  # Number of spin components in fxc
    @assert 1 ≤ n_spin ≤ 2

    kernel = similar(density.ρ_real, n_fxc, term.basis.fft_size...)
    kernel .= 0
    for xc in term.functionals
        xc.family == :lda || error("compute_kernel only implemented for LDA")
        terms = evaluate(xc, density; derivatives=2:2)  # only valid for LDA
        kernel .+= terms.v2rho2
    end

    fac = term.scaling_factor
    if n_spin == 1
        Diagonal(vec(fac .* kernel))
    else
        # Blocks in the kernel matrix mapping (ρα, ρβ) ↦ (Vα, Vβ)
        Kαα = @view kernel[1, :, :, :]
        Kαβ = @view kernel[2, :, :, :]
        Kβα = Kαβ
        Kββ = @view kernel[3, :, :, :]

        # Blocks in the kernel matrix mapping (ρtot, ρspin) ↦ (Vα, Vβ)
        K_αtot  = Diagonal(vec(fac * (Kαα + Kαβ) / 2))
        K_αspin = Diagonal(vec(fac * (Kαα - Kαβ) / 2))
        K_βtot  = Diagonal(vec(fac * (Kβα + Kββ) / 2))
        K_βspin = Diagonal(vec(fac * (Kβα - Kββ) / 2))

        [K_αtot K_αspin;
         K_βtot K_βspin]
    end
end

function apply_kernel(term::TermXc, dρ::RealFourierArray, dρspin=nothing;
                      ρ::RealFourierArray, ρspin=nothing, kwargs...)
    basis  = term.basis
    T      = eltype(basis)
    n_spin = basis.model.n_spin_components
    isempty(term.functionals) && return nothing
    @assert all(xc.family in (:lda, :gga) for xc in term.functionals)
    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)

    # Take derivatives of the density and the perturbation if needed.
    max_ρ_derivs = maximum(max_required_derivative, term.functionals)
    density      = LibxcDensity(basis, max_ρ_derivs, ρ, ρspin)
    perturbation = LibxcDensity(basis, max_ρ_derivs, dρ, dρspin)

    result = similar(ρ.real, basis.fft_size..., n_spin)
    result .= 0
    for xc in term.functionals
        # TODO LDA actually only needs the 2nd derivatives for this ... could be optimised
        terms = evaluate(xc, density, derivatives=1:2)

        # Accumulate LDA and GGA terms in result
        @views for σ in 1:n_spin, τ in 1:n_spin
            στ = libxc_symmetric_index(σ, τ)
            result[:, :, :, σ] .+= terms.v2rho2[στ, :, :, :] .* perturbation.ρ_real[τ, :, :, :]
        end
        if haskey(terms, :v2rhosigma)
            @assert basis.model.spin_polarization in (:none, :spinless)
            result[1] .+= apply_kernel_term_gga(terms, density, perturbation)
        end
    end

    [from_real(basis, term.scaling_factor .* result[:, :, :, σ]) for σ in 1:n_spin]
end

#=
The total energy is
Etot = ∫ ρ E(ρ,σ), where σ = |∇ρ|^2 and E(ρ,σ) is the energy per unit particle.
libxc provides the scalars
    Vρ = ∂(ρ E)/∂ρ
    Vσ = ∂(ρ E)/∂σ

Consider a variation dϕi of an orbital (considered real for
simplicity), and let dEtot be the corresponding variation of the
energy. Then the potential Vxc is defined by
    dEtot = 2 ∫ Vxc ϕi dϕi

    dρ = 2 ϕi dϕi
    dσ = 2 ∇ρ ⋅ ∇dρ = 4 ∇ρ ⋅ ∇(ϕi dϕi)
    dEtot = ∫ Vρ dρ + Vσ dσ
          = 2 ∫ Vρ ϕi dϕi + 4 ∫ Vσ ∇ρ ⋅ ∇(ϕi dϕi)
          = 2 ∫ Vρ ϕi dϕi - 4 ∫ div(Vσ ∇ρ) ϕi dϕi
where we performed an integration by parts in the last equation
(boundary terms drop by periodicity). Therefore,
    Vxc = Vρ - 2 div(Vσ ∇ρ)
See eg Richard Martin, Electronic stucture, p. 158

A more algebraic way of deriving these for GGA, which helps see why
this is correct at the discrete level, is the following. Let

B the (Fourier, spherical) space of orbitals
C the (Fourier, cube) space of Fourier densities/potentials
C* the (real-space, cube) space of real-space densities/potentials

X : B -> C be the extension operator
R : C -> B the restriction (adjoint of X)
IF : C -> C* be the IFFT
F : C* -> be the FFT, adjoint of IF
K : C -> C be the gradient operator (multiplication by iK, assuming dimension 1 to
simplify notations)

Then
    ρ = IF(Xϕ).^2
    ρf = F ρ
    ∇ρf = K ρf
    ∇ρ = IF ∇ρf

    dρ  = 2 IF(Xϕ) .* IF(X dϕ)
    d∇ρ = IF K F dρ
        = 2 IF K F (IF(Xϕ) .* IF(X dϕ))

    Etot  = sum(ρ . * E(ρ, ∇ρ.^2))
    dEtot = sum(Vρ .* dρ + 2 (Vσ .* ∇ρ) .* d∇ρ)
          = sum(Vρ .* dρ + 2 (Vσ .* ∇ρ) .* [(IF K F) dρ)]
          = sum(Vρ .* dρ - 2 [(IF K F) (Vσ .* ∇ρ)] .* dρ)
where we used that (IF K F) is anti-self-adjoint and thought of sum(A .* op B )
like ⟨A | op B⟩; this is the discrete integration by parts. Now note that

    sum(V .* dρ) = 2 sum(V .* IF(Xϕ) .* IF(X dϕ))
                 = 2 sum(R(F[ V .* IF(Xϕ) ]) .* dϕ)

where we took the adjoint (IF X)^† = (R F) and the result follows.
=#

#=  Spin-polarised calculations

These expressions can be generalised for spin-polarised calculations.
In this case for example the energy per unit particle becomes
E(ρ_α, ρ_β, σ_αα, σ_αβ, σ_βα, σ_ββ), where σ_ij = ∇ρ_i ⋅ ∇ρ_j
and the XC potential is analogously
    Vxc_s = Vρ_s - 2 ∑_t div(Vσ_{st} ∇ρ_t)
where s, t ∈ {α, β} are the spin components and we understand
    Vρ_s     = ∂(ρ E)/∂(ρ_s)
    Vσ_{s,t} = ∂(ρ E)/∂(σ_{s,t})

Now, in contrast to this libxc explicitly uses the symmetry σ_αβ = σ_βα and sets σ
to be a vector of the three independent components only
    σ = [σ_αα, σ_x, σ_ββ]  where     σ_x = (σ_αβ + σ_βα)/2
Accordingly Vσ has the compoments
    [∂(ρ E)/∂σ_αα, ∂(ρ E)/∂σ_x, ∂(ρ E)/∂σ_ββ]
where in particular ∂(ρ E)/∂σ_x = (1/2) ∂(ρ E)/∂σ_αβ = (1/2) ∂(ρ E)/∂σ_βα.
This explains the extra factor (1/2) needed in the GGA term of the XC potential.
=#

function max_required_derivative(functional)
    functional.family == :lda && return 0
    functional.family == :gga && return 1
    error("Functional family $(functional.family) not known.")
end

"""
Translation from a double index (i, j) of two symmetric axes
to the linear index over the non-redundant components used in libxc
"""
function libxc_symmetric_index(i, j)
    i ≤ j || return libxc_symmetric_index(j, i)
    i + div(j * (j - 1), 2)
end

struct LibxcDensity
    basis
    max_derivative::Int
    ρ_real    # density on real-space grid
    ∇ρ_real   # density gradient on real-space grid
    σ_real    # contracted density gradient on real-space grid
end

"""
Compute density in real space and its derivatives starting from ρ
"""
function LibxcDensity(basis, max_derivative::Integer, ρ::RealFourierArray, ρspin=nothing)
    model = basis.model
    @assert model.spin_polarization in (:collinear, :none, :spinless)
    @assert max_derivative in (0, 1)
    model.spin_polarization == :collinear && (@assert !isnothing(ρspin))
    ifft(x) = real_checked(G_to_r(basis, clear_without_conjugate!(x)))

    n_spin    = model.n_spin_components
    σ_real    = nothing
    ∇ρ_real   = nothing
    if model.spin_polarization == :collinear
        @assert n_spin == 2
        ρ_real    = similar(ρ.real,    n_spin, basis.fft_size...)
        ρ_real[1, :, :, :] .= @. (ρ.real + ρspin.real) / 2
        ρ_real[2, :, :, :] .= @. (ρ.real - ρspin.real) / 2

        if max_derivative > 0
            ρ_fourier = similar(ρ.fourier, n_spin, basis.fft_size...)
            ρ_fourier[1, :, :, :] .= @. (ρ.fourier + ρspin.fourier) / 2
            ρ_fourier[2, :, :, :] .= @. (ρ.fourier - ρspin.fourier) / 2
        end
    else
        @assert n_spin == 1
        ρ_real    = reshape(ρ.real,    1, basis.fft_size...)
        ρ_fourier = reshape(ρ.fourier, 1, basis.fft_size...)
    end

    if max_derivative > 0
        n_spin_σ = div((n_spin + 1) * n_spin, 2)
        ∇ρ_real = similar(ρ.real,   n_spin, basis.fft_size..., 3)
        σ_real  = similar(ρ.real, n_spin_σ, basis.fft_size...)

        for α = 1:3
            iGα = [im * G[α] for G in G_vectors_cart(basis)]
            for σ = 1:n_spin
                ∇ρ_real[σ, :, :, :, α] .= ifft(iGα .* @view ρ_fourier[σ, :, :, :])
            end
        end

        σ_real .= 0
        @views for α in 1:3
            # Does the index transformation (σ, τ) => στ like in libxc_symmetric_index
            σ_real[1, :, :, :] .+= ∇ρ_real[1, :, :, :, α] .* ∇ρ_real[1, :, :, :, α]
            if n_spin > 1
                σ_real[2, :, :, :] .+= ∇ρ_real[1, :, :, :, α] .* ∇ρ_real[2, :, :, :, α]
                σ_real[3, :, :, :] .+= ∇ρ_real[2, :, :, :, α] .* ∇ρ_real[2, :, :, :, α]
            end
        end
    end

    LibxcDensity(basis, max_derivative, ρ_real, ∇ρ_real, σ_real)
end

function Libxc.evaluate(xc::Functional, density::LibxcDensity; kwargs...)
    if xc.family == :lda
        evaluate(xc; rho=density.ρ_real, kwargs...)
    elseif xc.family == :gga
        evaluate(xc; rho=density.ρ_real, sigma=density.σ_real, kwargs...)
    else
        error("Not implemented for functional familiy $(xc.family)")
    end
end

"""
Compute divergence of an operand function, which returns the cartesian x,y,z
components in real space when called with the arguments 1 to 3.
The divergence is also returned as a real-space array.
"""
function divergence_real(operand, basis)
    gradsum = sum(1:3) do α
        operand_α = r_to_G(basis, complex(operand(α)))
        del_α = im * [G[α] for G in G_vectors_cart(basis)]
        del_α .* operand_α
    end
    real(G_to_r(basis, gradsum))
end

function apply_kernel_term_gga(terms, density, perturbation)
    error("apply_kernel_term_gga is not yet working")
    # Follows DOI 10.1103/PhysRevLett.107.216402
    #
    # For GGA V = Vρ - 2 ∇⋅(Vσ ∇ρ) = (∂ε/∂ρ) - 2 ∇⋅((∂ε/∂σ) ∇ρ)
    #
    # dV(r) = f(r,r') dρ(r') = (∂V/∂ρ) dρ + (∂V/∂σ) dσ
    #
    # therefore
    # dV(r) = (∂^2ε/∂ρ^2) dρ - 2 ∇⋅[(∂^2ε/∂σ∂ρ) ∇ρ + (∂ε/∂σ) (∂∇ρ/∂ρ)] dρ
    #       + (∂^2ε/∂ρ∂σ) dσ - 2 ∇⋅[(∂^ε/∂σ^2) ∇ρ  + (∂ε/∂σ) (∂∇ρ/∂σ)] dσ
    #
    # Note dσ = 2∇ρ d∇ρ = 2∇ρ ∇dρ, therefore
    #    - 2 ∇⋅((∂ε/∂σ) (∂∇ρ/∂σ)) dσ = - 2 ∇⋅((∂ε/∂σ) ∇dρ)
    # and (because assumed independent variables): (∂∇ρ/∂ρ) = 0.
    #
    # Note that below the LDA term (∂^2ε/∂ρ^2) dρ is ignored (already dealt with)
    ρ   = density.ρ_real
    ∇ρ  = density.∇ρ_real
    dρ  = perturbation.ρ_real
    ∇dρ = perturbation.∇ρ_real
    dσ = 2sum(∇ρ[α] .* ∇dρ[α] for α in 1:3)

    @views begin
        terms.v2rhosigma[1, :, :, :] .* dσ + divergence_real(density.basis) do α
            @. begin
                -2 * terms.v2rhosigma[1, :, :, :] * ∇ρ[α] * dρ
                -2 * terms.v2sigma2[1, :, :, :]   * ∇ρ[α] * dσ
                -2 * terms.vsigma[1, :, :, :]     * ∇dρ[α]
            end
        end
    end
end


text:
using StatsDiscretizations
using Test
using Infinity

@test cdf(Binomial(5,0.5), Interval{:open,:closed}(0, ∞)) == 1.0
@test cdf(Binomial(5,0.7), Interval{:open,:closed}(-∞,0)) ≈ pdf(Binomial(5,0.7), 0)

#pdf(Normal(), NormalSample(IntervalSets.Interval(0.0,Inf),1.0))


text:
module PointsReconstruction

using LinearAlgebra: norm, diagind, I, dot
using Distances: Euclidean, pairwise, PeriodicEuclidean, colwise
# using Plots: plot, scatter, plot!, scatter!
using Statistics: mean, var
using StaticArrays: SVector
using StatsBase: sample, sqL2dist
using FFTW: ifft, ifftshift, fftshift, fftfreq, fft
# using Interpolations: interpolate, Gridded, Linear, extrapolate
using Distributions: IsoNormal, MvNormal, pdf
# using ThreadTools: tmap1

include("generate_processes.jl")
include("summary_characteristics.jl")
# include("plot_functions.jl")
# include("analyze.jl")
include("wavelets.jl")
include("wavelet_sc.jl")

export uniform, true_uniform, net
export circle_cluster, normal_cluster, cross_cluster
export clusters, voronoi

export distance, distances, K, D, T, H

# export plotfdistr, plot_point_process

# export analyze, lossK, lossD, lossH, default_loss_params

export WaveletParams
export lossW, lossW_params

end


text:
function doubleIt_unmemoized(x)
    println("Evaluating...")
    2x
end

empty!(ObjectIdDict())

doubleIt(x)::Any = begin
    if haskey(ObjectIdDict(), (x,))
        (ObjectIdDict())[(x,)]::Core.Inference.return_type(doubleIt_unmemoized, typeof((x,)))
    else
        (ObjectIdDict())[(x,)] = doubleIt_unmemoized(x)
    end
end

max_stars_repo_name: gcalderone/JuliaForAstronomers

max_issues_repo_name: gcalderone/JuliaForAstronomers

max_forks_repo_name: gcalderone/JuliaForAstronomers


text:
using BoardGameGeek
using Test

@testset "BoardGameGeek.jl" begin
    @testset "Get game" begin
        include("test_game.jl")
    end
    @testset "Get user" begin
        include("test_user.jl")
    end
end


text:
import Base64
import Markdown

export Resource, RemoteResource, LocalResource, DownloadButton

"""
    Resource(src::String, mime=mime_from_filename(src)[, html_attributes::Pair...])

A container for a URL-addressed resource that displays correctly in rich IDEs.

# Examples
```
Resource("https://julialang.org/assets/infra/logo.svg")
```

```
Resource("https://interactive-examples.mdn.mozilla.net/media/examples/flower.webm", :width => 100)
```

```
md\"\"\"
This is what a duck sounds like: \$(Resource("https://interactive-examples.mdn.mozilla.net/media/examples/t-rex-roar.mp3"))
md\"\"\"
```

"""
struct Resource
    src::AbstractString
    mime::Union{Nothing,MIME}
    html_attributes::NTuple{N,Pair} where N
end
Resource(src::AbstractString, html_attributes::Pair...) = Resource(src, mime_fromfilename(src), html_attributes)

function Base.show(io::IO, ::MIME"text/html", r::Resource)
    tag = if r.mime ∈ imagemimes
        :img
    elseif r.mime isa MIME"text/javascript"
        :script
    elseif r.mime ∈ audiomimes
        :audio
    elseif r.mime ∈ videomimes
        :video
    else
        :data
    end
    type_attr = r.mime isa MIME ? [:type => string(r.mime)] : []

    Markdown.withtag(() -> (), io, tag, :src => r.src, :controls => "", type_attr..., (k => string(v) for (k, v) in r.html_attributes)...)
end

RemoteResource = Resource

"""
Create a `Resource` for a local file (a base64 encoded data URL is generated).

# WARNING
`LocalResource` **will not work** when you share the script/notebook with someone else, _unless they have those resources at exactly the same location on their file system_.

## Recommended alternatives (images)
1. Go to [imgur.com](https://imgur.com) and drag&drop the image to the page. Right click on the image, and select "Copy image location". You can now use the image like so: `PlutoUI.Resource("https://i.imgur.com/SAzsMMA.jpg")`.
2. If your notebook is part of a git repository, place the image in the repository and use a relative path: `PlutoUI.LocalResource("../images/cat.jpg")`.

# Examples
```
Resource("./cat.jpg")
```

```
Resource("/home/fons/Videos/nijmegen.mp4", :width => 200)
```

```
md\"\"\"
This is what a duck sounds like: \$(Resource("../data/hannes.mp3"))
md\"\"\"
```
"""
function LocalResource(path::AbstractString, html_attributes::Pair...)
    @warn """`LocalResource` **will not work** when you share the script/notebook with someone else, _unless they have those resources at exactly the same location on their file system_.

    ## Recommended alternatives (images)
    1. Go to [imgur.com](https://imgur.com) and drag&drop the image to the page. Right click on the image, and select "Copy image location". You can now use the image like so: `PlutoUI.Resource("https://i.imgur.com/SAzsMMA.jpg")`.
    2. If your notebook is part of a git repository, place the image in the repository and use a relative path: `PlutoUI.LocalResource("../images/cat.jpg")`."""
    mime = mime_fromfilename(path; default="")
    src = join([
		"data:",
		string(mime),
		";base64,",
		Base64.base64encode(read(path))
	])
	Resource(src, mime, html_attributes)
end

###
# DOWNLOAD BUTTON
###

"""
Button to download a Julia object as a file from the browser.

See [`DownloadButton`](@ref) to do the opposite.

# Examples

```julia
DownloadButton("Roses are red,", "novel.txt")
```
```julia
DownloadButton(UInt8[0xff, 0xd8, 0xff, 0xe1, 0x00, 0x69], "raw.dat")
```
```julia
import JSON
DownloadButton(JSON.json(Dict("name" => "merlijn", "can_cook" => true)), "staff.json")
```

"""
struct DownloadButton
    data
    filename::AbstractString
end
DownloadButton(data) = DownloadButton(data, "result")

function Base.show(io::IO, m::MIME"text/html", db::DownloadButton)
	write(io, "<button>Download...</button> ")
	write(io, db.filename)
	write(io, "")
end

###
# MIMES
###

"Attempt to find the MIME pair corresponding to the extension of a filename. Defaults to `text/plain`."
function mime_fromfilename(filename; default=nothing, filename_maxlength=2000)
	if length(filename) > filename_maxlength
		default
    else
		get(mimepairs, '.' * split(split(split(filename, '?')[1], '#')[1], '.')[end], default)
	end
end

const mimepairs = let
	# This bad boy is from: https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
	originals = Dict(".aac" => "audio/aac", ".bin" => "application/octet-stream", ".bmp" => "image/bmp", ".css" => "text/css", ".csv" => "text/csv", ".eot" => "application/vnd.ms-fontobject", ".gz" => "application/gzip", ".gif" => "image/gif", ".htm" => "text/html", ".html" => "text/html", ".ico" => "image/vnd.microsoft.icon", ".jpeg" => "image/jpeg", ".jpg" => "image/jpeg", ".js" => "text/javascript", ".json" => "application/json", ".jsonld" => "application/ld+json", ".mjs" => "text/javascript", ".mp3" => "audio/mpeg", ".mpeg" => "video/mpeg", ".mp4" => "video/mp4", ".oga" => "audio/ogg", ".ogg" => "audio/ogg", ".ogv" => "video/ogg", ".ogx" => "application/ogg", ".opus" => "audio/opus", ".otf" => "font/otf", ".png" => "image/png", ".pdf" => "application/pdf", ".rtf" => "application/rtf", ".sh" => "application/x-sh", ".svg" => "image/svg+xml", ".tar" => "application/x-tar", ".tif" => "image/tiff", ".tiff" => "image/tiff", ".ttf" => "font/ttf", ".txt" => "text/plain", ".wav" => "audio/wav", ".weba" => "audio/webm", ".webm" => "video/webm", ".webp" => "image/webp", ".woff" => "font/woff", ".woff2" => "font/woff2", ".xhtml" => "application/xhtml+xml", ".xml" => "application/xml", ".xul" => "application/vnd.mozilla.xul+xml", ".zip" => "application/zip")
	Dict((k => MIME(v)) for (k, v) in originals)
end

const imagemimes = [MIME"image/svg+xml"(), MIME"image/png"(), MIME"image/webp"(), MIME"image/tiff"(), MIME"image/jpg"(), MIME"image/jpeg"(), MIME"image/bmp"(), MIME"image/gif"()]
const audiomimes = [MIME"audio/mpeg"(), MIME"audio/wav"(), MIME"audio/aac"(), MIME"audio/ogg"(), MIME"audio/opus"(), MIME"audio/webm"()]
const videomimes = [MIME"video/mpeg"(), MIME"video/ogg"(), MIME"video/webm"(), MIME"video/mp4"()]


text:
using DSGE
using Base.Test
using DataFrames: DataFrame
using DataArrays: @data

@testset "Miscellaneous data handling functions" begin
    # Previous and next quarter arithmetic
    q = Date(1913,12,13)
    DSGE.prev_quarter()
    DSGE.next_quarter()
    @test DSGE.prev_quarter(q) == Date(1913,09,30)
    @test DSGE.next_quarter(q) == Date(1914,03,31)

    start_date = Date(2000,01,01)
    end_date   = Date(2010,01,01)
    DSGE.get_quarter_ends(start_date,end_date)
    DSGE.subtract_quarters(end_date, start_date)

    df1 = DataFrame(a = [1])
    df2 = DataFrame()
    df1, df2 = DSGE.reconcile_column_names(df1, df2)
    @test names(df1) == names(df2)
    @test size(df2, 1) == 0

    quartertodate("11q3")
    quartertodate("1997q4")
    quartertodate("1985-Q1")
    @test_throws ParseError quartertodate("2005q9")
    @test_throws ParseError quartertodate("12345")

    df = DataFrame(date = ["1913-12-23", "1992-11-14", "2002-01-01", "2014-12-19"],
                   x = @data([1, 2, NA, 4]))
    DSGE.format_dates!(:date, df)
end

#DSGE.na2nan!(df)

# Test hpfilter, ensuring missing data are handled correctly
nanfront = [NaN; NaN]
nanback  = [NaN]
y  = [0.; 1.; 0.; 1.]
yn = [nanfront; y; nanback]
yf, yt   = hpfilter(y, 1600)
yfn, ytn = hpfilter(yn, 1600)

@testset "HP Filter and missing data handling" begin
    @test isequal(nanfront, yfn[1:length(nanfront)])
    @test isequal(nanfront, ytn[1:length(nanfront)])
    @test isequal(nanback, yfn[(end-length(nanback)+1):end])
    @test isequal(nanback, ytn[(end-length(nanback)+1):end])
end

nothing