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

"""
    create_replication_set(regions)
    create_replication_set(regions, params::Dict{String,<:Any})

A replication set replicates and encrypts your data to the provided Regions with the
provided KMS key.

# Arguments
- `regions`: The Regions that Incident Manager replicates your data to. You can have up to
  three Regions in your replication set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
"""
create_replication_set(regions; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createReplicationSet", Dict{String, Any}("regions"=>regions, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_replication_set(regions, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createReplicationSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("regions"=>regions, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_response_plan(incident_template, name)
    create_response_plan(incident_template, name, params::Dict{String,<:Any})

Creates a response plan that automates the initial response to incidents. A response plan
engages contacts, starts chat channel collaboration, and initiates runbooks at the
beginning of an incident.

# Arguments
- `incident_template`: Details used to create an incident when using this response plan.
- `name`: The short format name of the response plan. Can't include spaces.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"actions"`: The actions that the response plan starts at the beginning of an incident.
- `"chatChannel"`: The AWS Chatbot chat channel used for collaboration during an incident.
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
- `"displayName"`: The long format of the response plan name. This field can contain spaces.
- `"engagements"`: The contacts and escalation plans that the response plan engages during
  an incident.
- `"tags"`: A list of tags that you are adding to the response plan.
"""
create_response_plan(incidentTemplate, name; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createResponsePlan", Dict{String, Any}("incidentTemplate"=>incidentTemplate, "name"=>name, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_response_plan(incidentTemplate, name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createResponsePlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("incidentTemplate"=>incidentTemplate, "name"=>name, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    create_timeline_event(client_token, event_data, event_time, event_type, incident_record_arn)
    create_timeline_event(client_token, event_data, event_time, event_type, incident_record_arn, params::Dict{String,<:Any})

Creates a custom timeline event on the incident details page of an incident record.
Timeline events are automatically created by Incident Manager, marking key moment during an
incident. You can create custom timeline events to mark important events that are
automatically detected by Incident Manager.

# Arguments
- `client_token`: A token ensuring that the action is called only once with the specified
  details.
- `event_data`: A short description of the event.
- `event_time`: The time that the event occurred.
- `event_type`: The type of the event. You can create timeline events of type Custom Event.
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident record you are
  adding the event to.

"""
create_timeline_event(clientToken, eventData, eventTime, eventType, incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createTimelineEvent", Dict{String, Any}("clientToken"=>clientToken, "eventData"=>eventData, "eventTime"=>eventTime, "eventType"=>eventType, "incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
create_timeline_event(clientToken, eventData, eventTime, eventType, incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/createTimelineEvent", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>clientToken, "eventData"=>eventData, "eventTime"=>eventTime, "eventType"=>eventType, "incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)

"""
    delete_incident_record(arn)
    delete_incident_record(arn, params::Dict{String,<:Any})

Delete an incident record from Incident Manager.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the incident record you are deleting.

"""
delete_incident_record(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteIncidentRecord", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_incident_record(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteIncidentRecord", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_replication_set(arn)
    delete_replication_set(arn, params::Dict{String,<:Any})

Deletes all Regions in your replication set. Deleting the replication set deletes all
Incident Manager data.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the replication set you're deleting.

"""
delete_replication_set(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteReplicationSet", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_replication_set(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteReplicationSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_resource_policy(policy_id, resource_arn)
    delete_resource_policy(policy_id, resource_arn, params::Dict{String,<:Any})

Deletes the resource policy that AWS Resource Access Manager uses to share your Incident
Manager resource.

# Arguments
- `policy_id`: The ID of the resource policy you're deleting.
- `resource_arn`: The Amazon Resource Name (ARN) of the resource you're deleting the policy
  from.

"""
delete_resource_policy(policyId, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteResourcePolicy", Dict{String, Any}("policyId"=>policyId, "resourceArn"=>resourceArn); aws_config=aws_config)
delete_resource_policy(policyId, resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteResourcePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("policyId"=>policyId, "resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    delete_response_plan(arn)
    delete_response_plan(arn, params::Dict{String,<:Any})

Deletes the specified response plan. Deleting a response plan stops all linked CloudWatch
alarms and EventBridge events from creating an incident with this response plan.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the response plan.

"""
delete_response_plan(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteResponsePlan", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_response_plan(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteResponsePlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_timeline_event(event_id, incident_record_arn)
    delete_timeline_event(event_id, incident_record_arn, params::Dict{String,<:Any})

Deletes a timeline event from an incident.

# Arguments
- `event_id`: The ID of the event you are updating. You can find this by using
  ListTimelineEvents.
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident that the event is
  part of.

"""
delete_timeline_event(eventId, incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteTimelineEvent", Dict{String, Any}("eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
delete_timeline_event(eventId, incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/deleteTimelineEvent", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)

"""
    get_incident_record(arn)
    get_incident_record(arn, params::Dict{String,<:Any})

Returns the details of the specified incident record.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the incident record.

"""
get_incident_record(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getIncidentRecord", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_incident_record(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getIncidentRecord", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_replication_set(arn)
    get_replication_set(arn, params::Dict{String,<:Any})

Retrieve your Incident Manager replication set.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the replication set you want to retrieve.

"""
get_replication_set(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getReplicationSet", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_replication_set(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getReplicationSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_resource_policies(resource_arn)
    get_resource_policies(resource_arn, params::Dict{String,<:Any})

Retrieves the resource policies attached to the specified response plan.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the response plan with the attached
  resource policy.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of resource policies to display per page of results.
- `"nextToken"`: The pagination token to continue to the next page of results.
"""
get_resource_policies(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/getResourcePolicies", Dict{String, Any}("resourceArn"=>resourceArn); aws_config=aws_config)
get_resource_policies(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/getResourcePolicies", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    get_response_plan(arn)
    get_response_plan(arn, params::Dict{String,<:Any})

Retrieves the details of the specified response plan.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the response plan.

"""
get_response_plan(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getResponsePlan", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_response_plan(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getResponsePlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_timeline_event(event_id, incident_record_arn)
    get_timeline_event(event_id, incident_record_arn, params::Dict{String,<:Any})

Retrieves a timeline event based on its ID and incident record.

# Arguments
- `event_id`: The ID of the event. You can get an event's ID when you create it or by using
  ListTimelineEvents.
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident that the timeline
  event is part of.

"""
get_timeline_event(eventId, incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getTimelineEvent", Dict{String, Any}("eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
get_timeline_event(eventId, incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/getTimelineEvent", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)

"""
    list_incident_records()
    list_incident_records(params::Dict{String,<:Any})

Lists all incident records in your account. Use this command to retrieve the Amazon
Resource Name (ARN) of the incident record you want to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filters"`: Filter the list of incident records you are searching through. You can
  filter on the following keys:    creationTime     impact     status     createdBy
- `"maxResults"`: The maximum number of results per page.
- `"nextToken"`: The pagination token to continue to the next page of results.
"""
list_incident_records(; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listIncidentRecords"; aws_config=aws_config)
list_incident_records(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listIncidentRecords", params; aws_config=aws_config)

"""
    list_related_items(incident_record_arn)
    list_related_items(incident_record_arn, params::Dict{String,<:Any})

List all related items for an incident record.

# Arguments
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident record that you are
  listing related items for.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of related items per page.
- `"nextToken"`: The pagination token to continue to the next page of results.
"""
list_related_items(incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listRelatedItems", Dict{String, Any}("incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
list_related_items(incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listRelatedItems", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)

"""
    list_replication_sets()
    list_replication_sets(params::Dict{String,<:Any})

Lists details about the replication set configured in your account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results per page.
- `"nextToken"`: The pagination token to continue to the next page of results.
"""
list_replication_sets(; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listReplicationSets"; aws_config=aws_config)
list_replication_sets(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listReplicationSets", params; aws_config=aws_config)

"""
    list_response_plans()
    list_response_plans(params::Dict{String,<:Any})

Lists all response plans in your account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of response plans per page.
- `"nextToken"`: The pagination token to continue to the next page of results.
"""
list_response_plans(; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listResponsePlans"; aws_config=aws_config)
list_response_plans(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listResponsePlans", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags that are attached to the specified response plan.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the response plan.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    list_timeline_events(incident_record_arn)
    list_timeline_events(incident_record_arn, params::Dict{String,<:Any})

Lists timeline events of the specified incident record.

# Arguments
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident that the event is
  part of.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filters"`: Filters the timeline events based on the provided conditional values. You
  can filter timeline events using the following keys:    eventTime     eventType
- `"maxResults"`: The maximum number of results per page.
- `"nextToken"`: The pagination token to continue to the next page of results.
- `"sortBy"`: Sort by the specified key value pair.
- `"sortOrder"`: Sorts the order of timeline events by the value specified in the sortBy
  field.
"""
list_timeline_events(incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listTimelineEvents", Dict{String, Any}("incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
list_timeline_events(incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/listTimelineEvents", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)

"""
    put_resource_policy(policy, resource_arn)
    put_resource_policy(policy, resource_arn, params::Dict{String,<:Any})

Adds a resource policy to the specified response plan.

# Arguments
- `policy`: Details of the resource policy.
- `resource_arn`: The Amazon Resource Name (ARN) of the response plan you're adding the
  resource policy to.

"""
put_resource_policy(policy, resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/putResourcePolicy", Dict{String, Any}("policy"=>policy, "resourceArn"=>resourceArn); aws_config=aws_config)
put_resource_policy(policy, resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/putResourcePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("policy"=>policy, "resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    start_incident(response_plan_arn)
    start_incident(response_plan_arn, params::Dict{String,<:Any})

Used to start an incident from CloudWatch alarms, EventBridge events, or manually.

# Arguments
- `response_plan_arn`: The Amazon Resource Name (ARN) of the response plan that pre-defines
  summary, chat channels, SNS topics, runbooks, title, and impact of the incident.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
- `"impact"`: Defines the impact to the customers. Providing an impact overwrites the
  impact provided by a response plan.  Possible impacts:     1 - Critical impact, this
  typically relates to full application failure that impacts many to all customers.     2 -
  High impact, partial application failure with impact to many customers.    3 - Medium
  impact, the application is providing reduced service to customers.    4 - Low impact,
  customer might aren't impacted by the problem yet.    5 - No impact, customers aren't
  currently impacted but urgent action is needed to avoid impact.
- `"relatedItems"`: Add related items to the incident for other responders to use. Related
  items are AWS resources, external links, or files uploaded to an S3 bucket.
- `"title"`: Provide a title for the incident. Providing a title overwrites the title
  provided by the response plan.
- `"triggerDetails"`: Details of what created the incident record in Incident Manager.
"""
start_incident(responsePlanArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/startIncident", Dict{String, Any}("responsePlanArn"=>responsePlanArn, "clientToken"=>string(uuid4())); aws_config=aws_config)
start_incident(responsePlanArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/startIncident", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("responsePlanArn"=>responsePlanArn, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds a tag to a response plan.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the response plan you're adding the
  tags to.
- `tags`: A list of tags that you are adding to the response plan.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes a tag from a resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the response plan you're removing a tag
  from.
- `tag_keys`: The name of the tag you're removing from the response plan.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_deletion_protection(arn, deletion_protected)
    update_deletion_protection(arn, deletion_protected, params::Dict{String,<:Any})

Update deletion protection to either allow or deny deletion of the final Region in a
replication set.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the replication set you're updating.
- `deletion_protected`: Details if deletion protection is enabled or disabled in your
  account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
"""
update_deletion_protection(arn, deletionProtected; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateDeletionProtection", Dict{String, Any}("arn"=>arn, "deletionProtected"=>deletionProtected, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_deletion_protection(arn, deletionProtected, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateDeletionProtection", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn, "deletionProtected"=>deletionProtected, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_incident_record(arn)
    update_incident_record(arn, params::Dict{String,<:Any})

Update the details of an incident record. You can use this action to update an incident
record from the defined chat channel. For more information about using actions in chat
channels, see Interacting through chat.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the incident record you are updating.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"chatChannel"`: The AWS Chatbot chat channel for responders to collaborate in.
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
- `"impact"`: Defines the impact to customers and applications. Providing an impact
  overwrites the impact provided by the response plan.  Possible impacts:     1 - Critical
  impact, this typically relates to full application failure that impacts many to all
  customers.     2 - High impact, partial application failure with impact to many customers.
    3 - Medium impact, the application is providing reduced service to customers.    4 - Low
  impact, customer might aren't impacted by the problem yet.    5 - No impact, customers
  aren't currently impacted but urgent action is needed to avoid impact.
- `"notificationTargets"`: The SNS targets that AWS Chatbot uses to notify the chat channel
  of updates to an incident. You can also make updates to the incident through the chat
  channel using the SNS topics.  Using multiple SNS topics creates redundancy in the case
  that a Region is down during the incident.
- `"status"`: The status of the incident. An incident can be Open or Resolved.
- `"summary"`: The summary describes what has happened during the incident.
- `"title"`: The title of the incident is a brief and easily recognizable.
"""
update_incident_record(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateIncidentRecord", Dict{String, Any}("arn"=>arn, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_incident_record(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateIncidentRecord", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_related_items(incident_record_arn, related_items_update)
    update_related_items(incident_record_arn, related_items_update, params::Dict{String,<:Any})

Add or remove related items from the related items tab of an incident record.

# Arguments
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident record you are
  updating related items in.
- `related_items_update`: Details about the item you are adding or deleting.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
"""
update_related_items(incidentRecordArn, relatedItemsUpdate; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateRelatedItems", Dict{String, Any}("incidentRecordArn"=>incidentRecordArn, "relatedItemsUpdate"=>relatedItemsUpdate, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_related_items(incidentRecordArn, relatedItemsUpdate, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateRelatedItems", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("incidentRecordArn"=>incidentRecordArn, "relatedItemsUpdate"=>relatedItemsUpdate, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_replication_set(actions, arn)
    update_replication_set(actions, arn, params::Dict{String,<:Any})

Add or delete Regions from your replication set.

# Arguments
- `actions`: An action to add or delete a Region.
- `arn`: The Amazon Resource Name (ARN) of the replication set you're updating.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
"""
update_replication_set(actions, arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateReplicationSet", Dict{String, Any}("actions"=>actions, "arn"=>arn, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_replication_set(actions, arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateReplicationSet", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("actions"=>actions, "arn"=>arn, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_response_plan(arn)
    update_response_plan(arn, params::Dict{String,<:Any})

Updates the specified response plan.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the response plan.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"actions"`: The actions that this response plan takes at the beginning of an incident.
- `"chatChannel"`: The AWS Chatbot chat channel used for collaboration during an incident.
- `"clientToken"`: A token ensuring that the action is called only once with the specified
  details.
- `"displayName"`: The long format name of the response plan. Can't contain spaces.
- `"engagements"`: The contacts and escalation plans that Incident Manager engages at the
  start of the incident.
- `"incidentTemplateDedupeString"`: Used to create only one incident record for an incident.
- `"incidentTemplateImpact"`: Defines the impact to the customers. Providing an impact
  overwrites the impact provided by a response plan.  Possible impacts:     5 - Severe impact
     4 - High impact    3 - Medium impact    2 - Low impact    1 - No impact
- `"incidentTemplateNotificationTargets"`: The SNS targets that AWS Chatbot uses to notify
  the chat channels and perform actions on the incident record.
- `"incidentTemplateSummary"`: A brief summary of the incident. This typically contains
  what has happened, what's currently happening, and next steps.
- `"incidentTemplateTitle"`: The short format name of the incident. Can't contain spaces.
"""
update_response_plan(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateResponsePlan", Dict{String, Any}("arn"=>arn, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_response_plan(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateResponsePlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    update_timeline_event(client_token, event_id, incident_record_arn)
    update_timeline_event(client_token, event_id, incident_record_arn, params::Dict{String,<:Any})

Updates a timeline event. You can update events of type Custom Event.

# Arguments
- `client_token`: A token ensuring that the action is called only once with the specified
  details.
- `event_id`: The ID of the event you are updating. You can find this by using
  ListTimelineEvents.
- `incident_record_arn`: The Amazon Resource Name (ARN) of the incident that the timeline
  event is part of.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"eventData"`: A short description of the event.
- `"eventTime"`: The time that the event occurred.
- `"eventType"`: The type of the event. You can update events of type Custom Event.
"""
update_timeline_event(clientToken, eventId, incidentRecordArn; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateTimelineEvent", Dict{String, Any}("clientToken"=>clientToken, "eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn); aws_config=aws_config)
update_timeline_event(clientToken, eventId, incidentRecordArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ssm_incidents("POST", "/updateTimelineEvent", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>clientToken, "eventId"=>eventId, "incidentRecordArn"=>incidentRecordArn), params)); aws_config=aws_config)


text:
# hyperbolic_integral_expansion.jl

"""
Note: This is meant to be the full expansion for the hyperbolic integral using the correct arbitrary precision methods, but it currently does not converge to the correct values, so made a typo/error somewhere. It's not massively efficient either, so optimisation will be important!
"""

include("arb_hypgeom.jl")
import .one_F_two_fast

include("hypergeometric_expansion.jl")
import .hypergeom_exp

using ArbNumerics
using QuadGK

"""
Solving the hyperbolic integral with expansions and summations.
"""
function arb_binomial(x, y)
    x = ArbReal("$x")
    y = ArbReal("$y")
    one = ArbReal("1")
    ArbNumerics.gamma(x + one) / (ArbNumerics.gamma(y + one) * ArbNumerics.gamma(x - y + one))
end

function hyperbolic_expansion(Ω, β, α, v, w; prec = 64)

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    Ω = ArbReal("$Ω")
    β = ArbReal("$β")
    α = ArbReal("$α")
    v = ArbReal("$v")
    w = ArbReal("$w")

    R = ArbReal((v^2 - w^2) / (w^2 * v))
    a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
    b = ArbReal(R * β / sinh(β * v / 2))

    coefficient = ArbReal(α * β^(3/2) * v^3 / (3 * a * √π * w^3 * sinh(β / 2)))

    M_c(z, n) = ArbReal(hypergeometric_expansion.hypergeom_exp(β * z / 2, n, β, a, 0; prec = prec))
    M_s(z, n) = ArbReal(hypergeometric_expansion.hypergeom_exp(β * z / 2, n, β, a, 1; prec = prec))

    c = ArbReal(cosh(Ω * β / 2))
    s = ArbReal(sinh(Ω * β / 2))

    n = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(ArbReal(0, bits = prec))  # Machine accuracy of specified precision prec.

    while true

        previous_result = ArbReal("$result")

        if mod(n, 2) == 0
            even_term = ArbReal(arb_binomial(n, n/2) * M_c(1, n))
            @fastmath @inbounds @simd for z in [ArbReal(Ω + 1), ArbReal(Ω - 1)]
                even_term += ArbReal(arb_binomial(n, n/2) * (s * M_s(z, n) - c * M_c(z, n)) / 2)
            end
        else
            even_term = ArbReal("0.0")
        end

        k_term = ArbReal("0.0")
        @fastmath @inbounds @simd for k in ArbReal("0"):ArbReal(floor(n / 2 - 1 / 2))
            @fastmath @inbounds @simd for z in [ArbReal(1 + v * (n - 2 * k)), ArbReal(1 - v * (n - 2 * k))]
                k_term += ArbReal(arb_binomial(n, k) * M_c(z, n))
            end
            @fastmath @inbounds @simd for z in [ArbReal(Ω + 1 + v * (n - 2 * k)), ArbReal(Ω - 1 + v * (n - 2 * k)), ArbReal(Ω + 1 - v * (n - 2 * k)), ArbReal(Ω - 1 - v * (n - 2 * k))]
                k_term += ArbReal(arb_binomial(n, k) * (s * M_s(z, n) - c * M_c(z, n)) / 2)
            end
        end

        term = ArbReal(arb_binomial(-3/2, n) * (-b / (2 * a))^n * (even_term
         + k_term))
        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end
        result += ArbReal(term)
        println("term: n = ", n, "\nterm value: ", coefficient * term, "\ncurrent result: ", coefficient * result, "\n")
        n += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) * 2)

            println("Not precise enough. Required error < ", err, ". Increasing precision to ", precision(result) * 2, " bits.\n")

            Ω = ArbReal("$Ω")
            β = ArbReal("$β")
            α = ArbReal("$α")
            v = ArbReal("$v")
            w = ArbReal("$w")

            R = ArbReal((v^2 - w^2) / (w^2 * v))
            a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
            b = ArbReal(R * β / sinh(β * v / 2))

            coefficient = ArbReal(α * β^(3/2) * v^3 / (3 * a * √π * w^3 * sinh(β / 2)))

            M_c(z, n) = ArbReal(hypergeometric_expansion.hypergeom_exp(β * z / 2, n, β, a, 0; prec = precision(result)))
            M_s(z, n) = ArbReal(hypergeometric_expansion.hypergeom_exp(β * z / 2, n, β, a, 1; prec = precision(result)))

            c = ArbReal(cosh(Ω * β / 2))
            s = ArbReal(sinh(Ω * β / 2))

            n = ArbReal("$(n - 1)")
            result = ArbReal("$previous_result")
        end
    end
    result *= ArbReal(coefficient)
    println("Frequency: ", ArbReal(Ω, bits = prec + 8), ". Final result: ", ArbReal(result, bits = prec + 8))
    ArbReal(result, bits = prec + 8)
end

"""
Solving the hyperbolic integral with quadgk
"""
function hyperbolic_integral(Ω, β, α, v, w)

    # Initialise constants.
    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    coefficient = 2 * α * β^(3 / 2) * v^3 / (3 * sqrt(π) * sinh(β / 2) * w^3)

    integrand(x) = (1 - cosh(Ω * x)) * cosh(x - β / 2) / (a^2 - β^2 / 4 + x * (β - x) - b * cosh(v * (x - β / 2)))^(3 / 2)

    integral = quadgk(x -> integrand(x), 0, β / 2)
    println(Ω, " ", coefficient * integral[1])
    return coefficient * integral[1]
end

"""
Plot results for comparison.
"""

Ω_range = 7.01
β = 2.0
α = 7.0
v = 5.8
w = 1.6
R_int = [hyperbolic_integral(Ω, β, α, v, w) for Ω in Ω_range]
# R_osc = [oscillatory_integral_expansion(Ω, β, α, v, w; prec = 32) for Ω in Ω_range]
R_exp = [hyperbolic_expansion(Ω, β, α, v, w; prec = 32) for Ω in Ω_range]
# @show(R_exp)
# p = plot(Ω_range, abs.(R_exp), yaxis = :log)
# plot!(Ω_range, abs.(R_int))
# display(p)

ext: jl


text:
"""
    center_points_simulation(dcm_array, output, header, tmp_center, CCI_slice)

Function ...
"""
function center_points_simulation(dcm_array, output, header, tmp_center, CCI_slice)
	PixelSpacing = PhantomSegmentation.get_pixel_size(header)
	rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])
    sizes = []
    for row in eachrow(output[3])
		area = row[:area]
		append!(sizes, area)
	end

	centroids = output[4]
    largest = Dict()
    for index in 1:length(centroids)
		x = centroids[index][1]
		y = centroids[index][2]
		dist_loc = sqrt((tmp_center[2] - x)^2 + (tmp_center[1] - y)^2)
        dist_loc *= PixelSpacing[1]
        if dist_loc > 31
            largest[index] = [round(y), round(x)]
        else
            nothing
		end
	end

    max_dict = Dict()
	radius = round(2.5 / PixelSpacing[1], RoundUp)
    for key in largest
        tmp_arr = create_circular_mask(rows, cols, (key[2][2], key[2][1]), radius)
        tmp_arr = @. abs(tmp_arr * dcm_array[:,:,CCI_slice]) + abs(tmp_arr * dcm_array[:,:,CCI_slice - 1]) + abs(tmp_arr * dcm_array[:,:,CCI_slice + 1])
        tmp_arr = @. ifelse(tmp_arr == 0, missing, tmp_arr)
        max_dict[key[1]] = median(skipmissing(tmp_arr))
	end
    _, large1_key = maximum(zip(values(max_dict), keys(max_dict)))
    center1 = largest[large1_key]

	center = vec(tmp_center')
	p1 = vec(center1')
	offset = -2π/180
	offset2 = -5π/180
	p2, p3 = find_triangle_points(p1, center; offset=offset, offset2=offset2)
	cent1, cent2, cent3 = vec(p1'), vec(p2'), vec(p3')

    center = find_circle(cent1, cent2, cent3)
	return center, cent1, cent2, cent3
end

"""
    calc_centers_simulation(dcm_array, output, header, tmp_center, CCI_slice)

Function ...
"""
function calc_centers_simulation(dcm_array, output, header, tmp_center, CCI_slice)
    PixelSpacing = PhantomSegmentation.get_pixel_size(header)
    center, center1, center2, center3 = center_points(
        dcm_array, output, header, tmp_center, CCI_slice
    )
    centers = Dict()
    for center_index in (center1, center2, center3)
        side_x = abs(center[1] - center_index[1])
        side_y = abs(center[2] - center_index[2])
        angle = angle_calc(side_x, side_y)
        if (center_index[1] < center[1] && center_index[2] < center[2])
            medium_calc = [
                center_index[1] + (10.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (10.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] + (17 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (17 / PixelSpacing[2]) * cos(angle)),
            ]

        elseif (center_index[1] < center[1] && center_index[2] > center[2])
            medium_calc = [
                center_index[1] + (10.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (10.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] + (17 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (17 / PixelSpacing[2]) * cos(angle)),
            ]

        elseif (center_index[1] > center[1] && center_index[2] < center[2])
            medium_calc = [
                center_index[1] - (10.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (10.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] - (17 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (17 / PixelSpacing[2]) * cos(angle)),
            ]

        elseif (center_index[1] > center[1] && center_index[2] > center[2])
            medium_calc = [
                center_index[1] - (10.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (10.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] - (17 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (17 / PixelSpacing[2]) * cos(angle)),
            ]

        elseif (side_x == 0 && center_index[2] < center[2])
            medium_calc = [center_index[1], center_index[2] + (10.5 / PixelSpacing[2])]
            low_calc = [center_index[1], center_index[2] + (17 / PixelSpacing[2])]

        elseif (side_x == 0 && center_index[2] > center[2])
            medium_calc = [center_index[1], center_index[2] - (10.5 / PixelSpacing[2])]
            low_calc = [center_index[1], center_index[2] - (17 / PixelSpacing[2])]

        elseif (center_index[1] > center[1] && side_y == 0)
            medium_calc = [center_index[1] - (10.5 / PixelSpacing[1]), center_index[2]]
            low_calc = [center_index[1] - (17 / PixelSpacing[1]), center_index[2]]

        elseif (center_index[1] > center[1] && side_y == 0)
            medium_calc = [center_index[1] + (10.5 / PixelSpacing[1]), center_index[2]]
            low_calc = [(center_index[1] + (17 / PixelSpacing[1])), center_index[1]]

        else
            error("unknown angle")
        end

        if center_index == center1
            centers[:Large_HD] = Int.(round.(center_index))
            centers[:Medium_HD] = Int.(round.(medium_calc))
            centers[:Small_HD] = Int.(round.(low_calc))

        elseif center_index == center2
            centers[:Large_MD] = Int.(round.(center_index))
            centers[:Medium_MD] = Int.(round.(medium_calc))
            centers[:Small_MD] = Int.(round.(low_calc))

        elseif center_index == center3
            centers[:Large_LD] = Int.(round.(center_index))
            centers[:Medium_LD] = Int.(round.(medium_calc))
            centers[:Small_LD] = Int.(round.(low_calc))

        else
            nothing
        end
    end
    return centers
end

"""
    mask_inserts_simulation(
        dcm_array,
        masked_array,
        header,
        CCI_slice,
        center_insert;
        calcium_threshold=130,
        comp_connect=trues(3, 3),
    )

Function ...
"""
function mask_inserts_simulation(
    dcm_array,
    masked_array,
    header,
    CCI_slice,
    center_insert;
    calcium_threshold=130,
    comp_connect=trues(3, 3),
)
    output = calc_output(masked_array, header, CCI_slice, calcium_threshold, comp_connect)
    insert_centers = calc_centers_simulation(
        dcm_array, output, header, center_insert, CCI_slice
    )

    PixelSpacing = PhantomSegmentation.get_pixel_size(header)
    rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])

    mask_L_HD = create_circular_mask(
        cols, rows, insert_centers[:Large_HD], (round(5 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_L_MD = create_circular_mask(
        cols, rows, insert_centers[:Large_MD], (round(5 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_L_LD = create_circular_mask(
        cols, rows, insert_centers[:Large_LD], (round(5 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_M_HD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_HD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_M_MD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_MD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_M_LD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_LD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_S_HD = create_circular_mask(
        cols, rows, insert_centers[:Small_HD], (round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_S_MD = create_circular_mask(
        cols, rows, insert_centers[:Small_MD], (round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_S_LD = create_circular_mask(
        cols, rows, insert_centers[:Small_LD], (round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )

    return transpose(mask_L_HD),
    transpose(mask_M_HD),
    transpose(mask_S_HD),
    transpose(mask_L_MD),
    transpose(mask_M_MD),
    transpose(mask_S_MD),
    transpose(mask_L_LD),
    transpose(mask_M_LD),
    transpose(mask_S_LD)
end


text:
###############################################################################
#
# Initialize SPACMono structure
#
###############################################################################
"""
    initialize_spac_canopy!(node::SPACMono{FT}) where {FT<:AbstractFloat}

Initialize the RT parameters for a given
- `node` [`SPACMono`](@ref) type struct
"""
function initialize_spac_canopy!(
            node::SPACMono{FT}
) where {FT<:AbstractFloat}
    # 0.1 create variables required
    @unpack angles, envirs, in_rad, leaves_rt, n_canopy, plant_ps, photo_set,
            rt_con, soil_opt, wl_set = node;
    canopy_rt = node.canopy_rt;
    can_opt = node.can_opt;
    can_rad = node.can_rad;
    plant_hs = node.plant_hs;
    fraction_sl::Array{FT,1} = repeat(canopy_rt.lidf, outer=[canopy_rt.nAzi]) /
                               length(canopy_rt.lazitab);
    n_sl = length(canopy_rt.lidf) * length(canopy_rt.lazitab);

    # fluspect the canopy layers
    for ican in 1:n_canopy
        fluspect!(leaves_rt[ican], wl_set);
    end

    # Four Different steps to compute Short-Wave RT
    canopy_geometry!(canopy_rt, angles, can_opt, rt_con)
    canopy_matrices!(leaves_rt, can_opt);
    short_wave!(canopy_rt, can_opt, can_rad, in_rad,
                soil_opt, rt_con);
    canopy_fluxes!(canopy_rt, can_opt, can_rad, in_rad,
                   soil_opt, leaves_rt, wl_set, rt_con);

    # Compute Long Wave (Last term is LW incoming in W m^-2)
    thermal_fluxes!(leaves_rt, can_opt, can_rad, canopy_rt,
                    soil_opt, [FT(400.0)], wl_set);

    # update the canopy leaf area partition information
    for i_can in 1:n_canopy
        rt_layer = n_canopy + 1 - i_can;
        iPS      = plant_ps[i_can];

        # calculate the fraction of sunlit and shaded leaves
        f_view = (can_opt.Ps[rt_layer] + can_opt.Ps[rt_layer+1]) / 2;

        for iLF in 1:n_sl
            iPS.APAR[iLF] = can_rad.absPAR_sunCab[(rt_layer-1)*n_sl + iLF] *
                            FT(1e6);
            iPS.LAIx[iLF] = f_view * fraction_sl[iLF];
        end
        iPS.APAR[end] = can_rad.absPAR_shadeCab[rt_layer] * FT(1e6);
        iPS.LAIx[end] = 1 - f_view;

        update_leaf_TP!(photo_set, iPS, plant_hs.leaves[i_can], envirs[i_can]);

        envirs[i_can].t_air = T_25(FT);
        envirs[i_can].p_sat = saturation_vapor_pressure(T_25(FT));
        envirs[i_can].p_H₂O = envirs[i_can].p_sat / 2;
    end

    return nothing
end

###############################################################################
#
# As the land model is modularized, there could be multiple places that conatin
#     exact the same information. Thus, these information needs to be synced.
#     This type of information includes leaf area, hydraulic traits, and
#     photosynthetic traits.
#
# The following function is meant for leaf area
#
###############################################################################
"""
    update_LAI!(node::SPACMono{FT}, lai::FT) where {FT<:AbstractFloat}

Update SPAC leaf area index
"""
function update_LAI!(node::SPACMono{FT}, lai::FT) where {FT<:AbstractFloat}
    # 1. Update the general information of the SPAC
    node.la = node.ga * lai;

    # 2. Update the LAI for canopy radiation model
    node.canopy_rt.LAI = lai;

    # 3. Update the LA and LAI for stomatal and photosynthesis models
    for _iPS in node.plant_ps
        _iPS.LA   = node.la / node.n_canopy;
        _iPS.LAI  = lai / node.n_canopy;
        _iPS.tLAI = lai;
    end

    # 4. Update the LA for hydraulic model
    #    TODO abstract this to avoid memory allocation
    for _iHS in node.plant_hs.leaves
        _iHS.area = node.la / node.n_canopy;
    end

    return nothing
end

"""
    update_Kmax!(node::SPACMono{FT}, kmax::FT) where {FT<:AbstractFloat}

Update the maximal hydraulic conductance for SPAC

TODO need to abstractize this to PlantHydraulics.jl
"""
function update_Kmax!(node::SPACMono{FT}, kmax::FT) where {FT<:AbstractFloat}
    # Root:Stem:Leaf conductance ratio is set as 1:2:2.
    #    For trees, roots: 2X, trunk: 8X, branch: 8X, leaves: 4X
    #    For palms, roots: 2X, trunk: 4X, leaves: 4X
    #    For grasses, roots: 2X, leaves: 2X

    # 1. update the hydraulic conductance in Roots
    for root in node.plant_hs.roots
        root.k_max      = 2 * kmax / node.plant_hs.n_root;
        root.k_element .= root.k_max * root.N;
    end

    # 2. If is a tree
    if typeof(node.plant_hs) <: TreeLikeOrganism
        node.plant_hs.trunk.k_max      = 8 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for stem in node.plant_hs.branch
            stem.k_max      = 8 * kmax / node.plant_hs.n_canopy;
            stem.k_element .= stem.k_max * stem.N;
        end
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 3. If is a palm
    if typeof(node.plant_hs) <: PalmLikeOrganism
        node.plant_hs.trunk.k_max      = 4 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 4. If is a grass
    if typeof(node.plant_hs) <: GrassLikeOrganism
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 2 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    return nothing
end

"""
    update_VJRWW!(node::SPACMono{FT}, vcmax::FT) where {FT<:AbstractFloat}

Update Vcmax25(WW), Jmax25(WW), and Rd25(WW) from a given Vcmax25
"""
function update_VJRWW!(node::SPACMono{FT}, vcmax::FT) where {FT<:AbstractFloat}
    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)
    for _iPS in node.plant_ps
        _iPS.ps.Vcmax     = vcmax;
        _iPS.ps.Vcmax25   = vcmax;
        _iPS.ps.Vcmax25WW = vcmax;
        _iPS.ps.Jmax      = vcmax * 1.67;
        _iPS.ps.Jmax25    = vcmax * 1.67;
        _iPS.ps.Jmax25WW  = vcmax * 1.67;
        _iPS.ps.Rd        = vcmax * 0.015;
        _iPS.ps.Rd25      = vcmax * 0.015;
        _iPS.ps.Rd25WW    = vcmax * 0.015;
    end

    return nothing
end

"""
    update_VJR!(node::SPACMono{FT}, ratio::FT) where {FT<:AbstractFloat}

Update effective Vcmax25, Jmax25, and Rd25 from a given tune factor
"""
function update_VJR!(node::SPACMono{FT}, ratio::FT) where {FT<:AbstractFloat}
    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)
    for _iPS in node.plant_ps
        _iPS.ps.Vcmax   = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Vcmax25 = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Jmax    = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Jmax25  = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Rd      = _iPS.ps.Rd25WW * ratio;
        _iPS.ps.Rd25    = _iPS.ps.Rd25WW * ratio;
    end

    return nothing
end

"""
    update_Cab!(node::SPACMono{FT}, cab::FT) where {FT<:AbstractFloat}

Update leaf chlorophyll content, and then rerun `fluspect!`
"""
function update_Cab!(node::SPACMono{FT}, cab::FT) where {FT<:AbstractFloat}
    for _leaf in node.leaves_rt
        _leaf.Cab = cab;
        fluspect!(_leaf, node.wl_set);
    end

    return nothing
end

"""
    update_Weibull!(node::SPACMono{FT}, b::FT, c::FT) where {FT<:AbstractFloat}

Update Weibull B and C for all components
"""
function update_Weibull!(node::SPACMono{FT}, b::FT) where {FT<:AbstractFloat}
    # 1. update B and C for roots
    for _root in node.plant_hs.roots
        _root.vc.b = b;
    end

    # 2. update B and C for trunk (if exists)
    if !(typeof(node.plant_hs) <: GrassLikeOrganism)
        node.plant_hs.trunk.vc.b = b;
    end

    # 3. update B and C for branches (if exist)
    if typeof(node.plant_hs) <: TreeLikeOrganism
        for _stem in node.plant_hs.branch
            _stem.vc.b = b;
        end
    end

    # 4. update B and C for leaves
    for _leaf in node.plant_hs.leaves
        _leaf.vc.b = b;
    end

    return nothing
end

function update_Weibull!(node::SPACMono{FT}, b::FT, c::FT) where {FT<:AbstractFloat}
    # 1. update B and C for roots
    for _root in node.plant_hs.roots
        _root.vc.b = b;
        _root.vc.c = c;
    end

    # 2. update B and C for trunk (if exists)
    if !(typeof(node.plant_hs) <: GrassLikeOrganism)
        node.plant_hs.trunk.vc.b = b;
        node.plant_hs.trunk.vc.c = c;
    end

    # 3. update B and C for branches (if exist)
    if typeof(node.plant_hs) <: TreeLikeOrganism
        for _stem in node.plant_hs.branch
            _stem.vc.b = b;
            _stem.vc.c = c;
        end
    end

    # 4. update B and C for leaves
    for _leaf in node.plant_hs.leaves
        _leaf.vc.b = b;
        _leaf.vc.c = c;
    end

    return nothing
end

max_stars_repo_name: climate-machine/Land


text:
# Utility functions for compressed sensing
export CSImage, loadimage, imvec, vecim, randmask, measure, SNR, renormalize
export binvalue, contrast

"Compressed sensing image data"
struct CSImage{T}
    "The measurements of light intensity"
    measurements :: Array{T,1}
    "The patterns that the light was transmitted through, each as a row"
    patterns     :: Array{T,2}
end

"""
Read image data CSV file with `filename` as `(measurements, patterns)`.
"""
function loadimage(filename)
    d = readdlm(filename, ',')
    CSImage(d[:,1], d[:,2:end])
end

"Turn an image into a vector suitable for processing."
imvec(img, conversion=float) = vec(channelview(conversion.(Gray.(img))))

"Turn a processed vector `vec` into a grayscale image with given `size`."
vecim(vec, size = (Int.(sqrt(length(vec))), :)) = Gray.(reshape(vec, size))

"""
Create a random mask array of zero or one entries.

### Examples

```
julia> randmask(2, 3)
2×3 Array{Float64,2}:
 1.0  0.0  1.0
 1.0  0.0  0.0
```
"""
randmask(dims...) = (sign.(randn(dims...)) .+ 1.0) ./ 2.0

"Simulates measurements of `img` through `pattern` or through random masks."
function measure(img, m, patterns = randmask(m, length(img)))
    CSImage(patterns * img, patterns)
end

"""
Calculate the signal-to-noise ratio of `signal` (\$x\$) to `actual` (\$y\$)
using

\$\\textrm{SNR} = -20 \\log\\left(\\frac{\\|x - y\\|}{\\|y\\|}\\right).\$
"""
SNR(signal, actual) = -20.0 * log10(norm(signal - actual) / norm(actual))

"Rescale `data` to the range \$[0, 1]\$, then multiply by `scale` if given."
function renormalize(data, scale=oneunit(eltype(data)))
    (a, b) = extrema(data)
    if a == b
        fill(scale / 2, size(data))
    else
        @. (data - a) * (scale / (b - a))
    end
end

"""
Creates a function that bins the values of an image according to `values`.

The partition `values` is of the form `[a_0, a_1, ..., a_n]`, where
\$0 = a_0 < a_1 < \\cdots < a_n = 1.\$
By default, the value for a given bin is the midpoint of the values at the
endpoints of the bin, but this may be changed by specifying `bins`.

### Examples

```
julia> binvalue([0.0, 0.2, 0.8, 1.0]).([0.01, 0.3, 0.9])
3-element Array{Float64,1}:
 0.1
 0.5
 0.9

julia> binvalue([0.0, 0.5, 1.0], [0.0, 0.8]).([0.2, 0.6])
2-element Array{Float64,1}:
 0.0
 0.8
```
"""
function binvalue(values, bins = (values[2:end] + values[1:end-1]) ./ 2.0)
    v -> if v ≤ values[1]
        bins[1]
    else
        bins[(1:length(values))[v .≤ values][1] - 1]
    end
end

"Change contrast and brightness of `img` by a clamped linear transformation."
contrast(a, b, img) = Gray.(clamp01.(a .* img .+ b))


text:
part1(input) = compute(input, 10)
part2(input) = compute(input, 40)

function compute(input, num_iterations)
    n, state, table1, table2, first = parse_input(input)
    state2 = Vector{Int}(undef, n^2)
    for i = 1:num_iterations
        iterate!(state, state2, table1, table2)
        state, state2 = state2, state
    end
    k = 1
    freq = zeros(Int, n)
    freq[first + 1] = 1
    for i = 1:n, j = 1:n
        freq[j] += state[k]
        k += 1
    end
    a, b = extrema(freq)
    return b - a
end

function parse_input(input)
    char_mapping = fill(0xFF, 256)
    n = 0
    i = 0
    data = input.data
    while true
        i += 1
        c = data[i]
        c == UInt8('\n') && break
        if char_mapping[c] == 0xFF
            char_mapping[c] = n
            n += 1
        end
    end
    initial_state = zeros(Int, n^2)
    for j = 2:i-1
        a = char_mapping[data[j - 1]]
        b = char_mapping[data[j]]
        initial_state[n * a + b + 1] += 1
    end
    i += 2
    table1 = zeros(UInt8, n^2)
    table2 = zeros(UInt8, n^2)
    while i < length(data)
        a = char_mapping[data[i]]
        b = char_mapping[data[i + 1]]
        c = char_mapping[data[i + 6]]
        table1[n * a + b + 1] = n * a + c + 1
        table2[n * a + b + 1] = n * c + b + 1
        i += 8
    end
    first = char_mapping[data[1]]
    return n, initial_state, table1, table2, first
end

function iterate!(state, state2, table1, table2)
    fill!(state2, 0)
    for (i, n) in enumerate(state)
        state2[table1[i]] += n
        state2[table2[i]] += n
    end
end


text:
@defcomp comp1 begin
   par1 = Parameter(index=[time])
   var1 = Variable(index=[time])

   function run_timestep(p,v,d,t)
        v.var1[t] = p.par1[t]
   end
end

@defcomp comp2 begin
   par2 = Parameter(index=[time, country])
   var2 = Variable(index=[time,country])

   function run_timestep(p,v,d,t)
        v.var2[t] = p.par2[t]
   end
end

@defcomp comp3 begin
    par3 = Parameter(index=[time, town])
    var3 = Variable(index=[time,town])

    function run_timestep(p,v,d,t)
         v.var3[t,:] = p.par3[t,:]
    end
end

 m = Model()
set_dimension!(m, :time, 2000:2010)
set_dimension!(m, :country, [:A, :B, :C])
set_dimension!(m, :town, ["a", "b", "c", "d"])
add_comp!(m, comp2)
add_comp!(m, comp3)
set_param!(m, :comp2, :par2, zeros(11,3))
connect_param!(m, :comp3 => :par3, :comp2 => :var2)


text:
"""
`next!(p, [color])` reports that one unit of progress has been
made. Depending on the time interval since the last update, this may
or may not result in a change to the display.
You may optionally change the color of the display. See also `update!`.
"""
function next!(p::ProgressContinue; options...)
    p.counter += 1
    updateProgress!(p; options...)
end

"""
`cancel(p, [msg], [color=:red])` cancels the progress display
before all tasks were completed. Optionally you can specify the
message printed and its color.
See also `finish!`.
"""
function cancel(p::ProgressContinue, msg::AbstractString = "Aborted before all tasks were completed", color = :red; showvalues = Any[], valuecolor = :blue)
    if p.printed
        move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
        printover(p.output, msg, color)
        printvalues!(p, showvalues; color = valuecolor)
        println(p.output)
    end
    return
end

"""
`finish!(p)` indicates that all tasks have been completed.
See also `cancel`.
"""
function finish!(p::Progress; options...)
    update!(p; p.counter; options...)
end


text:
using MichiganAutonomousVehicles
using NLOptControl
using Parameters
using DataFrames

using Base.Test

include("AutonomousControl/hmmwv.jl")
include("VehicleModels/runtests.jl")


text:
### Global variables
global vf_iterator_maxiter = 400
global Ψ_precission = 10^-2
global Ψ_maxiter = 400
global Ψ_punish_degree = 1.0
global DR_punish_degree = 0.5
global GridDef_maxiter = 10^5

### Custom structs

struct GridDef

    VarNames::NTuple{N,Symbol} where N
    Grids::Vector{T} where T #Union{AbstractRange,AbstractVector{<:Real},,AbstractVector{<:Integer}}
    GridDef(VarNames,Grids) = (length(VarNames) != length(Grids)) ? error("variables and grid dimmension do not match") : new(VarNames,Grids)

end

mutable struct Utils

    f::Function
    qb::Union{Vector{Array{Float64}}, Nothing}
    Utils() = new() # This one is for Init phase
    Utils(f) = new(f,nothing)
    Utils(f,qb) = new(f,qb)

end

struct ModelInit{T<:Real}

    Prices::Vector{T}
    GetAP::Function
    GridDef::GridDef
    λ::Matrix{T}
    β::T
    # Statvar::AbstractString
    HAutility::Function
    HAconstraint::Function
    HAconstraint_Alt::Function
    # VFstart::Array{T}
    # VFalter::Vector{T}
    VF_afterprocess::Function
    # Ψᵢₙᵢₜ::Array{T}
    alternate_return::Function
    MCcondition::Function
    InitParams::Function

end

mutable struct Model{T<:Real}

    Prices::Vector{T}
    AddPrices::Vector{T}
    GetAP::Function
    GridDef::GridDef
    λ::Matrix{T}
    β::T
    # Statvar::Char # statvar must start with 1, Integer only
    HAutility::Function
    Utils::Utils
    HAconstraint::Function
    HAconstraint_Alt::Function
    # VFstart::Array{T}
    VF::Matrix{T}
    DR::Matrix{<:Integer}
    VFalter::Union{Vector{T},Matrix{T}}
    VF_afterprocess::Function
    Ψ::Matrix{T}
    # Ψᵢₙᵢₜ::Array{T}
    Λ::Vector{T}
    alternate_return::Function
    MCcondition::Function
    History::Vector{NTuple{N,Float64}} where N
    InitParams::Function

end


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.QUOTE}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    arg = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, arg, start_col)
    next(ps)

    # Construction
    ret = EXPR{Quote}(EXPR[kw, arg, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")

    return ret
end


text:
module SimpleSDMLayers

using ArchGDAL
using Downloads
using RecipesBase
using Colors, ColorBlendModes
using ZipFile
using Requires
using Distances
using Statistics
using GeometryBasics
export Point, Polygon
using PolygonOps
using StatsBase

# Basic types for the package
include(joinpath("lib", "types.jl"))
export SimpleSDMLayer, SimpleSDMResponse, SimpleSDMPredictor

# Main functions to match coordinates
include(joinpath("lib", "coordinateconversion.jl"))

# Implements a series of interfaces (AbstractArray, iteration, and indexing)
include(joinpath("interfaces", "common.jl"))
include(joinpath("interfaces", "iteration.jl"))
include(joinpath("interfaces", "indexing.jl"))
include(joinpath("interfaces", "broadcast.jl"))

# Additional overloads
include(joinpath("lib", "overloads.jl"))

# Raster clipping
include(joinpath("lib", "clip.jl"))

include(joinpath("lib", "generated.jl"))

include(joinpath("lib", "basics.jl"))
export latitudes, longitudes, boundingbox, grid

include(joinpath("datasets", "ascii.jl"))
include(joinpath("datasets", "geotiff.jl"))
export geotiff

include(joinpath("datasets", "types.jl"))
export WorldClim, CHELSA, EarthEnv
export BioClim, LandCover, HabitatHeterogeneity, Elevation, Topography
export CMIP6, SharedSocioeconomicPathway
export CMIP5, RepresentativeConcentrationPathway
for s in instances(CMIP5)
    @eval export $(Symbol(s))
end
for s in instances(CMIP6)
    @eval export $(Symbol(s))
end
for s in instances(RepresentativeConcentrationPathway)
    @eval export $(Symbol(s))
end
for s in instances(SharedSocioeconomicPathway)
    @eval export $(Symbol(s))
end

include(joinpath("datasets", "layernames.jl"))
export layernames

include(joinpath("datasets", "chelsa", "download.jl"))
include(joinpath("datasets", "chelsa", "bioclim.jl"))

include(joinpath("datasets", "worldclim", "download.jl"))
include(joinpath("datasets", "worldclim", "bioclim.jl"))
include(joinpath("datasets", "worldclim", "elevation.jl"))

include(joinpath("datasets", "earthenv", "download.jl"))
include(joinpath("datasets", "earthenv", "landcover.jl"))
include(joinpath("datasets", "earthenv", "habitatheterogeneity.jl"))
include(joinpath("datasets", "earthenv", "topography.jl"))

include(joinpath("pseudoabsences", "main.jl"))
include(joinpath("pseudoabsences", "radius.jl"))
include(joinpath("pseudoabsences", "randomselection.jl"))
include(joinpath("pseudoabsences", "surfacerangeenvelope.jl"))
export WithinRadius, RandomSelection, SurfaceRangeEnvelope

include(joinpath("operations", "coarsen.jl"))
include(joinpath("operations", "sliding.jl"))
include(joinpath("operations", "mask.jl"))
include(joinpath("operations", "rescale.jl"))
include(joinpath("operations", "mosaic.jl"))
export coarsen, slidingwindow, mask, rescale!, rescale, mosaic

include(joinpath("recipes", "recipes.jl"))
export bivariate

# This next bit is about being able to change the path for raster assets
# globally, which avoids duplication this argument across multiple functions.
_layers_assets_path = get(ENV, "SDMLAYERS_PATH", "assets")
isdir(_layers_assets_path) || mkpath(_layers_assets_path)

# Fixes the export of clip when GBIF or others are loaded
export clip

function __init__()
    @require GBIF="ee291a33-5a6c-5552-a3c8-0f29a1181037" begin
        @info "Loading GBIF support for SimpleSDMLayers.jl"
        include("integrations/GBIF.jl")
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        @info "Loading DataFrames support for SimpleSDMLayers.jl"
        include("integrations/DataFrames.jl")
    end

end

end # module


text:
export Root, Renderable
export draw, ColorRect, Sprite
export defaultDrawRect, defaultDrawSprite

draw(e) = nothing

mutable struct Root <: Entity end # used for root of ECS tree

abstract type Renderable <: Entity end

mutable struct ColorRect <: Renderable
  function ColorRect(w, h; color=colorant"white", kw...)
    instantiate!(new(); w=w, h=h, color=color, kw...)
  end
end

defaultDrawRect(r) = TS_VkCmdDrawRect(vulkan_colors(r.color)..., r.abs_pos.x, r.abs_pos.y, r.w, r.h)

draw(r::ColorRect) = defaultDrawRect(r)

mutable struct Sprite <: Renderable
  # allow textures larger than a single sprite
  # by supporting cell_size, cell_ind, and region
  # fields. cell_size turns img into a matrix of
  # adjacent regions of the same size starting in
  # the top left of the texture. cell_ind is a 0-indexed
  # (row,col) index into a particular cell. region
  # overrides both and denotes the top-left corner
  # of a rectangle, its width, and its height, and
  # uses only that region of the texture.
  function Sprite(img; cell_size=[0, 0], region=[0, 0, 0, 0], cell_ind=[0, 0],
    color=colorant"white", scale=XYZ(1,1,1), kw...)
    instantiate!(new(); img=img, cell_size=cell_size,
      region=region, cell_ind=cell_ind, color=color,
      scale=scale, kw...)
  end
end

defaultDrawSprite(s) = TS_VkCmdDrawSprite(s.img, vulkan_colors(s.color)...,
  s.region[1], s.region[2], s.region[3], s.region[4],
  s.cell_size[1], s.cell_size[2], s.cell_ind[1], s.cell_ind[2],
  s.abs_pos.x, s.abs_pos.y, s.scale.x, s.scale.y)

draw(s::Sprite) = defaultDrawSprite(s)


text:
"""
    propagate(F::judiPropagator{T, mode}, q)

Base propagation interfaces that calls the devito `mode` propagator (forward/adjoint/..)
with `q` as a source. The return type is infered from `F`.
"""
function propagate(F::judiPropagator{T, O}, q::AbstractArray{T}) where {T, O}
    srcGeometry, srcData, recGeometry, recData, dm = make_input(F, q)
    return time_modeling(F.model, srcGeometry, srcData, recGeometry, recData, dm, O, F.options)
end

function propagate(F::judiPropagator{T, :adjoint}, q::AbstractArray{T}) where {T}
    srcGeometry, srcData, recGeometry, recData, dm = make_input(F, q)
    return time_modeling(F.model, recGeometry, recData, srcGeometry, srcData, dm, :adjoint, F.options)
end

propagate(t::Tuple{judiPropagator, AbstractArray}) = propagate(t[1], t[2])

"""
    run_and_reduce(func, pool, nsrc, arg_func)

Runs the function `func` for indices `1:nsrc` within arguments `func(arg_func(i))`. If the
the pool is empty, a standard loop and accumulation is ran. If the pool is a julia WorkerPool or
any custom Distributed pool, the loop is distributed via `remotecall` followed by are binary tree remote reduction.
"""
function run_and_reduce(func, pool, nsrc, arg_func::Function)
    res = Vector{_TFuture}(undef, nsrc)
    for i = 1:nsrc
        args_loc = arg_func(i)
        res[i] = remotecall(func, pool, args_loc...)
    end
    res = reduce!(res)
    return res
end

run_and_reduce(func, ::Nothing, nsrc, arg_func::Function) = mapreduce(i -> func(arg_func(i)...), single_reduce!, 1:nsrc)

src_i(::judiAbstractJacobian{T, :born, FT}, q::dmType{T}, ::Integer) where {T<:Number, FT} = q
src_i(::judiPropagator{T, O}, q::judiMultiSourceVector{T}, i::Integer) where {T, O} = q[i]
src_i(::judiPropagator{T, O}, q::Vector{<:Array{T}}, i::Integer) where {T, O} = q[i]

get_nsrc(::judiPropagator, q::judiMultiSourceVector) = q.nsrc
get_nsrc(::judiPropagator, q::Vector{<:Array}) = length(q)
get_nsrc(J::judiAbstractJacobian, ::dmType{T}) where T<:Number = J.q.nsrc

"""
    multi_src_propagate(F::judiPropagator{T, O}, q::AbstractVector)

Propagates the source `q` with the `F` propagator. The return type is infered from `F` and the
propagation kernel is defined by `O` (forward, adjoint, born or adjoint_born).
"""
function multi_src_propagate(F::judiPropagator{T, O}, q::AbstractVector{T}) where {T<:Number, O}
    q = process_input_data(F, q)
    # Number of sources and init result
    nsrc = get_nsrc(F, q)
    pool = _worker_pool()
    arg_func = i -> (F[i], src_i(F, q, i))
    # Distribute source
    res = run_and_reduce(propagate, pool, nsrc, arg_func)
    return as_vec(res, Val(F.options.return_array))
end

"""
    multi_src_fg!(G, model, q, dobs, dm; options=Options(), nlind=false, lin=false)

This is the main multi-source wrapper function for `fwi_objective` and `lsrtm_objective`.
Computes the misifit and gradient (LSRTM if `lin` else FWI) for the given `q` source and `dobs` and
perturbation `dm`.
"""
function multi_src_fg!(G, model, q, dobs, dm; options=Options(), kw...)
    check_non_indexable(kw)
    # Number of sources and init result
    nsrc = try q.nsrc catch; dobs.nsrc end
    pool = _worker_pool()
    # Distribute source
    arg_func = i -> (model, q[i], dobs[i], dm, options[i], values(kw)...)
    # Distribute source
    res = run_and_reduce(multi_src_fg, pool, nsrc, arg_func)
    f, g = as_vec(res, Val(options.return_array))
    G .= g
    return f
end

check_non_indexable(d::Dict) = for (k, v) in d check_non_indexable(v) end
check_non_indexable(x) = false
check_non_indexable(::Bool) = true
check_non_indexable(::Number) = true
check_non_indexable(::judiMultiSourceVector) = throw(ArgumentError("Keyword arguments must not be source dependent"))
function check_non_indexable(::AbstractArray)
    @warn "keyword arguement Array considered source independent and copied to all workers"
    true
end


text:
using Publish
using DLPipelines

p = Publish.Project(DLPipelines)
rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)
deploy(DLPipelines; root = "/DLPipelines.jl", force = true, label = "dev")


text:
module RK

export ide_solve

# ============================================================================================== #
# ============================================================================================== #
function ide_solve(idefun, Core, delays_int, history, tspan, stepsize, delays=false; overlapping=false)
    #      idefun - right-hand side function (t - time, y - solution, z - dicrete delays, i - integral)
    #        Core - Kernel (integrated function)
    #  delays_int - distributed delays function (lower integration limit)
    #     History - history function
    #       tspan - solution interval
    #    stepsize - step of numerical Runge-Kutta method
    #      delays - (optional) discrete delays function (if idefun has 'z')
    # overlapping - (optional) if equation has overlapping in discrete delays. This option uses the 7-step method

    t0 = tspan[1]       # t begin
    tf = tspan[2]       # t end
    y0 = history(t0)    # initial solution
    neq = length(y0)    # number of equations
    htry = stepsize     # constant mesh step-size
    h = htry            # current step-size (for computing the last step)

    if delays != false
        d_t0 = delays(t0,y0)
        nz = length(d_t0)
    else
        nz = 1
    end
    # ============================================================== #
    if overlapping
        # VIDE Runge-Kutta Tavernini
        A = [0 1 3/8  1/2  5/24  1/6
            0 0 1/8  1/2     0    0
            0 0   0    0   1/3 -1/3
            0 0   0    0 -1/24  1/6
            0 0   0    0     0    1
            0 0   0    0     0    0]

        b = [1/6; 0; 0; 0; 2/3; 1/6]
        c = [0 1 1/2 1 1/2 1]
        d = [1/2; 1/2; 1/2; 1/2; 1]
        b4 = b4_Tav
        F_stages_calc = [2 3]
        method = "Tav"
    else
        # VIDE Runge-Kutta 4
        A = [0 1/2   0  0
            0   0 1/2  0
            0   0   0  1
            0   0   0  0 ]

        b = [1/6; 1/3; 1/3; 1/6]
        c = [0 1/2 1/2 1]
        d = [1/2; 1/2; 1]
        b4 = b4_RK4
        F_stages_calc = [2 4]
        method = "RK4"
    end

    s = length(b)
    # ============================================================== #
    nint = length(delays_int(t0))
    # Calculate integral (F) in history
    F = zeros(nint)
    F_1 = zeros(nint)
    F_half = zeros(nint)

    # Piece from delays(t0) to the next mesh point
    # tj, tj_half, tj_1, Core_tj, Core_tj_1, Core_tj_h are always needed
    # to use Simpson's Method for numerical approximations of integrals
    tj = delays_int(t0)     # Begin
    for ij = 1:nint
        step = trunc(Int, (t0 - tj[ij]) / htry)     # The number of memorized intervals of history
        tj_1 = t0 - step * h                        # End
        tj_half = (tj[ij] + tj_1) / 2               # Half segment

        # Calculate Kernel values at the nodes
        Core_tj = Core(t0, tj[ij], history(tj[ij]))
        Core_tj_1 = Core(t0, tj_1, history(tj_1))
        Core_tj_h = Core(t0, tj_half, history(tj_half))

        # Simpson's method
        F[ij] = F[ij] + int_simpson(tj_1 - tj[ij], Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

        # Main integral over the mesh points
        Core_tj = Core_tj_1
        for j = step - 1 : -1 : 0
            tj_1 = t0 - j * h
            tj_half = tj_1 - h / 2

            # Calculate Kernel values at the nodes
            Core_tj_h = Core(t0, tj_half, history(tj_half))
            Core_tj_1 = Core(t0, tj_1, history(tj_1))

            # Simpson's method
            F[ij] = F[ij] + int_simpson(h, Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

            # Kernel of tj_1 is equal to tj of the next step
            Core_tj = Core_tj_1
        end
    end
    # F = ones(neq).*(exp(t0)-exp(delays(t0)))
    # F = ones(neq).*(exp(t0*(1-t0))/(1-t0)-exp(delays(t0)*(1-t0))/(1-t0))
    # ============================================================== #
    # Initialization | First Step | Y | K
    t = t0
    y = zeros(neq, 1)
    y[:, 1] = ones(neq) .* y0
    k = 1   # step number

    z = zeros(neq, nz)
    if delays != false
        for j = 1 : nz
            z[:, j] = ones(neq) .* history(d_t0[j])
        end
    end

    Y = zeros(neq, s)
    K = zeros(neq, s, 1)

    K[:, 1, k] = ones(neq) .* idefun(t[k], y[:, 1], z, F)
    Core_di = zeros(nint, s)
    # ============================================================== #

    while t[k] < tf
        # ============================================================== #
        # Last step
        if t[k] + h > tf
            h = tf - t[k]
        end
        # ============================================================== #
        Z = zeros(nint, s)
        Y[:, 1] = y[:, k]

        # Runge-Kutta steps
        for i = 2 : s
            ti = t[k] + c[i] * h
            # ============================================================== #
            # Calculate integral (F)
            if i in F_stages_calc # c[3] = c[2] so the same F value is used
                F = zeros(nint)

                dtk_begin = delays_int(ti)  # lower integration limit
                for ij = 1 : nint
                    if dtk_begin[ij] < t0
                        # ============================================================== #
                        # Integral begins in the history

                        # Step of delays(ti) in the history
                        step = trunc(Int, (t0 - dtk_begin[ij]) / htry)

                        # Add piece from dtk_begin to the next mesh point in the history
                        tj = dtk_begin[ij]
                        tj_1 = t0 - step * htry
                        tj_half = (tj + tj_1) / 2

                        # Calculate Kernel values at the nodes
                        Core_tj = Core(ti, tj, history(tj))
                        Core_tj_1 = Core(ti, tj_1, history(tj_1))
                        Core_tj_h = Core(ti, tj_half, history(tj_half))

                        # Simpson's method
                        F[ij] = F[ij] + int_simpson(tj_1 - tj, Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

                        # Main integral in the history
                        Core_tj = Core_tj_1
                        for j = step - 1 : -1 : 0
                            tj_1 = t0 - j * htry
                            tj_half = tj_1 - htry / 2

                            # Calculate Kernel values at the nodes
                            Core_tj_1 = Core(ti, tj_1, history(tj_1))
                            Core_tj_h = Core(ti, tj_half, history(tj_half))

                            # Simpson's method
                            F[ij] = F[ij] + int_simpson(htry, Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

                            # Kernel of tj_1 is equal to tj of the next step
                            Core_tj = Core_tj_1
                        end

                        # Add integral in the solution to t(k)
                        for j = 2 : k
                            tj_half = t[j] - htry / 2
                            y_half = y[:, j - 1] + htry * ( K[:, :, j - 1] * b4(1/2) )

                            # Calculate Kernel values at the nodes
                            if neq == 1
                                Core_tj_h = Core(ti, tj_half, y_half[1])
                                Core_tj_1 = Core(ti, t[j], y[1, j])
                            else
                                Core_tj_h = Core(ti, tj_half, y_half)
                                Core_tj_1 = Core(ti, t[j], y[:, j])
                            end

                            # Simpson's method
                            F[ij] = F[ij] + int_simpson(htry, Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

                            # Kernel of tj_1 is equal to tj of the next step
                            Core_tj = Core_tj_1
                        end
                        # ============================================================== #
                    else
                        # ============================================================== #
                        # Integral only over the solution

                        # Step of delays(ti) in the solution
                        step = trunc(Int, (dtk_begin[ij] - t0) / htry + 1)

                        # Add piece from dtk_begin to the mesh point in the solution
                        tj_half = (t[step + 1] + dtk_begin[ij]) / 2

                        y_begin = y[:, step] + htry * ( K[:, :, step] * b4((dtk_begin[ij] - t[step]) / htry) )
                        y_begin_h = y[:, step] + htry * ( K[:, :, step] * b4((tj_half - t[step]) / htry) )

                        # Calculate Kernel values at the nodes
                        if neq == 1
                            Core_tj = Core(ti, dtk_begin[ij], y_begin[1])
                            Core_tj_1 = Core(ti, t[step + 1], y[1, step + 1])
                            Core_tj_h = Core(ti, tj_half, y_begin_h[1])
                        else
                            Core_tj = Core(ti, dtk_begin[ij], y_begin)
                            Core_tj_1 = Core(ti, t[step + 1], y[:, step + 1])
                            Core_tj_h = Core(ti, tj_half, y_begin_h)
                        end

                        # Simpson's method
                        F[ij] = F[ij] + int_simpson(t[step + 1] - dtk_begin[ij], Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

                        # Main integral to t(k)
                        Core_tj = Core_tj_1
                        for j = step + 2 : k
                            tj_half = t[j] - htry / 2
                            y_half = y[:, j - 1] + htry * ( K[:, :, j - 1] * b4(1/2) )

                            # Calculate Kernel values at the nodes
                            if neq == 1
                                Core_tj_h = Core(ti, tj_half, y_half[1])
                                Core_tj_1 = Core(ti, t[j], y[1, j])
                            else
                                Core_tj_h = Core(ti, tj_half, y_half)
                                Core_tj_1 = Core(ti, t[j], y[:, j])
                            end

                            # Simpson's method
                            F[ij] = F[ij] + int_simpson(htry, Core_tj[ij], Core_tj_1[ij], Core_tj_h[ij])

                            # Kernel of tj_1 is equal to tj of the next step
                            Core_tj = Core_tj_1
                        end
                        # ============================================================== #
                    end
                end
                if method == "Tav"
                    if i == 2
                        F_1 = F
                    else
                        F_half = F
                    end
                else
                    if i == 2
                        F_half = F
                    end
                end
            end
            if method == "Tav"
                if i == 4 || i == 6
                    F = F_1
                elseif i == 5
                    F = F_half
                end
            else
                if i == 3
                    F = F_half
                end
            end

            # F = ones(neq).*(exp(t[k])-exp(delays(ti)))
            # F = ones(neq).*(exp(t[k]*(1-ti))/(1-ti)-exp(delays(ti)*(1-ti))/(1-ti))
            # ============================================================== #
            # Y2-S
            Y[:, i] = y[:, k] + h * (K[:, 1:i-1, k] * A[1:i-1, i])

            # Z2-S
            if nint == 1
                Core_di[1, i - 1] = Core(t[k] + d[i - 1] * h, t[k] + c[i - 1] * h, Y[1, i - 1])
            else
                Core_di[:, i - 1] = Core(t[k] + d[i - 1] * h, t[k] + c[i - 1] * h, Y[:, i - 1])
            end
            Z[:, i] = h * (Core_di[:, 1:i-1] * A[1:i-1, i])

            # ============================================================== #
            # Finding delays Z
            if delays != false
                d_ti = delays(ti, Y[:, i])

                for kz = 1 : nz
                    if d_ti[kz] < t0
                        z[:, kz] = ones(neq) .* history(d_ti[kz])
                    elseif ti < d_ti[kz]
                        # wrong overlapping
                        error("Delays went ahead.")
                    elseif t[k] - d_ti[kz] <= 0
                        # overlapping
                        if method == "Tav"
                            teta = (d_ti[kz] - t[k]) / h
                            z[:, kz] = y[:, k] + h * (K[:, 1:i-1, k] * MatrixA(teta, i - 1))
                        else
                            error("Overlapping. Choose method Tavernini with key 'overlapping=true'")
                        end
                    else
                        # find t
                        # ============ Binary search algorithm =========== #
                        tcur = d_ti[kz]
                        iz = 1
                        jz = length(t)
                        nstep = trunc(Int, jz / 2)

                        while (((t[nstep + 1] < tcur || t[nstep] > tcur)) && iz < jz)
                            if tcur > t[nstep]
                                iz = nstep + 1
                            else
                                jz = nstep - 1
                            end
                            nstep = trunc(Int, (iz + jz) / 2)
                        end
                        # ================================================ #

                        # find y_t
                        theta = (tcur - t[nstep]) / htry
                        z[:, kz] = y[:, nstep] + htry * (K[:, :, nstep] * b4(theta))
                    end
                end
            end
            # ============================================================== #
            # K2-S
            K[:, i, k] = ones(neq) .* idefun(ti, Y[:, i], z, F + Z[:, i])
        end
        # ============================================================== #
        # Final approximation of RK Method
        t = [t; t[k] + h]
        y = [y y[:, k] + h * ( K[:, :, k] * b )]
        # ============================================================== #
        # Calculate K(1) for next step
        y_k_half = y[:, k] + h * ( K[:, :, k] * b4(1/2) )
        if neq == 1
            Core_tj = Core(t[k + 1], t[k], y[1, k])
            Core_tk_half = Core(t[k + 1], t[k] + h / 2, y_k_half[1])
            Core_tk = Core(t[k + 1], t[k + 1], y[1, k + 1])
        else
            Core_tj = Core(t[k + 1], t[k], y[:, k])
            Core_tk_half = Core(t[k + 1], t[k] + h / 2, y_k_half)
            Core_tk = Core(t[k + 1], t[k + 1], y[:, k + 1])
        end

        for ij = 1 : nint
            F[ij] = F[ij] + int_simpson(h, Core_tj[ij], Core_tk[ij], Core_tk_half[ij])
        end
        K1 = zeros(neq, s, 1)
        K = cat(K, K1, dims = 3)

        # F = ones(neq).*(exp(t[k+1])-exp(delays(t[k+1])))
        # F = exp(t[k+1]*(1-t[k+1]))/(1-t[k+1])-exp(delays(t[k+1])*(1-t[k+1]))/(1-t[k+1])

        K[:, 1, k + 1] = ones(neq) .* idefun(t[k + 1], y[:, k + 1], z, F)

        k = k + 1
    end
    return t, y'
end
# ============================================================================================== #
# ============================================================================================== #

function int_simpson(h, y_begin, y_end, y_half)
    # Simpson's method
    return h/6 * (y_begin + 4 * y_half + y_end)
end

function MatrixA(a, step)
    A = zeros(step);
    if step == 0
        A = 0
    elseif step == 1
        A = a
    elseif step == 2
        A = [a * (1 - a * 1/2); 1/2 * a^2]
    elseif step == 3
        A = [a * (1 - a * 1/2); 1/2 * a^2; 0]
    elseif step == 4
        sqrA = a^2
        A = [a * (1 + a * (-3/2 + a * 2/3)); 0;
            sqrA * (2 - a * 4/3); sqrA * (-1/2 + a * 2/3)]
    elseif step == 5
        sqrA = a^2
        A = [a * (1 + a * (-3/2 + a * 2/3)); 0;
            sqrA * (1 - a * 4/3); sqrA * (-1/2 + a * 2/3); sqrA]
    end
    return A
end

function b4_Tav(a)
    x = zeros(6,1)
    sqrA = a^2
    x[1,1] = a * (1 + a * (-3/2 + a * 2/3))
    x[2,1] = 0
    x[3,1] = 0
    x[4,1] = 0
    x[5,1] = sqrA * (2 + a * -4/3)
    x[6,1] = sqrA * (-1/2 + a * 2/3)
    return x
end

function b4_RK4(a)
    x = zeros(4,1)
    sqrA = a^2
    x[1,1] = a * (1 + a * (-3/2 + a * 2/3))
    x[2,1] = sqrA * (1 + a * -2/3)
    x[3,1] = sqrA * (1 + a * -2/3)
    x[4,1] = sqrA * (-1/2 + a * 2/3)
    return x
end

end

max_stars_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_issues_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_forks_repo_name: vitesempl/IDE-Delay-Runge-Kutta


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(
    joinpath(pkgpath, "test", "3_PowerSimulations_examples", "01_operations_problems.jl"),
)

sys_RT = build_system(PSITestSystems, "modified_RTS_GMLC_RT_sys")

#set_device_model!(template_ed, GenericBattery, BookKeeping)
template_uc

template_ed = template_economic_dispatch()

problems = SimulationProblems(
    UC = OperationsProblem(template_uc, sys, optimizer = solver),
    ED = OperationsProblem(
        template_ed,
        sys_RT,
        optimizer = solver,
        balance_slack_variables = true,
    ),
)

feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))

feedforward = Dict(
    ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
        binary_source_problem = PSI.ON,
        affected_variables = [PSI.ACTIVE_POWER],
    ),
)

intervals = Dict("UC" => (Hour(24), Consecutive()), "ED" => (Minute(15), Consecutive()))

DA_RT_sequence = SimulationSequence(
    problems = problems,
    intervals = intervals,
    ini_cond_chronology = InterProblemChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)

sim = Simulation(
    name = "rts-test",
    steps = 1,
    problems = problems,
    sequence = DA_RT_sequence,
    simulation_folder = dirname(dirname(pathof(SIIPExamples))),
)

build!(sim)

execute!(sim, enable_progress_bar = false)

results = SimulationResults(sim);
uc_results = get_problem_results(results, "UC"); # UC stage result metadata
ed_results = get_problem_results(results, "ED"); # ED stage result metadata

read_variables(uc_results, names = [:P__ThermalStandard, :P__RenewableDispatch])

read_parameter(
    ed_results,
    :P__max_active_power__RenewableFix,
    initial_time = DateTime("2020-01-01T06:00:00"),
    count = 5,
)

read_realized_variables(uc_results, names = [:P__ThermalStandard, :P__RenewableDispatch])

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module TestSchemalizer

using Test
using TSML

function generatedf()
  Random.seed!(123)
  gdate = DateTime(2015,1,1):Dates.Minute(30):DateTime(2015,1,2) |> collect
  len = length(gdate)
  sex = rand(['m','f'],len)
  x1=rand(1:100,len)
  x2=rand(1:100,len)
  x3=rand(1:1000,len)
  DataFrame(date=gdate,sex=sex,f1=x1,f2=x2,f3=x3)
end

function test_schemalizer()
  df = generatedf()
  m = Schemalizer(Dict(:schema=>
      Dict(:sex => ML.Categorical)))
  fit!(m,df)
  res=transform!(m,df)
  @test isapprox(sum(Matrix(res.^2)),193,atol=1e-2)
  m = Schemalizer(Dict())
  fit!(m,df)
  res = transform!(m,df)
  @test isapprox(sum(Matrix(res.^2)),144,atol=1e-2)
end
@testset "Schemalizer" begin
  test_schemalizer()
end

end


text:
# Global configuration settings

"Configuration options"
module config
    # Data directories
    "Location of the data files in this package"
    const datadir = normpath(joinpath(dirname(@__FILE__), "..", "data"))
end

include("config_defaults.jl")


text:
# Extracting variable values along the magnetopause Bz==0 contour from multiple frames.
#
# Usage:
#   julia -t 4 demo_magnetopause_2d_mt.jl
# or
#   JULIA_NUM_THREADS=4 julia demo_magnetopause_2d_mt.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob
using Vlasiator: RE # Earth radius, [m]
using JLD2: jldsave

"""
    extract_magnetopause_var(files; zmin=-4RE, zmax=4RE, verbose=true)

Extract variables on the magnetopause defined by ``B_z = 0``.
"""
function extract_magnetopause_var(files; zmin=-4RE, zmax=4RE, verbose=true)
   nfiles = length(files)

   verbose && println("Number of files: $nfiles")

   meta = load(files[1])

   x = LinRange{Float32}(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])
   z = LinRange{Float32}(meta.coordmin[3], meta.coordmax[3], meta.ncells[3])

   z_range = let zmin_ = searchsortedfirst(z, zmin), zmax_ = searchsortedlast(z, zmax)
      zmin_:zmax_
   end

   x_crossing = zeros(Float32, length(z_range))
   z_crossing = z[z_range]

   cellids = Vector{UInt64}(undef, length(z_range))

   x_crossings = zeros(length(z_range), nfiles)
   v = zeros(3, length(z_range), nfiles)
   ey = zeros(length(z_range), nfiles)

   Threads.@threads for ifile in eachindex(files)
      verbose && println("$(files[ifile]) on thread $(Threads.threadid())")
      meta = load(files[ifile])

      # Obtain thermal temperature
      b = meta["vg_b_vol"]
      b = reshape(b, 3, meta.ncells[1], meta.ncells[3])

      # Extract the last point from right to left which fulfills Bz < 0
      for (i,k) in enumerate(z_range) # scan in z direction
         ind_ = findlast(>(0), @view b[3,:,k]) + 1 # count from upstream
         isnothing(ind_) && (ind_ = 1) # if not found then set to 1
         x_crossing[i] = x[ind_]
      end

      cellids = [ getcell(meta, [x_crossing[i], 0, z_crossing[i]])
         for i in eachindex(x_crossing, z_crossing) ]

      x_crossings[:,ifile] = x_crossing
      v[:,:,ifile] = readvariable(meta, "proton/vg_v", cellids)
      ey[:,ifile] = readvariable(meta, "vg_e_vol", cellids)[2,:]
   end

   z_crossing, x_crossings, v, ey
end

######

files = glob("bulk*.vlsv", ".")

@time z,x,v,ey = extract_magnetopause_var(files)

jldsave("magnetopause.jld2"; z,x,v,ey)


text:
## Makes the model associated with a particle filter using the "locally optimal
## proposal". Note that this is only optimal in the specific and limited sense
## that the weights are a function of xprev only, and hence the simulated value
## of x is not a source of variability.

mutable struct LGLOPParticle
  x::Float64
  xprev::Float64
  LGLOPParticle() = new()
end

function makeLGLOPModel(theta::LGTheta, ys::Vector{Float64})
  n::Int64 = length(ys)

  RC2v0::Float64 = theta.R + theta.C * theta.C * theta.v0
  v0RoverRC2v0::Float64 = theta.v0*theta.R/RC2v0
  sqrtv0RoverRC2v0::Float64 = sqrt(v0RoverRC2v0)

  mu1::Float64 = v0RoverRC2v0*(theta.x0/theta.v0 + theta.C*ys[1]/theta.R)
  tmp1::Float64 = theta.C * theta.x0 - ys[1]
  lG1::Float64 = -0.5 * log(2 * π * RC2v0) - tmp1 * 0.5/RC2v0 * tmp1

  RC2Q::Float64 = theta.R + theta.C * theta.C * theta.Q
  invRC2Qover2::Float64 = 0.5/RC2Q
  QRoverRC2Q::Float64 = theta.Q * theta.R / RC2Q
  sqrtQRoverRC2Q::Float64 = sqrt(QRoverRC2Q)

  invRover2::Float64 = 0.5/theta.R
  logncG::Float64 = -0.5 * log(2 * π * RC2Q)

  @inline function lG(p::Int64, particle::LGLOPParticle, ::Nothing)
    if p == 1
      return lG1
    else
      v::Float64 = theta.C * theta.A * particle.xprev - ys[p]
      return logncG - v * invRC2Qover2 * v
    end
  end
  @inline function M!(newParticle::LGLOPParticle, rng::RNG, p::Int64,
    particle::LGLOPParticle, ::Nothing)
    if p == 1
      newParticle.x = mu1 + sqrtv0RoverRC2v0*randn(rng)
    else
      mu::Float64 = QRoverRC2Q*(theta.A*particle.x/theta.Q + theta.C*ys[p]/theta.R)
      newParticle.x = mu + sqrtQRoverRC2Q*randn(rng)
      newParticle.xprev = particle.x
    end
  end
  model::SMCModel = SMCModel(M!, lG, n, LGLOPParticle, Nothing)
  return model
end


text:
using Genie.Router
import CardgenController

route("/") do
  serve_static_file("welcome.html")
end

route("/gencard/:round::Int/:id::Int/card.ps",
      CardgenController.gencard_ps_r,
      named=:get_card_ps)
route("/gencard/:round::Int/:id::Int/card.pdf",
      CardgenController.gencard_pdf_r,
      named=:get_card_pdf)

route("/gencard/:round::Int/seq.ps",
      CardgenController.genseq_rg(:get_card_ps))
route("/gencard/:round::Int/seq.pdf",
      CardgenController.genseq_rg(:get_card_pdf))

route("/gencard/:round::Int/:id::Int/addcheck.ps",
      CardgenController.addcheck_rg(:get_card_ps))
route("/gencard/:round::Int/:id::Int/addcheck.pdf",
      CardgenController.addcheck_rg(:get_card_pdf))

max_stars_repo_name: tuckertwo/ECbingo

max_issues_repo_name: tuckertwo/ECbingo

max_forks_repo_name: tuckertwo/ECbingo