text:
using SuperFastGraphs
using Base.Test
using LightGraphs

# g = Graph(5)

# add_edge!(g, 1, 2)
# add_edge!(g, 1, 3)
# add_edge!(g, 1, 4)
# add_edge!(g, 1, 5)
# # add_edge!(g, 1, 6)
# # add_edge!(g, 1, 7)

# add_edge!(g, 2, 3)
# add_edge!(g, 2, 4)
# # add_edge!(g, 5, 8)
# # add_edge!(g, 2, 9)
# # add_edge!(g, 4, 5)
# # add_edge!(g, 8, 9)

g = Graph(6)

add_edge!(g, 1, 2)
add_edge!(g, 1, 3)
add_edge!(g, 1, 4)
add_edge!(g, 1, 5)
add_edge!(g, 2, 3)
add_edge!(g, 3, 6)
add_edge!(g, 4, 5)

g = loadgraph("./testdata/karate.lg", "graph")

@test sampleDistance(g, 1000) == 2.408199643493761

@test iFub(g, 1) == 5
@test iFub(g) == 5

cc1 = closeness_centrality(g)[1]
cc1Sample = ccSample(g, 10, 1)
@test (cc1Sample < (cc1 + 1)) & (cc1Sample > (cc1 - 1))

@test (triangleCountingDegree(g) == triangles(g))
@test (triangleForward(g) == sum(triangles(g))/3)

## TEST PREPROCESS
gDWeakly = DiGraph(6)
add_edge!(gDWeakly, 1, 2)
add_edge!(gDWeakly, 1, 3)
add_edge!(gDWeakly, 1, 4)
add_edge!(gDWeakly, 5, 1)
add_edge!(gDWeakly, 2, 3)
add_edge!(gDWeakly, 3, 6)
add_edge!(gDWeakly, 4, 5)
@test preProcess(gDWeakly) == ([[1, 4, 5], [2], [3], [6]], [6, 3, 2, 1], [10, 4, 2, 1])

gDStrongly = DiGraph(3)
add_edge!(gDStrongly, 1, 2)
add_edge!(gDStrongly, 2, 3)
add_edge!(gDStrongly, 3, 1)
@test preProcess(gDStrongly) == ([[1, 2, 3]], [3], [3])

gNotD_NotC = Graph(6)
add_edge!(gNotD_NotC, 1, 2)
add_edge!(gNotD_NotC, 1, 3)
add_edge!(gNotD_NotC, 2, 4)
add_edge!(gNotD_NotC, 2, 5)
add_edge!(gNotD_NotC, 5, 4)
@test preProcess(gNotD_NotC) == ([[1, 2, 3, 4, 5], [6]], [5, 1], [5, 1])

## TEST FASTCLOSENESSCENTRALITY
@test ((sort!(closeness_centrality(g)))[end] == topKcc(g, 1)[1][2])
@test ((sort!(closeness_centrality(g)))[end - 5] == topKcc(g, 5)[1][2])

## TEST TRIANGLE COUNTING

ext: jl

max_stars_repo_name: mameli/SuperFastGraphs

max_issues_repo_name: mameli/SuperFastGraphs

max_forks_repo_name: mameli/SuperFastGraphs


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

"""
    AssociateAdminAccount()

Sets the AWS Firewall Manager administrator account. AWS Firewall Manager must be associated with the master account of your AWS organization or associated with a member account that has the appropriate permissions. If the account ID that you submit is not an AWS Organizations master account, AWS Firewall Manager will set the appropriate permissions for the given member account. The account that you associate with AWS Firewall Manager is called the AWS Firewall Manager administrator account.

# Required Parameters
- `AdminAccount`: The AWS account ID to associate with AWS Firewall Manager as the AWS Firewall Manager administrator account. This can be an AWS Organizations master account or a member account. For more information about AWS Organizations and master accounts, see Managing the AWS Accounts in Your Organization.

"""
associate_admin_account(AdminAccount; aws_config::AWSConfig=global_aws_config()) = fms("AssociateAdminAccount", Dict{String, Any}("AdminAccount"=>AdminAccount); aws_config=aws_config)
associate_admin_account(AdminAccount, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("AssociateAdminAccount", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AdminAccount"=>AdminAccount), args)); aws_config=aws_config)

"""
    DeleteAppsList()

Permanently deletes an AWS Firewall Manager applications list.

# Required Parameters
- `ListId`: The ID of the applications list that you want to delete. You can retrieve this ID from PutAppsList, ListAppsLists, and GetAppsList.

"""
delete_apps_list(ListId; aws_config::AWSConfig=global_aws_config()) = fms("DeleteAppsList", Dict{String, Any}("ListId"=>ListId); aws_config=aws_config)
delete_apps_list(ListId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("DeleteAppsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ListId"=>ListId), args)); aws_config=aws_config)

"""
    DeleteNotificationChannel()

Deletes an AWS Firewall Manager association with the IAM role and the Amazon Simple Notification Service (SNS) topic that is used to record AWS Firewall Manager SNS logs.

"""
delete_notification_channel(; aws_config::AWSConfig=global_aws_config()) = fms("DeleteNotificationChannel"; aws_config=aws_config)
delete_notification_channel(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("DeleteNotificationChannel", args; aws_config=aws_config)

"""
    DeletePolicy()

Permanently deletes an AWS Firewall Manager policy.

# Required Parameters
- `PolicyId`: The ID of the policy that you want to delete. You can retrieve this ID from PutPolicy and ListPolicies.

# Optional Parameters
- `DeleteAllPolicyResources`: If True, the request performs cleanup according to the policy type.  For AWS WAF and Shield Advanced policies, the cleanup does the following:   Deletes rule groups created by AWS Firewall Manager   Removes web ACLs from in-scope resources   Deletes web ACLs that contain no rules or rule groups   For security group policies, the cleanup does the following for each security group in the policy:   Disassociates the security group from in-scope resources    Deletes the security group if it was created through Firewall Manager and if it's no longer associated with any resources through another policy   After the cleanup, in-scope resources are no longer protected by web ACLs in this policy. Protection of out-of-scope resources remains unchanged. Scope is determined by tags that you create and accounts that you associate with the policy. When creating the policy, if you specify that only resources in specific accounts or with specific tags are in scope of the policy, those accounts and resources are handled by the policy. All others are out of scope. If you don't specify tags or accounts, all resources are in scope.
"""
delete_policy(PolicyId; aws_config::AWSConfig=global_aws_config()) = fms("DeletePolicy", Dict{String, Any}("PolicyId"=>PolicyId); aws_config=aws_config)
delete_policy(PolicyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("DeletePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("PolicyId"=>PolicyId), args)); aws_config=aws_config)

"""
    DeleteProtocolsList()

Permanently deletes an AWS Firewall Manager protocols list.

# Required Parameters
- `ListId`: The ID of the protocols list that you want to delete. You can retrieve this ID from PutProtocolsList, ListProtocolsLists, and GetProtocolsLost.

"""
delete_protocols_list(ListId; aws_config::AWSConfig=global_aws_config()) = fms("DeleteProtocolsList", Dict{String, Any}("ListId"=>ListId); aws_config=aws_config)
delete_protocols_list(ListId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("DeleteProtocolsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ListId"=>ListId), args)); aws_config=aws_config)

"""
    DisassociateAdminAccount()

Disassociates the account that has been set as the AWS Firewall Manager administrator account. To set a different account as the administrator account, you must submit an AssociateAdminAccount request.

"""
disassociate_admin_account(; aws_config::AWSConfig=global_aws_config()) = fms("DisassociateAdminAccount"; aws_config=aws_config)
disassociate_admin_account(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("DisassociateAdminAccount", args; aws_config=aws_config)

"""
    GetAdminAccount()

Returns the AWS Organizations master account that is associated with AWS Firewall Manager as the AWS Firewall Manager administrator.

"""
get_admin_account(; aws_config::AWSConfig=global_aws_config()) = fms("GetAdminAccount"; aws_config=aws_config)
get_admin_account(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetAdminAccount", args; aws_config=aws_config)

"""
    GetAppsList()

Returns information about the specified AWS Firewall Manager applications list.

# Required Parameters
- `ListId`: The ID of the AWS Firewall Manager applications list that you want the details for.

# Optional Parameters
- `DefaultList`: Specifies whether the list to retrieve is a default list owned by AWS Firewall Manager.
"""
get_apps_list(ListId; aws_config::AWSConfig=global_aws_config()) = fms("GetAppsList", Dict{String, Any}("ListId"=>ListId); aws_config=aws_config)
get_apps_list(ListId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetAppsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ListId"=>ListId), args)); aws_config=aws_config)

"""
    GetComplianceDetail()

Returns detailed compliance information about the specified member account. Details include resources that are in and out of compliance with the specified policy. Resources are considered noncompliant for AWS WAF and Shield Advanced policies if the specified policy has not been applied to them. Resources are considered noncompliant for security group policies if they are in scope of the policy, they violate one or more of the policy rules, and remediation is disabled or not possible. Resources are considered noncompliant for Network Firewall policies if a firewall is missing in the VPC, if the firewall endpoint isn't set up in an expected Availability Zone and subnet, if a subnet created by the Firewall Manager doesn't have the expected route table, and for modifications to a firewall policy that violate the Firewall Manager policy's rules.

# Required Parameters
- `MemberAccount`: The AWS account that owns the resources that you want to get the details for.
- `PolicyId`: The ID of the policy that you want to get the details for. PolicyId is returned by PutPolicy and by ListPolicies.

"""
get_compliance_detail(MemberAccount, PolicyId; aws_config::AWSConfig=global_aws_config()) = fms("GetComplianceDetail", Dict{String, Any}("MemberAccount"=>MemberAccount, "PolicyId"=>PolicyId); aws_config=aws_config)
get_compliance_detail(MemberAccount, PolicyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetComplianceDetail", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MemberAccount"=>MemberAccount, "PolicyId"=>PolicyId), args)); aws_config=aws_config)

"""
    GetNotificationChannel()

Information about the Amazon Simple Notification Service (SNS) topic that is used to record AWS Firewall Manager SNS logs.

"""
get_notification_channel(; aws_config::AWSConfig=global_aws_config()) = fms("GetNotificationChannel"; aws_config=aws_config)
get_notification_channel(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetNotificationChannel", args; aws_config=aws_config)

"""
    GetPolicy()

Returns information about the specified AWS Firewall Manager policy.

# Required Parameters
- `PolicyId`: The ID of the AWS Firewall Manager policy that you want the details for.

"""
get_policy(PolicyId; aws_config::AWSConfig=global_aws_config()) = fms("GetPolicy", Dict{String, Any}("PolicyId"=>PolicyId); aws_config=aws_config)
get_policy(PolicyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("PolicyId"=>PolicyId), args)); aws_config=aws_config)

"""
    GetProtectionStatus()

If you created a Shield Advanced policy, returns policy-level attack summary information in the event of a potential DDoS attack. Other policy types are currently unsupported.

# Required Parameters
- `PolicyId`: The ID of the policy for which you want to get the attack information.

# Optional Parameters
- `EndTime`: The end of the time period to query for the attacks. This is a timestamp type. The request syntax listing indicates a number type because the default used by AWS Firewall Manager is Unix time in seconds. However, any valid timestamp format is allowed.
- `MaxResults`: Specifies the number of objects that you want AWS Firewall Manager to return for this request. If you have more objects than the number that you specify for MaxResults, the response includes a NextToken value that you can use to get another batch of objects.
- `MemberAccountId`: The AWS account that is in scope of the policy that you want to get the details for.
- `NextToken`: If you specify a value for MaxResults and you have more objects than the number that you specify for MaxResults, AWS Firewall Manager returns a NextToken value in the response, which you can use to retrieve another group of objects. For the second and subsequent GetProtectionStatus requests, specify the value of NextToken from the previous response to get information about another batch of objects.
- `StartTime`: The start of the time period to query for the attacks. This is a timestamp type. The request syntax listing indicates a number type because the default used by AWS Firewall Manager is Unix time in seconds. However, any valid timestamp format is allowed.
"""
get_protection_status(PolicyId; aws_config::AWSConfig=global_aws_config()) = fms("GetProtectionStatus", Dict{String, Any}("PolicyId"=>PolicyId); aws_config=aws_config)
get_protection_status(PolicyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetProtectionStatus", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("PolicyId"=>PolicyId), args)); aws_config=aws_config)

"""
    GetProtocolsList()

Returns information about the specified AWS Firewall Manager protocols list.

# Required Parameters
- `ListId`: The ID of the AWS Firewall Manager protocols list that you want the details for.

# Optional Parameters
- `DefaultList`: Specifies whether the list to retrieve is a default list owned by AWS Firewall Manager.
"""
get_protocols_list(ListId; aws_config::AWSConfig=global_aws_config()) = fms("GetProtocolsList", Dict{String, Any}("ListId"=>ListId); aws_config=aws_config)
get_protocols_list(ListId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetProtocolsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ListId"=>ListId), args)); aws_config=aws_config)

"""
    GetViolationDetails()

Retrieves violations for a resource based on the specified AWS Firewall Manager policy and AWS account.

# Required Parameters
- `MemberAccount`: The AWS account ID that you want the details for.
- `PolicyId`: The ID of the AWS Firewall Manager policy that you want the details for. This currently only supports security group content audit policies.
- `ResourceId`: The ID of the resource that has violations.
- `ResourceType`: The resource type. This is in the format shown in the AWS Resource Types Reference. Supported resource types are: AWS::EC2::Instance, AWS::EC2::NetworkInterface, AWS::EC2::SecurityGroup, AWS::NetworkFirewall::FirewallPolicy, and AWS::EC2::Subnet.

"""
get_violation_details(MemberAccount, PolicyId, ResourceId, ResourceType; aws_config::AWSConfig=global_aws_config()) = fms("GetViolationDetails", Dict{String, Any}("MemberAccount"=>MemberAccount, "PolicyId"=>PolicyId, "ResourceId"=>ResourceId, "ResourceType"=>ResourceType); aws_config=aws_config)
get_violation_details(MemberAccount, PolicyId, ResourceId, ResourceType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("GetViolationDetails", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MemberAccount"=>MemberAccount, "PolicyId"=>PolicyId, "ResourceId"=>ResourceId, "ResourceType"=>ResourceType), args)); aws_config=aws_config)

"""
    ListAppsLists()

Returns an array of AppsListDataSummary objects.

# Required Parameters
- `MaxResults`: The maximum number of objects that you want AWS Firewall Manager to return for this request. If more objects are available, in the response, AWS Firewall Manager provides a NextToken value that you can use in a subsequent call to get the next batch of objects. If you don't specify this, AWS Firewall Manager returns all available objects.

# Optional Parameters
- `DefaultLists`: Specifies whether the lists to retrieve are default lists owned by AWS Firewall Manager.
- `NextToken`: If you specify a value for MaxResults in your list request, and you have more objects than the maximum, AWS Firewall Manager returns this token in the response. For all but the first request, you provide the token returned by the prior request in the request parameters, to retrieve the next batch of objects.
"""
list_apps_lists(MaxResults; aws_config::AWSConfig=global_aws_config()) = fms("ListAppsLists", Dict{String, Any}("MaxResults"=>MaxResults); aws_config=aws_config)
list_apps_lists(MaxResults, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListAppsLists", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MaxResults"=>MaxResults), args)); aws_config=aws_config)

"""
    ListComplianceStatus()

Returns an array of PolicyComplianceStatus objects. Use PolicyComplianceStatus to get a summary of which member accounts are protected by the specified policy.

# Required Parameters
- `PolicyId`: The ID of the AWS Firewall Manager policy that you want the details for.

# Optional Parameters
- `MaxResults`: Specifies the number of PolicyComplianceStatus objects that you want AWS Firewall Manager to return for this request. If you have more PolicyComplianceStatus objects than the number that you specify for MaxResults, the response includes a NextToken value that you can use to get another batch of PolicyComplianceStatus objects.
- `NextToken`: If you specify a value for MaxResults and you have more PolicyComplianceStatus objects than the number that you specify for MaxResults, AWS Firewall Manager returns a NextToken value in the response that allows you to list another group of PolicyComplianceStatus objects. For the second and subsequent ListComplianceStatus requests, specify the value of NextToken from the previous response to get information about another batch of PolicyComplianceStatus objects.
"""
list_compliance_status(PolicyId; aws_config::AWSConfig=global_aws_config()) = fms("ListComplianceStatus", Dict{String, Any}("PolicyId"=>PolicyId); aws_config=aws_config)
list_compliance_status(PolicyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListComplianceStatus", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("PolicyId"=>PolicyId), args)); aws_config=aws_config)

"""
    ListMemberAccounts()

Returns a MemberAccounts object that lists the member accounts in the administrator's AWS organization. The ListMemberAccounts must be submitted by the account that is set as the AWS Firewall Manager administrator.

# Optional Parameters
- `MaxResults`: Specifies the number of member account IDs that you want AWS Firewall Manager to return for this request. If you have more IDs than the number that you specify for MaxResults, the response includes a NextToken value that you can use to get another batch of member account IDs.
- `NextToken`: If you specify a value for MaxResults and you have more account IDs than the number that you specify for MaxResults, AWS Firewall Manager returns a NextToken value in the response that allows you to list another group of IDs. For the second and subsequent ListMemberAccountsRequest requests, specify the value of NextToken from the previous response to get information about another batch of member account IDs.
"""
list_member_accounts(; aws_config::AWSConfig=global_aws_config()) = fms("ListMemberAccounts"; aws_config=aws_config)
list_member_accounts(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListMemberAccounts", args; aws_config=aws_config)

"""
    ListPolicies()

Returns an array of PolicySummary objects.

# Optional Parameters
- `MaxResults`: Specifies the number of PolicySummary objects that you want AWS Firewall Manager to return for this request. If you have more PolicySummary objects than the number that you specify for MaxResults, the response includes a NextToken value that you can use to get another batch of PolicySummary objects.
- `NextToken`: If you specify a value for MaxResults and you have more PolicySummary objects than the number that you specify for MaxResults, AWS Firewall Manager returns a NextToken value in the response that allows you to list another group of PolicySummary objects. For the second and subsequent ListPolicies requests, specify the value of NextToken from the previous response to get information about another batch of PolicySummary objects.
"""
list_policies(; aws_config::AWSConfig=global_aws_config()) = fms("ListPolicies"; aws_config=aws_config)
list_policies(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListPolicies", args; aws_config=aws_config)

"""
    ListProtocolsLists()

Returns an array of ProtocolsListDataSummary objects.

# Required Parameters
- `MaxResults`: The maximum number of objects that you want AWS Firewall Manager to return for this request. If more objects are available, in the response, AWS Firewall Manager provides a NextToken value that you can use in a subsequent call to get the next batch of objects. If you don't specify this, AWS Firewall Manager returns all available objects.

# Optional Parameters
- `DefaultLists`: Specifies whether the lists to retrieve are default lists owned by AWS Firewall Manager.
- `NextToken`: If you specify a value for MaxResults in your list request, and you have more objects than the maximum, AWS Firewall Manager returns this token in the response. For all but the first request, you provide the token returned by the prior request in the request parameters, to retrieve the next batch of objects.
"""
list_protocols_lists(MaxResults; aws_config::AWSConfig=global_aws_config()) = fms("ListProtocolsLists", Dict{String, Any}("MaxResults"=>MaxResults); aws_config=aws_config)
list_protocols_lists(MaxResults, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListProtocolsLists", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MaxResults"=>MaxResults), args)); aws_config=aws_config)

"""
    ListTagsForResource()

Retrieves the list of tags for the specified AWS resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the resource to return tags for. The AWS Firewall Manager resources that support tagging are policies, applications lists, and protocols lists.

"""
list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = fms("ListTagsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    PutAppsList()

Creates an AWS Firewall Manager applications list.

# Required Parameters
- `AppsList`: The details of the AWS Firewall Manager applications list to be created.

# Optional Parameters
- `TagList`: The tags associated with the resource.
"""
put_apps_list(AppsList; aws_config::AWSConfig=global_aws_config()) = fms("PutAppsList", Dict{String, Any}("AppsList"=>AppsList); aws_config=aws_config)
put_apps_list(AppsList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("PutAppsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AppsList"=>AppsList), args)); aws_config=aws_config)

"""
    PutNotificationChannel()

Designates the IAM role and Amazon Simple Notification Service (SNS) topic that AWS Firewall Manager uses to record SNS logs. To perform this action outside of the console, you must configure the SNS topic to allow the Firewall Manager role AWSServiceRoleForFMS to publish SNS logs. For more information, see Firewall Manager required permissions for API actions in the AWS Firewall Manager Developer Guide.

# Required Parameters
- `SnsRoleName`: The Amazon Resource Name (ARN) of the IAM role that allows Amazon SNS to record AWS Firewall Manager activity.
- `SnsTopicArn`: The Amazon Resource Name (ARN) of the SNS topic that collects notifications from AWS Firewall Manager.

"""
put_notification_channel(SnsRoleName, SnsTopicArn; aws_config::AWSConfig=global_aws_config()) = fms("PutNotificationChannel", Dict{String, Any}("SnsRoleName"=>SnsRoleName, "SnsTopicArn"=>SnsTopicArn); aws_config=aws_config)
put_notification_channel(SnsRoleName, SnsTopicArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("PutNotificationChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SnsRoleName"=>SnsRoleName, "SnsTopicArn"=>SnsTopicArn), args)); aws_config=aws_config)

"""
    PutPolicy()

Creates an AWS Firewall Manager policy. Firewall Manager provides the following types of policies:    An AWS WAF policy (type WAFV2), which defines rule groups to run first in the corresponding AWS WAF web ACL and rule groups to run last in the web ACL.   An AWS WAF Classic policy (type WAF), which defines a rule group.    A Shield Advanced policy, which applies Shield Advanced protection to specified accounts and resources.   A security group policy, which manages VPC security groups across your AWS organization.    An AWS Network Firewall policy, which provides firewall rules to filter network traffic in specified Amazon VPCs.   Each policy is specific to one of the types. If you want to enforce more than one policy type across accounts, create multiple policies. You can create multiple policies for each type. You must be subscribed to Shield Advanced to create a Shield Advanced policy. For more information about subscribing to Shield Advanced, see CreateSubscription.

# Required Parameters
- `Policy`: The details of the AWS Firewall Manager policy to be created.

# Optional Parameters
- `TagList`: The tags to add to the AWS resource.
"""
put_policy(Policy; aws_config::AWSConfig=global_aws_config()) = fms("PutPolicy", Dict{String, Any}("Policy"=>Policy); aws_config=aws_config)
put_policy(Policy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("PutPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Policy"=>Policy), args)); aws_config=aws_config)

"""
    PutProtocolsList()

Creates an AWS Firewall Manager protocols list.

# Required Parameters
- `ProtocolsList`: The details of the AWS Firewall Manager protocols list to be created.

# Optional Parameters
- `TagList`: The tags associated with the resource.
"""
put_protocols_list(ProtocolsList; aws_config::AWSConfig=global_aws_config()) = fms("PutProtocolsList", Dict{String, Any}("ProtocolsList"=>ProtocolsList); aws_config=aws_config)
put_protocols_list(ProtocolsList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("PutProtocolsList", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ProtocolsList"=>ProtocolsList), args)); aws_config=aws_config)

"""
    TagResource()

Adds one or more tags to an AWS resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the resource to return tags for. The AWS Firewall Manager resources that support tagging are policies, applications lists, and protocols lists.
- `TagList`: The tags to add to the resource.

"""
tag_resource(ResourceArn, TagList; aws_config::AWSConfig=global_aws_config()) = fms("TagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagList"=>TagList); aws_config=aws_config)
tag_resource(ResourceArn, TagList, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagList"=>TagList), args)); aws_config=aws_config)

"""
    UntagResource()

Removes one or more tags from an AWS resource.

# Required Parameters
- `ResourceArn`: The Amazon Resource Name (ARN) of the resource to return tags for. The AWS Firewall Manager resources that support tagging are policies, applications lists, and protocols lists.
- `TagKeys`: The keys of the tags to remove from the resource.

"""
untag_resource(ResourceArn, TagKeys; aws_config::AWSConfig=global_aws_config()) = fms("UntagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(ResourceArn, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = fms("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys), args)); aws_config=aws_config)


text:
# ---
# title: 1443. Minimum Time to Collect All Apples in a Tree
# id: problem1443
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/minimum-time-to-collect-all-apples-in-a-tree/description/>
# hidden: true
# ---
#
# Given an undirected tree consisting of `n` vertices numbered from `0` to
# `n-1`, which has some apples in their vertices. You spend 1 second to walk
# over one edge of the tree. _Return the minimum time in seconds you have to
# spend to collect all apples in the tree, starting at **vertex 0** and coming
# back to this vertex._
#
# The edges of the undirected tree are given in the array `edges`, where
# `edges[i] = [ai, bi]` means that exists an edge connecting the vertices `ai`
# and `bi`. Additionally, there is a boolean array `hasApple`, where
# `hasApple[i] = true` means that vertex `i` has an apple; otherwise, it does
# not have any apple.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/23/min_time_collect_apple_1.png)**
#
#
#
#     Input: n = 7, edges = [[0,1],[0,2],[1,4],[1,5],[2,3],[2,6]], hasApple = [false,false,true,false,true,true,false]
#     Output: 8
#     Explanation: The figure above represents the given tree where red vertices have an apple. One optimal path to collect all apples is shown by the green arrows.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/23/min_time_collect_apple_2.png)**
#
#
#
#     Input: n = 7, edges = [[0,1],[0,2],[1,4],[1,5],[2,3],[2,6]], hasApple = [false,false,true,false,false,true,false]
#     Output: 6
#     Explanation: The figure above represents the given tree where red vertices have an apple. One optimal path to collect all apples is shown by the green arrows.
#
#
# **Example 3:**
#
#
#
#     Input: n = 7, edges = [[0,1],[0,2],[1,4],[1,5],[2,3],[2,6]], hasApple = [false,false,false,false,false,false,false]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `edges.length == n - 1`
#   * `edges[i].length == 2`
#   * `0 <= ai < bi <= n - 1`
#   * `fromi < toi`
#   * `hasApple.length == n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# 1D diffusion problem

# Packages and inclusions
using ModelingToolkit, MethodOfLines, LinearAlgebra, Test
using ModelingToolkit: Differential
using NonlinearSolve
using DomainSets

# Handrolled discrete laplace problem
@testset "Hand rolled 1D Laplace" begin
    @variables a b c d
    eqs = [0 ~ a - 1,
        0 ~ a + c - 2 * b,
        0 ~ b + d - 2 * c,
        0 ~ d - 1]
    @named ns = NonlinearSystem(eqs, [a, b, c, d], [])
    f = eval(generate_function(ns, [a, b, c, d])[2])
    prob = NonlinearProblem(ns, zeros(4), [])
    sol = NonlinearSolve.solve(prob, NewtonRaphson())

    @test sol.u ≈ ones(4)
end

# Laplace's Equation, same as above but with MOL discretization
@testset "1D Laplace - constant solution" begin
    @parameters x
    @variables u(..)
    Dxx = Differential(x)^2

    eq = Dxx(u(x)) ~ 0
    dx = 1 / 3

    bcs = [u(0) ~ 1,
        u(1) ~ 1]

    # Space and time domains
    domains = [x ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem([eq], bcs, domains, [x], [u(x)])
    discretization = MOLFiniteDifference([x => dx], nothing, approx_order=2)
    prob = discretize(pdesys, discretization)
    grid = get_discrete(pdesys, discretization)
    sol = NonlinearSolve.solve(prob, NewtonRaphson())
    sol′ = map(d -> sol[d], grid[u(x)])
    @test sol′ ≈ ones(4)
end

# Laplace's Equation, linear solution
@testset "1D Laplace - linear solution" begin
    @parameters x
    @variables u(..)
    Dxx = Differential(x)^2

    eq = Dxx(u(x)) ~ 0
    dx = 0.1

    bcs = [u(0) ~ 1,
        u(1) ~ 2]

    # Space and time domains
    domains = [x ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem([eq], bcs, domains, [x], [u(x)])
    discretization = MOLFiniteDifference([x => dx], nothing, approx_order=2)
    grid = get_discrete(pdesys, discretization)

    prob = discretize(pdesys, discretization)
    sol = NonlinearSolve.solve(prob, NewtonRaphson())
    sol′ = map(d -> sol[d], grid[u(x)])

    @test sol′ ≈ 1.0:0.1:2.0
end

# 2D heat
@testset "2D heat equation" begin
    @parameters x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2

    eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
    dx = 0.1
    dy = 0.1

    bcs = [u(0, y) ~ x * y,
        u(1, y) ~ x * y,
        u(x, 0) ~ x * y,
        u(x, 1) ~ x * y]

    # Space and time domains
    domains = [x ∈ Interval(0.0, 1.0),
        y ∈ Interval(0.0, 1.0)]

    @named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])

    # Note that we pass in `nothing` for the time variable `t` here since we
    # are creating a stationary problem without a dependence on time, only space.
    discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)

    prob = discretize(pdesys, discretization)
    grid = get_discrete(pdesys, discretization)
    sol = NonlinearSolve.solve(prob, NewtonRaphson())
    xs, ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
    u_sol = map(d -> sol[d], grid[u(x, y)])

    # test boundary
    @test all(abs.(u_sol[:, 1]) .< eps(Float32))
    @test all(abs.(u_sol[1, :]) .< eps(Float32))
    @test u_sol[2:end-1, end] ≈ (0.0:dy:1.0)[2:end-1]
    @test u_sol[end, 2:end-1] ≈ (0.0:dx:1.0)[2:end-1]

    # test interior with finite differences
    interior = CartesianIndices((axes(xs)[1], axes(ys)[1]))[2:end-1, 2:end-1]
    fd = map(interior) do I
        abs(u_sol[(I-CartesianIndex(1, 0))] + u_sol[(I+CartesianIndex(1, 0))] - 2 * u_sol[I]) < 0.01
    end#
    @test all(fd)
end


text:
# using Plots
# pyplot()

using RecipesBase

"""
`Basics.plot()`  ... plots various quantities, often in a new window.

  + `("radial potentials", potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N::Int64 = 0)`
    ... to plot one or more radial potentials, and where N::Int64 describes the number of grid points to be considered.
        call:  using Plots; pyplot()    ... to access this method by plot(...)
"""
function Basics.plot(sa::String, potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N::Int64 = 0)
    error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...)")
end

@recipe function f(sa::String, potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N = 0)
    !(sa == "radial potentials")   &&   error("Unsupported keystring = $sa")
    wa = Float64[];   wc = [NaN for i=1:N];   labels = String[];   np = length(potentials)
    for  pot in potentials
        wb = wc
        nx = min(length(pot.Zr), N);   wb[1:nx] = pot.Zr[1:nx]
        append!(wa, wb)
        push!(labels, pot.name)
    end
    x = grid.r[1:N];     y = reshape(wa, (N, np))
    label --> reshape(labels, (1,np))
    x, y
end

"""
  + `("radial orbitals: large", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
    ... to plot the large component of one or more radial orbitals, and where N::Int64 describes the number of grid
        points to be considered.
  + `("radial orbitals: small", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
    ... to plot the small component of one or more radial orbitals.
  + `("radial orbitals: both",  orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
    ... to plot the large and small component of one or more radial orbitals.

    call:  using Plots; pyplot()    ... to access this method by plot(...)
"""
function Basics.plot(sa::String, orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N = 0)
    error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...)")
end

@recipe function f(sa::String, orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N = 0)
    wa = Float64[];   wc = [NaN for i=1:N];   labels = String[];   np = length(orbitals)

    if       sa == "radial orbitals: large"
        for  orb in orbitals
            wb = wc
            nx = min(length(orb.P), N);   wb[1:nx] = orb.P[1:nx]
            append!(wa, wb)
            push!(labels, "$(orb.subshell):large")
        end
        ##x print("labels = $labels ")
        x = grid.r[1:N];     y = reshape(wa, (N, np))

    elseif   sa == "radial orbitals: small"
        for  orb in orbitals
            wb = wc
            nx = min(length(orb.Q), N);   wb[1:nx] = orb.Q[1:nx]
            append!(wa, wb)
            push!(labels, "$(orb.subshell):small")
        end
        ##x print("labels = $labels ")
        x = grid.r[1:N];     y = reshape(wa, (N, np))

    elseif   sa == "radial orbitals: both"
        for  orb in orbitals
            wb = wc;    nx = min(length(orb.P), N);   wb[1:nx] = orb.P[1:nx];    append!(wa, wb)
            wb = wc;    nx = min(length(orb.Q), N);   wb[1:nx] = orb.Q[1:nx];    append!(wa, wb)
            push!(labels, "$(orb.subshell):large");   push!(labels, "$(orb.subshell):small")
        end
        ##x print("labels = $labels ")
        x = grid.r[1:N];     y = reshape(wa, (N, 2np))

    else   error("Unsupported keystring = $sa")
    end
    label --> permutedims(labels)
    x, y
end

"""
  + `("spectrum: transition rates over energy", lines::Array{PhotoEmission.Line,1})`
    ... to plot the transition rates of all lines as function of their transition energies. The plot is shown in a new
        window but nothing is returned otherwise. **Not yet implemented !**

  + `("spectrum: oscillator strength over energy, emission", lines::Array{PhotoEmission.Line,1})` or
    `("spectrum: oscillator strength over energy, absorption", lines::Array{PhotoEmission.Line,1})`
    ... to plot the absorption oscillator strength of all lines as function of their transition energies. Again, a new
        window is opened but nothing returned by this method. **Not yet implemented !**
"""
function Basics.plot(sa::String, lines::Array{PhotoEmission.Line,1})
    error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...) ... not yet implemented !")
end

"""
  + `("spectrum: transition rates over energy, Gaussian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)` or
    `("spectrum: transition rates over energy, Lorentzian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)`
    ... to plot the transition rates of all lines as function of their transition energies but with a Gaussian or Lorentzian
        distribution. Again, a new window is opened but nothing returned by this method. It still need to be decided how
        the widths (and, perhaps, other parameters) will be communicated to the method. **Not yet implemented !**
"""
function Basics.plot(sa::String, lines::Array{PhotoEmission.Line,1}, widths::Float64)
    error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...) ... not yet implemented !")
end


text:
I1 = 1 - exp(-1)
f1(x) = (1-x)*exp(-x) # Gradshteyn 4.351(1)

for endpt in [neither, both, left, right]
    x, w = logweight(10, 0, endpt)
    @test I1 ≈ sum(w .* f1.(x))
    x, w = logweight(10, 0.0, endpt)
    @test I1 ≈ sum(w .* f1.(x))
end

I2 = (π/3)^2 * 2 * sqrt(3)
f2(x) = (1-x^2) / (1+x^3) # Gradshteyn 4.255(1), p=3/2, q=3
for endpt in [neither, both, left, right]
    x, w = logweight(10, -1/2, endpt)
    @test I2 ≈ sum(w .* f2.(x))
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --compile=min -O0 --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using GuillotineModels.Data: read_from_file, write_to_file, is_collection

function number_path(s, i)
	if match(r"\.", s) !== nothing
		replace(s, '.' => "_$i."; count = 1)
	else
		s * "_$i"
	end
end

function main(ARGS = ARGS)
	if length(ARGS) != 4
		println("USAGE: ./instance_converter.jl <input format> <input file> <output format> <output file>")
		println("If the input format has multiple instances, and the output format just support a single instance, then an output file will be created for each instance and '_X' (where X is the instance number) will be appended before the first dot in the output filename (or just before the end, if necessary).")
		exit()
	end
	in_format, in_file, out_format, out_file = ARGS
	in_format_val = Val(Symbol(in_format))
	instance = read_from_file(in_format_val, in_file)
	out_format_val = Val(Symbol(out_format))
	if is_collection(in_format_val) && !is_collection(out_format_val)
		for (index, inst) in enumerate(instance)
			write_to_file(out_format_val, inst, number_path(out_file, index))
		end
	else
		write_to_file(out_format_val, instance, out_file)
	end
end

function furini2016_convert(input_path, output_path)
	THOMOPULOS_THESIS_INSTANCES = vcat(
		# unweighted
		"gcut" .* string.(1:12)
		, String.(split("wang20 2s 3s A1s A2s STS2s STS4s"))
		, ["OF1", "OF2", "W", "CHL1s", "CHL2s"]
		, "A" .* string.(3:5)
		, "CHL" .* string.(5:7)
		, ["CU1", "CU2"]
		, "Hchl" .* split("3s 4s 6s 7s 8s")
		# weighted
		, "cgcut" .* string.(1:3)
		, "okp" .* string.(1:5)
		, String.(split("HH 2 3 A1 A2 STS2 STS4 CHL1 CHL2"))
		, "CW" .* string.(1:3)
		, "Hchl" .* ["2", "9"]
	) :: Vector{String}
	for inst_name in THOMOPULOS_THESIS_INSTANCES
		main([
			"Classic_G2KP",
			joinpath(input_path, inst_name),
			"Simple_CPG_SLOPP",
			joinpath(output_path, inst_name)
		])
	end
end

main()

# Comment the `main()` call and uncomment the below to have the script
# making all 59 Furini2016 instances (from Classic_G2KP format) into the
# Simple_CPG_SLOPP format.

#furini2016_convert(ARGS...)


text:
#include("../src/LatticeDiracOperators.jl")

using Gaugefields
using LinearAlgebra
using InteractiveUtils
using Random
using LatticeDiracOperators
import Gaugefields:Initialize_4DGaugefields
import LatticeDiracOperators.SSmodule:shiftedbicg_inSS,shiftedbicgstab_inSS,shiftedbicgstab,
            shiftedbicg_Frommer2003,shiftedbicg_Frommer2003_seed,shiftedbicg_Frommer2003_G_seed

function MDtest!(gauge_action,U,Dim,fermi_action,η,ξ)
    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients.
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 10
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0
    Random.seed!(123)

    numtrj = 6
    for itrj = 1:numtrj
       #@code_warntype MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
        #error("cc")
        @time accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("$itrj plaq_t = $plaq_t")
        println("acceptance ratio ",numaccepted/itrj)
    end
end

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC #evaluate_GaugeAction(gauge_action,U) = tr(evaluate_GaugeAction_untraced(gauge_action,U))
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
    Δτ = 1/MDsteps
    NC,_,NN... = size(U[1])

    for μ=1:Dim
        pwork = gauss_distribution(prod(NN)*(NC^2-1))
        substitute_U!(p[μ],pwork)
    end

    #println(p[1][1,1,1,1,1,1])

    #gauss_distribution!(p)

    substitute_U!(Uold,U)
    gauss_sampling_in_action!(ξ,U,fermi_action)
    sample_pseudofermions!(η,U,fermi_action,ξ)
    Sfold = real(dot(ξ,ξ))
    println("Sfold = $Sfold")

    #@code_warntype calc_action(gauge_action,U,p)

    Sold = calc_action(gauge_action,U,p) + Sfold
    println("Sold = ",Sold)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action)
        #println(" U1 = ", U[1][1,1,1,1,1,1])
        #        println(" p = ", p[1][1,1,1,1,1])
        P_update_fermion!(U,p,1.0,Δτ,Dim,gauge_action,fermi_action,η)
        #error("dd")

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end
    Sfnew = evaluate_FermiAction(fermi_action,U,η)
    println("Sfnew = $Sfnew")
    Snew = calc_action(gauge_action,U,p) + Sfnew

    println("Sold = $Sold, Snew = $Snew")
    println("Snew - Sold = $(Snew-Sold)")

    accept = exp(Sold - Snew) >= rand()

    #ratio = min(1,exp(Snew-Sold))
    if accept != true #rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,gauge_action,μ,U)
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function P_update_fermion!(U,p,ϵ,Δτ,Dim,gauge_action,fermi_action,η)  # p -> p +factor*U*dSdUμ
    #NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    UdSfdUμ = temps[1:Dim]
    factor =  -ϵ*Δτ

    #calc_p_UdSfdU!(p,fermi_action,U,η,factor)

    calc_UdSfdU!(UdSfdUμ,fermi_action,U,η)

    for μ=1:Dim
        Traceless_antihermitian_add!(p[μ],factor,UdSfdUμ[μ])
        #println(" p[μ] = ", p[μ][1,1,1,1,1])
    end

end

function test1()
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0
    Dim = 4
    NC = 3

    U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    #U  =Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 5.5/2
    push!(gauge_action,β,plaqloop)

    show(gauge_action)

    #Initialize_4DWilsonFermion(U[1])
    #error("init")
    #@code_warntype Initialize_pseudofermion_fields(U[1],"Wilson")
    #error("init")

    #x = Initialize_pseudofermion_fields(U[1],params)

    #x = Initialize_4DWilsonFermion(U[1])

    #=
    xplus = shift_fermion(x,1)

    ix = 1
    iy = 1
    iz = 1
    it = 1
    ialpha  =1
    xplus[1,ix,iy,iz,it,ialpha]
    println("dd")

    @inbounds for ix= 1:2
         @code_llvm x.f[1,ix,iy,iz,it,ialpha]
    end

    vv = rand(ComplexF64,2,2,2,2,2,2)
    @inbounds for ix= 1:2
        @code_llvm vv[1,ix,iy,iz,it,ialpha]
    end

    #error("err")
    @code_warntype shift_fermion(x,1)
    #error("2")
    ix = 1
    iy = 1
    iz = 1
    it = 1
    ialpha  =1
    @time @inbounds x.f[1,ix,iy,iz,it,ialpha]
    @code_llvm   @inbounds x.f[1,ix,iy,iz,it,ialpha]
    println("dd")
    @code_lowered   Base.getindex(x,1,ix,iy,iz,it,ialpha)
    println("dt")
    println(x[1,ix,iy,iz,it,ialpha])
    @code_llvm   @inbounds  x[1,ix,iy,iz,it,ialpha]

    @code_llvm  x[1,ix,iy,iz,it,ialpha]

    @code_llvm xplus[1,ix,iy,iz,it,ialpha]
    println(xplus[1,ix,iy,iz,it,ialpha])
    @code_llvm xplus[1,ix,iy,iz,it,ialpha]

    @code_warntype xplus[1,ix,iy,iz,it,ialpha]

    @time xplus[1,ix,iy,iz,it,ialpha]
    for ix = 1:4
        @time xplus.parent[1,ix,iy,iz,it,ialpha]
       # @time xplus[1,ix,iy,iz,it,ialpha]
    end
    #error("dd")
    #return

    =#

    x = Initialize_pseudofermion_fields(U[1],"Wilson",nowing = true)

    params = Dict()
    params["Dirac_operator"] = "Wilson"
    M = -1
    params["κ"] = 1/(2*Dim+M) #0.141139/2
    params["eps_CG"] = 1.0e-19
    params["verbose_level"] = 2
    #params["method_CG"] = "preconditiond_bicgstab"
    #params["method_CG"] = "bicgstab"
    params["faster version"] = true
    params["method_CG"] = "bicg"
    D = Dirac_operator(U,x,params)

    parameters_action = Dict()
    fermi_action = FermiAction(D,parameters_action)
    #gauss_sampling_in_action!(x,U,fermi_action)
    #println("Sfold = ", dot(x,x))
    y = similar(x)

    MDtest!(gauge_action,U,Dim,fermi_action,x,y)

    #g5D = γ5D(D(U))
    g5D = D(U)

    radius = 0.2
    origin = 0
    #es,residuals,eigenvectors,is = eigensystem_old(g5D,x,length(x),radius,origin)
    es,residuals,eigenvectors,is = eigensystem(g5D,x,radius,origin,Z2_distribution_fermion!,
                uniform_distribution_fermion!,α=1,epsCG=1e-20,κ=1.5)
    #println(result[1])
    fp = open("eigenvalues_GGwilson_SS.txt","w")
    println("es")
    println(es)
    println("es")
    println(residuals)
    for i = 1:length(es)
        println(fp,real(es[i]),"\t",imag(es[i]),"\t",real(residuals[i]))
        println(real(es[i]),"\t",imag(es[i]),"\t",real(residuals[i]))
    end
    close(fp)

    Dmat = construct_sparsematrix(g5D)
    e2,v = eigen(Matrix(Dmat))
    fp = open("eigenvalues_GGwilson_2.txt","w")
    for i = 1:length(e2)
        println(fp,real(e2[i]),"\t",imag(e2[i]))
        println(real(e2[i]),"\t",imag(e2[i]))
    end
    close(fp)

    return

    b = similar(x)
    for i in 1:length(b.f)
        b.f[i] = rand()
    end

    #=
    y = zero(x)

    for i= 1:length(e2)

        c = 0
        for j=1:length(e2)
            c += v[j,i]'*b.f[j]
        end

        for j= 1:length(e2)
            y.f[j] += sign(e2[i])*v[j,i]*c
        end
    end
    =#

    Nq = 150
    α = 0.1
    radius = 0.995
    origin_p = 1.01
    origin_m = -1.01

    vec_z = zeros(ComplexF64,2Nq)
    vec_w = zeros(ComplexF64,2Nq)
    for j in 1:Nq
        θj = (2π/Nq)*(j-1/2) # π/Nq   2π-π/Nq
        vec_z[j] = origin_p + radius*(cos(θj)+im*α*sin(θj))
        vec_w[j] = α*cos(θj)+im*sin(θj)

        vec_z[j+Nq] = origin_m + radius*(cos(θj)+im*α*sin(θj))
        vec_w[j+Nq] = α*cos(θj)+im*sin(θj)
    end

    y2 = zero(x)
    zseedin = origin_p+0.05im
    eps = 1e-20
    @time vec_x = shiftedbicg_Frommer2003_G_seed(vec_z,g5D,b,eps =eps,zseedin = zseedin) #(z I - A) x = b

    for j in 1:Nq
        coeff = radius*vec_w[j]/Nq
        axpy!(coeff,vec_x[j],y2)
        coeff = -radius*vec_w[j+Nq]/Nq
        axpy!(coeff,vec_x[j+Nq],y2)
        #vec_Sk[:,ik] += ρ*vec_w[j]*vec_z[j]^k*vec_y[:,j]/Nq
    end

    println(dot(b,b),"\t",dot(y2,y2))
    value = abs(dot(b,b) - dot(y2,y2))/abs(dot(b,b))
    println("value = $value")
    return

    fp = open("signHx_$(Nq).txt","w")

    for i= 1:length(e2)
        println(fp,real(y.f[i]),"\t",imag(y.f[i]),"\t",
                    real(y2.f[i]),"\t",imag(y2.f[i])
                    )
    end
    close(fp)

    println(dot(b,b),"\t",dot(y2,y2))
    value = abs(dot(b,b) - dot(y2,y2))/abs(dot(b,b))
    println("value = $value")

    axpby!(1,y,-1,y2)
    println(sqrt(real(dot(y2,y2)))/sqrt(real(dot(y,y))))

    return

    σi = zeros(ComplexF64,Nq)
    rmin = -2.0
    rmax = -rmin
    delta = 0.05
    dr = (rmax-rmin)/(div(Nq,2)-1)
    for j = 1:div(Nq,2)
        σi[j] =  (j-1)*dr + rmin + im*delta
        σi[j+div(Nq,2)] = conj(σi[j])
    end

    y2 = zero(x)
    zseedin = 1+0.05im

    vec_x = shiftedbicg_Frommer2003_G_seed(σi,g5D,b,eps =1e-20,zseedin = zseedin) #(z I - A) x = b

    for j = 1:div(Nq,2)
        si = ifelse( real(σi[j]) > 0,1,-1)
        axpy!(-si*dr/(2π*im),vec_x[j],y2)
        axpy!(si*dr/(2π*im),vec_x[j+div(Nq,2)],y2)
    end

    fp = open("signHx.txt","w")

    for i= 1:length(e2)
        println(fp,real(y.f[i]),"\t",imag(y.f[i]),"\t",
                    real(y2.f[i]),"\t",imag(y2.f[i])
                    )
    end
    close(fp)

    axpby!(1,y,-1,y2)
    println(sqrt(real(dot(y2,y2)))/sqrt(real(dot(y,y))))

    return

    a = zero(x)
    for i= 1:length(e2)
        for j= 1:length(e2)
            a.f[i] += conj(v[j,i])*b.f[j]
        end
    end

    y2 = zero(x)
    zseedin = 1+0.05im

    α = 0.1
    radius = 0.99
    origin = 1

    vec_z = zeros(ComplexF64,Nq)
    vec_w = zeros(ComplexF64,Nq)
    for j in 1:Nq
        θj = (2π/Nq)*(j-1/2) # π/Nq   2π-π/Nq
        vec_z[j] = origin + radius*(cos(θj)+im*α*sin(θj))
        vec_w[j] = α*cos(θj)+im*sin(θj)
    end

    y_z = zero(x)
    a_z = zero(x)

    for i=1:length(e2)
        for j in 1:Nq
            coeff = radius*vec_w[j]/Nq
            #coeff = radius*vec_w[j]*sign(vec_z[j])/Nq
            a_z.f[i] += coeff*a.f[i]/(vec_z[j]-e2[i])
        end
    end

    for α= 1:length(e2)
        for i= 1:length(e2)
            y_z.f[α] += v[α,i]*a_z.f[i]
        end
    end

    radius = 0.99
    origin = -1

    vec_z = zeros(ComplexF64,Nq)
    vec_w = zeros(ComplexF64,Nq)
    for j in 1:Nq
        θj = (2π/Nq)*(j-1/2) # π/Nq   2π-π/Nq
        vec_z[j] = origin + radius*(cos(θj)+im*α*sin(θj))
        vec_w[j] = α*cos(θj)+im*sin(θj)
    end

    a_z = zero(x)

    for i=1:length(e2)
        for j in 1:Nq
            #coeff = radius*vec_w[j]*sign(vec_z[j])/Nq
            coeff = radius*vec_w[j]/Nq
            a_z.f[i] += coeff*a.f[i]/(vec_z[j]-e2[i])
        end
    end

    for α= 1:length(e2)
        for i= 1:length(e2)
            y_z.f[α] += -v[α,i]*a_z.f[i]
        end
    end

#=
    fp = open("signHx_z.txt","w")
    for i= 1:length(e2)
        println(fp,real(y.f[i]),"\t",imag(y.f[i]),"\t",real(y_z.f[i]),"\t",imag(y_z.f[i]))
        println(real(y.f[i]),"\t",imag(y.f[i]),"\t",real(y_z.f[i]),"\t",imag(y_z.f[i]))

    end
    close(fp)

    axpby!(1,y,-1,y_z)
    println(sqrt(real(dot(y_z,y_z)))/sqrt(real(dot(y,y))))
    return
    =#

    #vec_x = shiftedbicg_Frommer2003_G_seed(σi,g5D,b,eps =1e-20,zseedin = zseedin) #(z I - A) x = b
    #(z I - A) x = b -> (-z I + A)(-x) = b
    vec_x = shiftedbicg_Frommer2003_G_seed(vec_z,g5D,b,eps =1e-20,zseedin = zseedin) #(z I - A) x = b

    y3 = zero(x)

    for i=1:length(σi)
        println("i = $i σi = $(vec_z[i])")
        mul!(y3,g5D,vec_x[i])
        axpby!(vec_z[i],vec_x[i], -1, y3) #Overwrite Y with X*a + Y*b
        #axpby!(σi[i],vec_x[i], -1, y3) #Overwrite Y with X*a + Y*b

        axpby!(-1,b, 1, y3) #Overwrite Y with X*a + Y*b
        #add_fermion!(y,-σi[i],vec_x[i])
        println(dot(y3,y3))
    end

    #=
    for j = 1:div(Nq,2)
        axpy!(-sign(σi[j])*dr/(2π*im),vec_x[j],y2)
        axpy!(sign(σi[j+div(Nq,2)])*dr/(2π*im),vec_x[j+div(Nq,2)],y2)
    end
    =#

    for j in 1:Nq
        coeff = radius*vec_w[j]*sign(vec_z[j])/Nq
        axpy!(coeff,vec_x[j],y2)
        #vec_Sk[:,ik] += ρ*vec_w[j]*vec_z[j]^k*vec_y[:,j]/Nq
    end

    fp = open("signHx.txt","w")

    for i= 1:length(e2)

        #println(y.f[i],"\t",y2.f[i])
        println(fp,real(y.f[i]),"\t",imag(y.f[i]),"\t",
                    real(y2.f[i]),"\t",imag(y2.f[i]),"\t",
                    real(y_z.f[i]),"\t",imag(y_z.f[i])
                    )
    end
    close(fp)

    axpby!(1,y,-1,y2)
    println(sqrt(real(dot(y2,y2)))/sqrt(real(dot(y,y))))

    #=
    fp = open("signHx_z.txt","w")
    for i= 1:length(e2)
        println(fp,real(y.f[i]),"\t",imag(y.f[i]),"\t",real(y_z.f[i]),"\t",imag(y_z.f[i]))
        println(real(y.f[i]),"\t",imag(y.f[i]),"\t",real(y_z.f[i]),"\t",imag(y_z.f[i]))

    end
    close(fp)
    =#

    axpby!(1,y,-1,y_z)
    println(sqrt(real(dot(y_z,y_z)))/sqrt(real(dot(y,y))))
    return

    return

    for j = 1:div(Nq,2)
        axpy!(-sign(σi[j])*dr/(2π*im),vec_x[j],y2)
        axpy!(sign(σi[j+div(Nq,2)])*dr/(2π*im),vec_x[j+div(Nq,2)],y2)
    end
    #println(dot(y,y2))

    for i= 1:length(e2)
        println(y.f[i],"\t",y2.f[i])
    end
    println(dot(y,y2))
    return

    radius = 0.1
    origin = 0.75
    result = eigensystem(g5D,x,radius,origin,Z2_distribution_fermion!,
                uniform_distribution_fermion!,α=1,epsCG=1e-20,κ=1.5)
    println(result[1])

    return

    radius = 0.03
    origin = 1
    result = eigensystem(D(U),x,radius,origin,Z2_distribution_fermion!,
                uniform_distribution_fermion!,α=1,epsCG=1e-20,κ=1.5)
    println(result[1])

    b = similar(x)
    y = similar(x)
    rtin = similar(x)

    for i in 1:length(b.f)
        b.f[i] = rand([-1,1])
        rtin.f[i] = rand([-1,1])
    end
    #b[1,1,1,1,1,1] = 1
    #b[1,1,1,1,1,2] = 1

    #shiftedbicgstab(x,D,b)
    #mul!(y,D,x)
    #axpby!(-1,b, 1, y)
    #println(dot(y,y))
    #return

    N = 10
    Nq = N
    zseed = 0.1
    γ = 1
    ρ = 0.1
    α= 1

    σi = zeros(ComplexF64,Nq)

    for j = 1:Nq
        θ = (2π/Nq)*(j-1/2)
        σi[j] = γ + ρ*(cos(θ)+im*sin(θ))
    end

    #σi = rand(ComplexF64,N)

    vec_x = shiftedbicg_Frommer2003_seed(σi,D(U),b,zseed=σi[1],eps =1e-20) #(zI - D)*x = b
    #vec_x = shiftedbicg_Frommer2003(σi,D(U),b)
    for i=1:length(σi)
        println("i = $i σi = $(σi[i])")
        mul!(y,D,vec_x[i])
        axpby!(σi[i],vec_x[i], 1, y) #Overwrite Y with X*a + Y*b
        axpby!(-1,b, 1, y) #Overwrite Y with X*a + Y*b
        #add_fermion!(y,-σi[i],vec_x[i])
        println(dot(y,y))
    end
    return

    #=
    for i=1:length(σi)
        println("i = $i σi = $(σi[i])")
        mul!(y,D,vec_x[i])
        axpby!(σi[i],vec_x[i], -1, y) #Overwrite Y with X*a + Y*b
        axpby!(-1,b, 1, y) #Overwrite Y with X*a + Y*b
        #add_fermion!(y,-σi[i],vec_x[i])
        println(dot(y,y))
    end

    return

    x,vec_x = shiftedbicg_inSS(σi,D,b,eps=1e-28,zseedin=zseed)
    mul!(y,D,x)
    axpy!(zseed,x,y)
    axpby!(-1,b, 1, y)
    println(dot(y,y))
    #println(dot(x,vec_x[1]))
    for i=1:length(σi)
        println("i = $i σi = $(σi[i])")
        mul!(y,D,vec_x[i])
        axpby!(σi[i],vec_x[i], 1, y) #Overwrite Y with X*a + Y*b
        axpby!(-1,b, 1, y) #Overwrite Y with X*a + Y*b
        #add_fermion!(y,-σi[i],vec_x[i])
        println(dot(y,y))
    end
    return

    =#

    result = eigensystem(D(U),x,length(x),0.05,1,α=1,epsCG=1e-20,κ=1.5)
    println(result[1])
    e = result[1]

    fp = open("eigenvalues_SS_wilson.txt","w")
    for i = 1:length(e)
        println(fp,real(e[i]),"\t",imag(e[i]))
        println(real(e[i]),"\t",imag(e[i]))
    end
    close(fp)

    Dmat = construct_sparsematrix(D(U))
    e2,v = eigen(Matrix(Dmat))
    fp = open("eigenvalues_wilson.txt","w")
    for i = 1:length(e2)
        println(fp,real(e2[i]),"\t",imag(e2[i]))
        println(real(e2[i]),"\t",imag(e2[i]))
    end
    close(fp)

    return

    parameters_action = Dict()
    fermi_action = FermiAction(D,parameters_action)
    #gauss_sampling_in_action!(x,U,fermi_action)
    #println("Sfold = ", dot(x,x))
    y = similar(x)

    MDtest!(gauge_action,U,Dim,fermi_action,x,y)

end

function gauss_distribution(nv)
    variance = 1
    nvh = div(nv,2)
    granf = zeros(Float64,nv)
    for i=1:nvh
        rho = sqrt(-2*log(rand())*variance)
        theta = 2pi*rand()
        granf[i] = rho*cos(theta)
        granf[i+nvh] = rho*sin(theta)
    end
    if 2*nvh == nv
        return granf
    end

    granf[nv] = sqrt(-2*log(rand())*variance) * cos(2pi*urand())
    return granf
end

test1()


text:
module ExampleTests
include("test_nooptask.jl")
end # module Examples

max_forks_repo_name: seung-lab/SimpleTasks


text:
# Computes transit timing variations to linear order
# in eccentricity.  Please cite Agol & Deck (2015) if
# you make use of this in published research.

module TTVFaster

VERSION < v"0.4-dev" && using Docile

export Planet_plane_hk, compute_ttv!
export Planet_plane                   # Deprecated, only exported to make the error message work

#include("ttv_coeff.jl")    # TODO: Any reason not to delete this?
include("ttv_succinct.jl")

include("laplace_coefficients_initialize.jl")
laplace_coefficients_initialize(jmax::Integer,alpha::Number) = LaplaceCoefficients.initialize(jmax,alpha)
export laplace_coefficients_initialize

immutable Planet_plane
# Parameters of a planet in a plane-parallel system
  # Mass ratio of the planet to the star:
  mass_ratio :: Float64
  # Initial time of transit:
  period   :: Float64
  trans0   :: Float64
  eccen    :: Float64
# longitude of periastron measured from line of sight, in radians:
  omega    :: Float64
end

type Planet_plane_hk{T<:Number} # Parameters of a planet in a plane-parallel system
  # Mass ratio of the planet to the star:
  mass_ratio :: T
  # Initial time of transit:
  period   :: T
  trans0   :: T
  # e times cos or sin of longitude of periastron measured from line of sight, in radians:
  ecosw    :: T
  esinw    :: T
end

"""
# Error message to explain to anyone who tries to use the old version
"""
function compute_ttv!(jmax::Integer,p1::Planet_plane,p2::Planet_plane,time1::Vector,time2::Vector,ttv1::Vector,ttv2::Vector,f1::Array,f2::Array,b::Array,alpha0::Number,b0::Array)
  error("The Planet_plane data structure has been deprecated in favor of Planet_plane_hk")
end

"""
# Computes transit-timing variations to linear order in
# eccentricity for non-resonant, plane-parallel planets.
# Input:
#   jmax:  Maximum j over which to sum the TTV calculation for both planets
#     p1:  Planet type for inner planet
#     p2:  Planet type for outer planet
#  time1:  Transit times for inner planet
#  time2:  Transit times for outer planet
# alpha0:  Initial alpha for Taylor expansion of coefficients
#     b0:  Laplace coefficients and derivatives for use in Taylor expansion
# Output:
#   ttv1: TTVs of the inner planet
#   ttv2: TTVs of the outer planet
#     f1: TTV coefficients for inner planet
#     f2: TTV coefficients for outer planet
#      b: Laplace coefficients (& derivatives) for outer planet
"""
function compute_ttv!(jmax::Integer,p1::Planet_plane_hk,p2::Planet_plane_hk,time1::Vector,time2::Vector,ttv1::Vector,ttv2::Vector,f1::Array,f2::Array,b::Array,alpha0::Number,b0::Array)

# Compute the semi-major axis ratio of the planets:
# println(p1.period,p2.period)
const alpha = (p1.period/p2.period)^(2//3)  # Julia supports rational numbers!
# Number of times:
const ntime1 = length(time1)
const ntime2 = length(time2)
# Compute the coefficients:
ttv_succinct!(jmax+1,alpha,f1,f2,b,alpha0,b0)  # I need to compute coefficients one higher than jmax
# Compute TTVs for inner planet (equation 33):
# Compute since of \pomegas:
e1 = sqrt(p1.esinw*p1.esinw+p1.ecosw*p1.ecosw)
e2 = sqrt(p2.esinw*p2.esinw+p2.ecosw*p2.ecosw)
sin1om=p1.esinw/e1
sin2om=p2.esinw/e2
cos1om=p1.ecosw/e1
cos2om=p2.ecosw/e2
# Compute mean motions:
n1=2pi/p1.period
n2=2pi/p2.period
# Compute initial longitudes:
lam10=-n1*p1.trans0 + 2*p1.esinw # 2*p1.eccen*sin1om
lam20=-n2*p2.trans0 + 2*p2.esinw # 2*p2.eccen*sin2om
@inbounds for i=1:ntime1
# Compute the longitudes of the planets at times of transit of planet 1 (equation 49):
  lam11 = n1*time1[i]+lam10
  lam21 = n2*time1[i]+lam20
  psi1  = lam11-lam21 # Compute difference in longitudes at times of transit of planet 1
  sinpsi1=sin(psi1)
  cospsi1=cos(psi1)
  sinlam11 = sin(lam11)
  coslam11 = cos(lam11)
  sinlam1om1=sinlam11*cos1om-coslam11*sin1om
  coslam1om1=coslam11*cos1om+sinlam11*sin1om
  sinlam1om2=sinlam11*cos2om-coslam11*sin2om
  coslam1om2=coslam11*cos2om+sinlam11*sin2om
  ttv1[i]=zero(p1.period) #0.0
  sinjm1psi1=zero(p1.period) #0.0
  cosjm1psi1=one(p1.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi1=sinjm1psi1*cospsi1+cosjm1psi1*sinpsi1
    cosjpsi1=cosjm1psi1*cospsi1-sinjm1psi1*sinpsi1
    ttv1[i] += f1[j+1,1]*sinjpsi1
    ttv1[i] += f1[j+1,2]*e1*(sinjpsi1*coslam1om1-cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j+1,3]*e1*(sinjpsi1*coslam1om1+cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j  ,4]*e2*(sinjpsi1*coslam1om2-cosjpsi1*sinlam1om2)
    ttv1[i] += f1[j+2,5]*e2*(sinjpsi1*coslam1om2+cosjpsi1*sinlam1om2)
    sinjm1psi1=sinjpsi1
    cosjm1psi1=cosjpsi1
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv1[i] = ttv1[i]*p1.period*p2.mass_ratio/(2pi)
end
# Compute TTVs for outer planet (equation 33):
@inbounds for i=1:ntime2
# Compute the longitudes of the planets at times of transit of planet 2:
  lam12 = n1*time2[i]+lam10
  lam22 = n2*time2[i]+lam20
  sinlam22 = sin(lam22)
  coslam22 = cos(lam22)
  psi2  = lam12-lam22 # Compute difference in longitudes at times of transit of planet 2
  sinpsi2=sin(psi2)
  cospsi2=cos(psi2)
  sinlam2om1=sinlam22*cos1om-coslam22*sin1om
  coslam2om1=coslam22*cos1om+sinlam22*sin1om
  sinlam2om2=sinlam22*cos2om-coslam22*sin2om
  coslam2om2=coslam22*cos2om+sinlam22*sin2om
  ttv2[i]=zero(p2.period) #0.0
  sinjm1psi2=zero(p2.period) #0.0
  cosjm1psi2=one(p2.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi2=sinjm1psi2*cospsi2+cosjm1psi2*sinpsi2
    cosjpsi2=cosjm1psi2*cospsi2-sinjm1psi2*sinpsi2
    ttv2[i] += f2[j+1,1]*sinjpsi2
    ttv2[i] += f2[j+1,2]*e2*(sinjpsi2*coslam2om2-cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+1,3]*e2*(sinjpsi2*coslam2om2+cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+2,4]*e1*(sinjpsi2*coslam2om1-cosjpsi2*sinlam2om1)
    ttv2[i] += f2[j  ,5]*e1*(sinjpsi2*coslam2om1+cosjpsi2*sinlam2om1)
    sinjm1psi2=sinjpsi2
    cosjm1psi2=cosjpsi2
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv2[i] = ttv2[i]*p2.period*p1.mass_ratio/(2pi)
end
# Finished!
return
end  # compute_ttv!

end # module

max_stars_repo_name: ericagol/TTVFaster

max_issues_repo_name: ericagol/TTVFaster

max_forks_repo_name: ericagol/TTVFaster


text:
## Logical operators for (Sym,Sym)

# http://docs.sympy.org/dev/_modules/sympy/core/relational.html
logical_sympy_methods = (:Eq, :Ne, :Lt, :Le, :Gt, :Ge)

## these are defined for Reals, not Sym, Real...
## Necessary to lambdify Indicators as of v0.6.0-dev
## not exported
const relational_sympy_values = (:GreaterThan, :LessThan,
                                 :StrictGreaterThan, :StrictLessThan,
                                 :Equality, :Unequality)
for meth in relational_sympy_values
    meth_name = string(meth)
    @eval begin
        @doc """
`$($meth_name)`: a SymPy function. [cf.](http://docs.sympy.org/dev/_modules/sympy/core/relational.html)
""" ->
        ($meth)(a::Real, b::Real) = sympy_meth($meth_name,a, b)
    end
#    eval(Expr(:export, meth))
end

## XXX Experimental! Not sure these are such a good idea ...
## but used with piecewise
Base.:&(x::Sym, y::Sym) = PyCall.pycall(PyObject(x)["__and__"], Sym, y)
Base.:|(x::Sym, y::Sym) =  PyCall.pycall(PyObject(x)["__or__"], Sym, y)
!(x::Sym)         =       PyCall.pycall(x.x["__invert__"], Sym)::Sym

## use ∨, ∧, ¬ for |,&,! (\vee<tab>, \wedge<tab>, \neg<tab>)
∨(x::Sym, y::Sym) = x | y
∧(x::Sym, y::Sym) = x & y
¬(x::Sym) = !x

## In SymPy, symbolic equations are not represented by `=` or `==`
## rather ther function `Eq` is used. Here we use the unicode
## `\Equal<tab>` for an infix operator. There are also unicode values to represent analogs of `<`, `<=`, `>=`. `>`. These are

## * `<`  is `\ll<tab>`
## * `<=` is `\leqq<tab>
## * `==` is `\Equal<tab>`
## * `>=` is `\geqq<tab>`
## * `>`  is `\gg<tab>`

## Instead we have:
## We use unicode for visual appeal of infix operators, but the Lt, Le, Eq, Ge, Gt are the proper way:

"This is `\\ll<tab>` mapped as an infix operator to `Lt`"
(≪)(a::Sym, b::Sym) = Lt(a,b)  # \ll<tab>
(≪)(a::Sym, b::Number) = Lt(a,Sym(b)) # \ll<tab>
(≪)(a::Number, b::Sym) = Lt(Sym(a),b) # \ll<tab>

## could just do this, but it would interfere with other uses outside of SymPy
## (≪)(a::Number, b::Number) = Lt(promote(a,b)...)  # \ll<tab>

"This is `\\leqq<tab>` mapped as an infix operator to `Le`"
(≦)(a::Sym, b::Sym) = Le(a,b)   # \ll<tab>
(≦)(a::Sym, b::Number) = Le(a,Sym(b))  # \ll<tab>
(≦)(a::Number, b::Sym) = Le(Sym(a),b)   # \ll<tab>

"This is `\\gg<tab>` mapped as an infix operator to `Gt`"
(≫)(a::Sym, b::Sym) = Gt(a,b) |> asBool
(≫)(a::Sym, b::Number) = Gt(a,Sym(b))
(≫)(a::Number, b::Sym) = Gt(Sym(a),b)

"This is `\\geqq<tab>` mapped as an infix operator to `Ge`"
(≧)(a::Sym, b::Sym) = Ge(a,b) |> asBool
(≧)(a::Sym, b::Number) = Ge(a,Sym(b))
(≧)(a::Number, b::Sym) = Ge(Sym(a),b)

"For infix `Eq` one can use \\Equal<tab> unicode operator"
(⩵)(a::Sym, b::Sym) = Eq(a,b)  # \Equal<tab>
(⩵)(a::Sym, b::Number) = Eq(a,Sym(b))  # \Equal<tab>
(⩵)(a::Number, b::Sym) = Eq(Sym(a),b)  # \Equal<tab>

export ≪,≦,⩵,≧,≫

"For hashing, we use equality at the python level."
==(x::Sym, y::Sym) = PyObject(x) == PyObject(y)  # but isequal is not the same!

## Question: what to do with comparisons between Sym and Sym which *can* be answered true or false?
## Here we throw an error when that is not the case in hopes that generic code can be called
## when decisions are possible.
## This is experimental

function asBool(x::Sym)
    x == SympyTRUE && return true
    x == !SympyTRUE && return false
    funcname(x) == "Equality" && return ==(args(x)...)

    convert(Bool, x)
#    throw(DomainError("Non Boolean value"))
end

Base.isless(a::Sym, b::Sym) = asBool(Lt(a,b))
Base.isless(a::Sym, b::Number) = asBool(Lt(promote(a,b)...))
Base.isless(a::Number, b::Sym) = asBool(Lt(promote(a,b)...))

Base.isequal(a::Sym, b::Sym) = asBool(Eq(a,b))
Base.isequal(a::Sym, b::Number) = asBool(Eq(promote(a,b)...))
Base.isequal(a::Number, b::Sym) = asBool(Eq(promote(a,b)...))

function !=(x::Sym, y::T)  where {T <: Real}
    try
        x = convert(Float64, x)
        x != y
    catch
        true
    end
end
function !=(x::Sym, y::T) where {T <: Complex}
    try
        x = complex(x)
        x != y
    catch
        true
    end
end

function init_logical()
    global const SympyTRUE = sympy_meth(:Lt, 0,1)
end


text:
const KMORAL = 5
const NSOC   = 50
const ρDEF   = 0.3
const λDEF	 = 2.0
const εDEF   = 0.1
const μDEF   = 0.
const s2DEF  = 4.0
const βDEF   = 15.0


text:
# distritbutions
@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_var(m_in2)
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - min2, vout + vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_cov(m_in2)
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - min2, vout + vin2)
end

# PointMass
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::PointMass) = PointMass(mean(m_out) - mean(m_in2))

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::NormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return NormalMeanPrecision(mean(m_out) - min2, win2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::NormalMeanPrecision, m_in2::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return NormalMeanPrecision(mout - mean(m_in2), wout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MvNormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return MvNormalMeanPrecision(mean(m_out) - min2, win2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MvNormalMeanPrecision, m_in2::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return MvNormalMeanPrecision(mout - mean(m_in2), wout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::UnivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(mean(m_out) - min2, vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(mout - mean(m_in2), vout)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(mean(m_out) - min2, vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(mout - mean(m_in2), vout)
end


text:
function variable_converter_filter_voltage(pm::AbstractBFModel; kwargs...)
    variable_converter_filter_voltage_magnitude_sqr(pm; kwargs...)
    variable_conv_transformer_current_sqr(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::AbstractBFModel; kwargs...)
    variable_converter_internal_voltage_magnitude_sqr(pm; kwargs...)
    variable_conv_reactor_current_sqr(pm; kwargs...)
end
"""
Converter transformer constraints
```
p_tf_fr + ptf_to ==  rtf*itf
q_tf_fr + qtf_to ==  xtf*itf
p_tf_fr^2 + qtf_fr^2 <= w/tm^2 * itf
wf == w/tm^2 -2*(rtf*ptf_fr + xtf*qtf_fr) + (rtf^2 + xtf^2)*itf
```
"""
function constraint_conv_transformer(pm::AbstractBFQPModel, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    w = PowerModels.var(pm, n, cnd, :w, acbus)
    itf = PowerModels.var(pm, n, cnd, :itf_sq, i)
    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)

    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    qtf_fr = PowerModels.var(pm, n, cnd, :qconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    if transformer
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model,   ptf_fr + ptf_to ==  rtf*itf)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model,   qtf_fr + qtf_to ==  xtf*itf)
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model,   ptf_fr^2 + qtf_fr^2 <= w/tm^2 * itf)
        PowerModels.con(pm, n, cnd, :conv_tf_q_to)[i] = @constraint(pm.model,   wf == w/tm^2 -2*(rtf*ptf_fr + xtf*qtf_fr) + (rtf^2 + xtf^2)*itf)
    else
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, ptf_fr + ptf_to == 0)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model, qtf_fr + qtf_to == 0)
        @constraint(pm.model, wf == w )
    end
end

function constraint_conv_transformer(pm::AbstractBFConicModel, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    w = PowerModels.var(pm, n, cnd, :w, acbus)
    itf = PowerModels.var(pm, n, cnd, :itf_sq, i)
    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)

    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    qtf_fr = PowerModels.var(pm, n, cnd, :qconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    if transformer
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model,   ptf_fr + ptf_to ==  rtf*itf)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model,   qtf_fr + qtf_to ==  xtf*itf)
        # PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model,   norm([2*ptf_fr; 2*qtf_fr; w/tm^2 - itf]) <= w/tm^2 + itf)
        # PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model,   norm([2*ptf_fr; 2*qtf_fr; w/tm - itf/tm]) <= w/tm + itf/tm)
        # PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model,   [2*ptf_fr, 2*qtf_fr, w/tm - itf/tm, w/tm + itf/tm] in JuMP.RotatedSecondOrderCone())
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model,   [w/(sqrt(2)*tm), itf/(sqrt(2)*tm), ptf_fr, qtf_fr] in JuMP.RotatedSecondOrderCone())
        PowerModels.con(pm, n, cnd, :conv_tf_q_to)[i] = @constraint(pm.model,   wf == w/tm^2 -2*(rtf*ptf_fr + xtf*qtf_fr) + (rtf^2 + xtf^2)*itf)
    else
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, ptf_fr + ptf_to == 0)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model, qtf_fr + qtf_to == 0)
        @constraint(pm.model, wf == w/tm^2 )
    end
end

"""
Converter reactor constraints
```
p_pr_fr + ppr_to == rc*ipr
q_pr_fr + qpr_to == xc*ipr
p_pr_fr^2 + qpr_fr^2 <= wf * ipr
wc == wf -2*(rc*ppr_fr + xc*qpr_fr) + (rc^2 + xc^2)*ipr
```
"""
function constraint_conv_reactor(pm::AbstractBFQPModel, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)
    ipr = PowerModels.var(pm, n, cnd, :irc_sq, i)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    ppr_to = -PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qpr_to = -PowerModels.var(pm, n, cnd, :qconv_ac, i)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)

    if reactor
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == rc*ipr)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == xc*ipr)
        @constraint(pm.model, ppr_fr^2 + qpr_fr^2 <= wf * ipr)
        @constraint(pm.model, wc == wf -2*(rc*ppr_fr + xc*qpr_fr) + (rc^2 + xc^2)*ipr)

    else
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == 0)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == 0)
        @constraint(pm.model, wc == wf)
    end
end

function constraint_conv_reactor(pm::AbstractBFConicModel, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    wf = PowerModels.var(pm, n, cnd, :wf_ac, i)
    ipr = PowerModels.var(pm, n, cnd, :irc_sq, i)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    ppr_to = -PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qpr_to = -PowerModels.var(pm, n, cnd, :qconv_ac, i)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)

    if reactor
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == rc*ipr)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == xc*ipr)
        # @constraint(pm.model, norm([2*ppr_fr; 2*qpr_fr; wf - ipr]) <= wf + ipr)
        # @constraint(pm.model, norm([2*ppr_fr; 2*qpr_fr; wf - ipr]) <= wf + ipr)
        # @constraint(pm.model, [2*ppr_fr, 2*qpr_fr, wf - ipr, wf + ipr] in JuMP.RotatedSecondOrderCone())
        @constraint(pm.model, [wf/sqrt(2), ipr/sqrt(2), ppr_fr, qpr_fr] in JuMP.RotatedSecondOrderCone())
        @constraint(pm.model, wc == wf -2*(rc*ppr_fr + xc*qpr_fr) + (rc^2 + xc^2)*ipr)

    else
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == 0)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == 0)
        @constraint(pm.model, wc == wf)
    end
end
"""
Links converter power & current
```
pconv_ac[i]^2 + pconv_dc[i]^2 <= wc[i] * iconv_ac_sq[i]
pconv_ac[i]^2 + pconv_dc[i]^2 <= (Umax)^2 * (iconv_ac[i])^2
```
"""
function constraint_converter_current(pm::AbstractBFQPModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)
    iconv_sq = PowerModels.var(pm, n, cnd, :iconv_ac_sq, i)

    PowerModels.con(pm, n, cnd, :conv_i)[i] = @constraint(pm.model,      pconv_ac^2 + qconv_ac^2 <=  wc * iconv_sq)
    PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @constraint(pm.model, pconv_ac^2 + qconv_ac^2 <= (Umax)^2 * iconv^2)
    @constraint(pm.model, iconv^2 <= iconv_sq)
    @constraint(pm.model, iconv_sq <= iconv*Imax)
end

function constraint_converter_current(pm::AbstractBFConicModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    wc = PowerModels.var(pm, n, cnd, :wc_ac, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)
    iconv_sq = PowerModels.var(pm, n, cnd, :iconv_ac_sq, i)

    # Normal form
    # PowerModels.con(pm, n, cnd, :conv_i)[i] = @constraint(pm.model, [2*pconv_ac, 2*qconv_ac,  wc - iconv_sq, wc + iconv_sq] in JuMP.RotatedSecondOrderCone())
    PowerModels.con(pm, n, cnd, :conv_i)[i] = @constraint(pm.model, [wc/sqrt(2), iconv_sq/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @constraint(pm.model, [Umax * iconv/sqrt(2), Umax * iconv/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    # PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @constraint(pm.model, norm([2*pconv_ac; 2*qconv_ac; Umax^2 - iconv_sq]) <= (Umax)^2 + iconv_sq)
    # @constraint(pm.model, [2*iconv, 1 - iconv_sq, 1 + iconv_sq] in JuMP.RotatedSecondOrderCone())
    @constraint(pm.model, [iconv_sq, iconv/sqrt(2), iconv/sqrt(2)] in JuMP.RotatedSecondOrderCone())
    @constraint(pm.model, iconv_sq <= iconv*Imax)
end


text:
##########################
# LU Factorization w/ PP #
##########################
function lufact(C::KronProdMat)
    Fs = [lufact(term) for term in terms(C)]
    ls,us,ps = zip([(F[:L],F[:U],F[:P]) for F in Fs]...)
    #note the the following only returns permutation matrices, as opposed to pivot indexes
    return [:L=>KronProdMat(ls...),:U=>KronProdMat(us...),:P=>KronProdMat(ps...)]
end

##########################
# Cholesky Factorization #
##########################
chol(C::KronProdMat) = KronProdMat([chol(term) for term in terms(C)])

###################
#QR Factorization #
###################
function qrfact(C::KronProdMat)
    qs,rs = zip([qr(term) for term in terms(C)])
    return [:Q=>KronProdMat(qs...),:R=>KronProdMat(rs...)]
end

######################
# Eigendecomposition #
######################
function eigfact(C::KronProdMat)
    issquare(C) || throw(DimensionMismatch("matrix not square"))
    all(issquare,terms(C)) || throw(DimensionMismatch("currently only supported for square inner and outer"))
    λs,us = zip([eig(term) for term in terms(C)]...)
    return [:vectors=>KronProdMat(us...),:values=>KronProdVec(λs...)]
end

eigvals(C::KronProdMat) = eigfact(C)[:values]

#######
# SVD #
#######
function svdfact(C::KronProdMat)
    us,σs,vs = zip([svd(term) for term in terms(C)]...)
    return [:U=>KronProdMat(us...),:S=>KronProdMat(σs...),:V=>KronProdMat(vs...)]
end

function svd(C::KronProdMat)
    F = svdfact(C)
    return (F[:U],F[:S],F[:V])
end

svdvals{T}(C::KronProdMat{T}) = KronProdVec([svdvals(term) for term in terms(C)]...)


text:
using JSON
include("../bike/bike.jl")
include("../bike/bike-initializer.jl")
include("../rider/rider.jl")

function initialize_simulation(number_of_bikes, number_of_seconds) :: Simulation

    bikes = Vector{Bike}(undef, number_of_bikes)
    for i in 1:number_of_bikes
        bikes[i] = initialize_random_bike()
    end
    timed_bike_positions = zeros(Int64, (number_of_seconds, 2, number_of_bikes))
    return Simulation(bikes, 1, timed_bike_positions, number_of_seconds)

end

function run_simulation_step(simulation :: Simulation)

    for bike_index in 1:length(simulation.bikes)
        bike = simulation.bikes[bike_index]
        move_bike(bike, simulation)
        simulation.timed_bike_positions[simulation.current_time, 1, bike_index] = bike.x_position
        simulation.timed_bike_positions[simulation.current_time, 2, bike_index] = bike.y_position

    end
    simulation.current_time += 1
end

function print_results(simulation :: Simulation)
    for t in 1:simulation.simulation_run_time
        println("time: $(t)")
        for bi in 1:length(simulation.bikes)
            bike = simulation.bikes[bi]
            println("bike $(bike.bike_attributes.make) $(bike.bike_attributes.model): x: $(simulation.timed_bike_positions[t, 1, bi]) y: $(simulation.timed_bike_positions[t, 2, bi])")
        end
    end
end

function run_simulation(number_of_bikes, number_of_seconds)

    simulation = initialize_simulation(number_of_bikes, number_of_seconds)

    while simulation.current_time <= number_of_seconds
        run_simulation_step(simulation)
    end
    println(JSON.json(simulation))
    print_results(simulation)
end


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
end
description = """
Showing how one can color mesh vertices.
"""

# one could also create a GLNormalAttributeMesh manually:
sphere = Sphere{Float32}(Point3f0(0), 2f0)
# decompose decomposes a mesh or geometry into the primitive of the first argument
vertices = decompose(Point3f0, sphere, 50)
faces = decompose(GLTriangle, sphere, 50)
# create some colors for each vertex, must be Float32
colors = RGBA{Float32}[RGBA{Float32}(rand(), rand(), rand(), 1.) for i=1:length(vertices)]

sphere_mesh = GLNormalVertexcolorMesh(
    vertices = vertices, faces = faces,
    color = colors
)
_view(visualize(sphere_mesh), window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
"""
    circuit_extraction(zxg::ZXGraph)

Extract circuit from a graph-like ZX-diagram.
"""
function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
    nzxg = copy(zxg)
    nbits = nqubits(nzxg)
    gads = Set{T}()
    for v in spiders(nzxg)
        if spider_type(nzxg, v) == SpiderType.Z && degree(nzxg, v) == 1
            push!(gads, v, neighbors(zxg, v)[1])
        end
    end

    Outs = get_outputs(nzxg)
    Ins = get_inputs(nzxg)
    if nbits == 0
        nbits = length(Outs)
    end
    cir = Chain()
    if length(Outs) != length(Ins)
        return cir
    end
    for v1 in Ins
        @inbounds v2 = neighbors(nzxg, v1)[1]
        if !is_hadamard(nzxg, v1, v2)
            insert_spider!(nzxg, v1, v2)
        end
    end
    @inbounds frontier = [neighbors(nzxg, v)[1] for v in Outs]
    qubit_map = Dict(zip(frontier, 1:nbits))

    for i = 1:nbits
        @inbounds w = neighbors(nzxg, Outs[i])[1]
        @inbounds if is_hadamard(nzxg, w, Outs[i])
            pushfirst_gate!(cir, Val{:H}(), i)
        end
        if phase(nzxg, w) != 0
            pushfirst_gate!(cir, Val{:Rz}(), i, phase(nzxg, w))
            set_phase!(nzxg, w, zero(P)) end
        @inbounds rem_edge!(nzxg, w, Outs[i])
    end
    for i = 1:nbits
        for j = i+1:nbits
            @inbounds if has_edge(nzxg, frontier[i], frontier[j])
                if is_hadamard(nzxg, frontier[i], frontier[j])
                    pushfirst_gate!(cir, Val{:CZ}(), i, j)
                    rem_edge!(nzxg, frontier[i], frontier[j])
                end
            end
        end
    end

    old_frontier = copy(frontier)
    max_iter = 1000
    current_iter = 1
    while !isempty(frontier)
        update_frontier!(nzxg, gads, frontier, qubit_map, cir)
        if frontier != old_frontier
            old_frontier = copy(frontier)
            current_iter = 1
        else
            current_iter += 1
            if current_iter > max_iter
                error("Circuit extraction failed!")
            end
        end
    end

    frontier = T[]
    for v in Ins
        nb = neighbors(nzxg, v)
        if length(nb) > 0
            push!(frontier, nb[])
        end
    end
    sort!(frontier, by = (v->qubit_map[v]))
    M = biadjacency(nzxg, frontier, Ins)
    M, steps = gaussian_elimination(M)
    for step in steps
        if step.op == :addto
            ctrl = step.r2
            loc = step.r1
            pushfirst_gate!(cir, Val{:CNOT}(), loc, ctrl)
        elseif step.op == :swap
            q1 = step.r1
            q2 = step.r2
            pushfirst_gate!(cir, Val{:SWAP}(), q1, q2)
        end
    end

    simplify_swap!(cir)
    return cir
end

"""
    update_frontier!(zxg, frontier, qubit_map, cir)

Update frontier. This is a important step in the circuit extraction algorithm.
For more detail, please check the paper [arXiv:1902.03178](https://arxiv.org/abs/1902.03178).
"""
function update_frontier!(zxg::ZXGraph{T, P}, gads::Set{T}, frontier::Vector{T}, qubit_map::Dict{T, Int}, cir) where {T, P}
    # TODO: use inplace methods
    deleteat!(frontier, [spider_type(zxg, f) != SpiderType.Z || (degree(zxg, f)) == 0 for f in frontier])

    for i = 1:length(frontier)
        v = frontier[i]
        nb_v = neighbors(zxg, v)
        u = findfirst([u in gads for u in nb_v])
        if u !== nothing
            u = nb_v[u]
            gad_u = zero(T)
            for w in neighbors(zxg, u)
                if w in gads
                    gad_u = w
                    break
                end
            end
            rewrite!(Rule{:pivot}(), zxg, [u, gad_u, v])
            pop!(gads, u)
            pop!(gads, gad_u)
            frontier[i] = u
            qubit_map[u] = qubit_map[v]
            pushfirst_gate!(cir, Val(:H), qubit_map[u])
            delete!(qubit_map, v)
            for j = 1:length(frontier)
                for k = j+1:length(frontier)
                    if is_hadamard(zxg, frontier[j], frontier[k])
                        pushfirst_gate!(cir, Val(:CZ), qubit_map[frontier[j]], qubit_map[frontier[k]])
                        rem_edge!(zxg, frontier[j], frontier[k])
                    end
                end
            end

            return frontier
        end
    end

    SetN = Set{T}()
    for f in frontier
        union!(SetN, neighbors(zxg, f))
    end
    N = collect(SetN)

    M = biadjacency(zxg, frontier, N)
    M0, steps = gaussian_elimination(M)
    ws = T[]
    @inbounds for i = 1:length(frontier)
        if sum(M0[i,:]) == 1
            push!(ws, N[findfirst(isone, M0[i,:])])
        end
    end
    # M1 = biadjacency(zxg, frontier, ws)
    @inbounds for e in findall(M .== 1)
        if has_edge(zxg, frontier[e[1]], N[e[2]])
            rem_edge!(zxg, frontier[e[1]], N[e[2]])
        end
    end
    @inbounds for e in findall(M0 .== 1)
        add_edge!(zxg, frontier[e[1]], N[e[2]])
    end

    @inbounds for step in steps
        if step.op == :addto
            ctrl = qubit_map[frontier[step.r2]]
            loc = qubit_map[frontier[step.r1]]
            pushfirst_gate!(cir, Val{:CNOT}(), loc, ctrl)
        elseif step.op == :swap
            q1 = qubit_map[frontier[step.r1]]
            q2 = qubit_map[frontier[step.r2]]

            pushfirst_gate!(cir, Val{:SWAP}(), q1, q2)
        end
    end

    for i in 1:length(frontier)
        v = frontier[i]
        if degree(zxg, v) > 1
            continue
        end
        w = neighbors(zxg, v)[1]
        if is_hadamard(zxg, v, w)
            pushfirst_gate!(cir, Val(:H), qubit_map[v])
        end
        if spider_type(zxg, w) == SpiderType.Z
            qubit_map[w] = qubit_map[v]
            if phase(zxg, w) != 0
                pushfirst_gate!(cir, Val{:Rz}(), qubit_map[w], phase(zxg, w))
                set_phase!(zxg, w, zero(P))
            end
            rem_edge!(zxg, v, w)
        else
            rem_edge!(zxg, v, w)
            add_edge!(zxg, v, w, EdgeType.SIM)
        end
        frontier[i] = w
    end

    @inbounds for i1 = 1:length(frontier)
        for i2 = i1+1:length(frontier)
            if has_edge(zxg, frontier[i1], frontier[i2])
                pushfirst_gate!(cir, Val{:CZ}(), qubit_map[frontier[i1]],
                    qubit_map[frontier[i2]])
                rem_edge!(zxg, frontier[i1], frontier[i2])
            end
        end
    end
    return frontier
end

"""
    biadjacency(zxg, F, N)

Return the biadjacency matrix of `zxg` from vertices in `F` to vertices in `N`.
"""
function biadjacency(zxg::ZXGraph{T, P}, F::Vector{T}, N::Vector{T}) where {T, P}
    M = zeros(Int, length(F), length(N))

    for i = 1:length(F)
        for v2 in neighbors(zxg, F[i])
            if v2 in N
                M[i, findfirst(isequal(v2), N)] = 1
            end
        end
    end
    return M
end

"""
    GEStep

A struct for representing steps in the Gaussian elimination.
"""
struct GEStep
    op::Symbol
    r1::Int
    r2::Int
end

"""
    gaussian_elimination(M[, steps])

Return result and steps of Gaussian elimination of matrix `M`. Here we assume
that the elements of `M` is in binary field F_2 = {0,1}.
"""
function gaussian_elimination(M::Matrix{T}, steps::Vector{GEStep} = Vector{GEStep}(); rev = false) where {T<:Integer}
    M = copy(M)
    nr, nc = size(M)
    current_col = 1
    for i = 1:nr
        if sum(M[i,:]) == 0
            continue
        end
        while current_col <= nc
            rs = findall(!iszero, M[i:nr, current_col])
            if length(rs) > 0
                sort!(rs, by = k -> sum(M[k,:]), rev = rev)
                r0 = rs[1]
                r0 += i - 1
                r0 == i && break
                M_r0 = M[r0,:]
                M[r0,:] = M[i,:]
                M[i,:] = M_r0
                step = GEStep(:swap, r0, i)
                push!(steps, step)
                break
            else
                current_col += 1
            end
        end
        current_col > nc && break
        for j = 1:nr
            j == i && continue
            if M[j, current_col] == M[i, current_col]
                M[j,:] = M[j,:] .⊻ M[i,:]
                step = GEStep(:addto, i, j)
                push!(steps, step)
            end
        end
        current_col += 1
    end
    return M, steps
end

# """
#     reverse_gaussian_elimination(M, steps)

# Apply back the operations in `steps` to `M`.
# """
# function reverse_gaussian_elimination(M, steps)
#     for i = length(steps):-1:1
#         s = steps[i]
#         op = s.op
#         r1 = s.r1
#         r2 = s.r2
#         if op == :addto
#             M[r2,:] = M[r2,:] .⊻ M[r1,:]
#         else
#             r_temp = M[r1,:]
#             M[r1,:] = M[r2,:]
#             M[r2,:] = r_temp
#         end
#     end
#     return M
# end


text:
"""
    NormalInverseGamma{T<:Real} <: ContinuousUnivariateDistribution

Normal Inverse Gamma distribution. Adapted from the `NormalInverseGamma` in the
[ConjugatePriors.jl](https://github.com/JuliaStats/ConjugatePriors.jl/blob/master/src/normalinversegamma.jl)
package.
"""
struct NormalInverseGamma{T <: Real} <: ContinuousUnivariateDistribution
    μ::T
    v::T
    α::T # shape
    θ::T # scale
    function NormalInverseGamma(μ::Real, v::Real, α::Real, θ::Real)
    	v > zero(v) && α > zero(α) && θ > zero(θ) || error("v, α and θ must be positive")
        T = promote_type(typeof(μ), typeof(v), typeof(α), typeof(θ))
    	return new{T}(T(μ), T(v), T(α), T(θ))
    end
end

function rand(d::NormalInverseGamma, numsamples::Int)
    sig2_list = rand(InverseGamma(d.α, d.θ), numsamples)
    mu_list = [rand(Normal(d.μ, sqrt(sig2 * d.v))) for sig2 in sig2_list]
    return mu_list, sig2_list
end


text:
# The paper:
#
#  Yi Sun, F. Gomez, T. Schaul, J. Schmidhuber. A Linear Time Natural Evolution
#  Strategy for Non-Separable Functions. In Proceedings of the Genetic and
#  Evolutionary Computation Conference (GECCO), Amsterdam, 2013.
#
# proposes a variant of a Natural Evolution Strategy which can successfully
# solve very high-dimensional problems. In particular they show state-of-the-art
# performance on the 512-dimensional Rosenbrock function.
#
# Here we try our DE implementations on this problem to see what kind of
# performance we get.

using BlackBoxOptim
using BlackBoxOptim.Problems

function fitness_for_opt(problem, numDimensions, populationSize, numSteps,
  optFunc = de_rand_1_bin_radiuslimited)
  problem = BlackBoxOptim.Problems.set_numdims!(numDimensions, problem)

  ss = search_space(problem)

  pop = BlackBoxOptim.rand_individuals_lhs(ss, populationSize)

  opt = optFunc(ss; population = pop)

  println("\n$(problem.name), n = $(numdims(problem)), optimizer = $(opt.name)")

  best, fitness = BlackBoxOptim.run_optimizer_on_problem(opt, problem, numSteps)
  fitness
end

problem = "Rosenbrock"
p = BlackBoxOptim.Problems.examples[problem]

# Call once to compile everything
fitness_for_opt(p, 2, 10, 100, adaptive_de_rand_1_bin)

for( (dims, steps) in [(128, 1e6), (256, 2e6), (512, 4e6)])
  println("Num dims = ", dims, ", Num steps = ", steps)
  fitness = fitness_for_opt(p, dims, 50, steps, adaptive_de_rand_1_bin)
end