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

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

Returns a list of accounts in the organization from Organizations that are affected by the
provided event. For more information about the different types of Health events, see Event.
 Before you can call this operation, you must first enable Health to work with
Organizations. To do this, call the EnableHealthServiceAccessForOrganization operation from
your organization's management account.  This API operation uses pagination. Specify the
nextToken parameter in the next request to return more results.

# Arguments
- `event_arn`: The unique identifier for the event. The event ARN has the
  arn:aws:health:event-region::event/SERVICE/EVENT_TYPE_CODE/EVENT_TYPE_PLUS_ID  format. For
  example, an event ARN might look like the following:
  arn:aws:health:us-east-1::event/EC2/EC2_INSTANCE_RETIREMENT_SCHEDULED/EC2_INSTANCE_RETIREMEN
  T_SCHEDULED_ABC123-DEF456

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_affected_accounts_for_organization(
    eventArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeAffectedAccountsForOrganization",
        Dict{String,Any}("eventArn" => eventArn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_affected_accounts_for_organization(
    eventArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return health(
        "DescribeAffectedAccountsForOrganization",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("eventArn" => eventArn), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns a list of entities that have been affected by the specified events, based on the
specified filter criteria. Entities can refer to individual customer resources, groups of
customer resources, or any other construct, depending on the Amazon Web Services service.
Events that have impact beyond that of the affected entities, or where the extent of impact
is unknown, include at least one entity indicating this. At least one event ARN is
required.    This API operation uses pagination. Specify the nextToken parameter in the
next request to return more results.   This operation supports resource-level permissions.
You can use this operation to allow or deny access to specific Health events. For more
information, see Resource- and action-based conditions in the Health User Guide.

# Arguments
- `filter`: Values to narrow the results returned. At least one event ARN is required.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_affected_entities(
    filter; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeAffectedEntities",
        Dict{String,Any}("filter" => filter);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_affected_entities(
    filter, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeAffectedEntities",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("filter" => filter), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns a list of entities that have been affected by one or more events for one or more
accounts in your organization in Organizations, based on the filter criteria. Entities can
refer to individual customer resources, groups of customer resources, or any other
construct, depending on the Amazon Web Services service. At least one event Amazon Resource
Name (ARN) and account ID are required. Before you can call this operation, you must first
enable Health to work with Organizations. To do this, call the
EnableHealthServiceAccessForOrganization operation from your organization's management
account.    This API operation uses pagination. Specify the nextToken parameter in the next
request to return more results.   This operation doesn't support resource-level
permissions. You can't use this operation to allow or deny access to specific Health
events. For more information, see Resource- and action-based conditions in the Health User
Guide.

# Arguments
- `organization_entity_filters`: A JSON set of elements including the awsAccountId and the
  eventArn.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_affected_entities_for_organization(
    organizationEntityFilters; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeAffectedEntitiesForOrganization",
        Dict{String,Any}("organizationEntityFilters" => organizationEntityFilters);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_affected_entities_for_organization(
    organizationEntityFilters,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return health(
        "DescribeAffectedEntitiesForOrganization",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("organizationEntityFilters" => organizationEntityFilters),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns the number of entities that are affected by each of the specified events.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"eventArns"`: A list of event ARNs (unique identifiers). For example:
  \"arn:aws:health:us-east-1::event/EC2/EC2_INSTANCE_RETIREMENT_SCHEDULED/EC2_INSTANCE_RETIREM
  ENT_SCHEDULED_ABC123-CDE456\",
  \"arn:aws:health:us-west-1::event/EBS/AWS_EBS_LOST_VOLUME/AWS_EBS_LOST_VOLUME_CHI789_JKL101
  \"
"""
function describe_entity_aggregates(; aws_config::AbstractAWSConfig=global_aws_config())
    return health(
        "DescribeEntityAggregates"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_entity_aggregates(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEntityAggregates",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns the number of events of each event type (issue, scheduled change, and account
notification). If no filter is specified, the counts of all events in each category are
returned.  This API operation uses pagination. Specify the nextToken parameter in the next
request to return more results.

# Arguments
- `aggregate_field`: The only currently supported value is eventTypeCategory.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: Values to narrow the results returned.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_event_aggregates(
    aggregateField; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventAggregates",
        Dict{String,Any}("aggregateField" => aggregateField);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_event_aggregates(
    aggregateField,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return health(
        "DescribeEventAggregates",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("aggregateField" => aggregateField), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns detailed information about one or more specified events. Information includes
standard event data (Amazon Web Services Region, service, and so on, as returned by
DescribeEvents), a detailed event description, and possible additional metadata that
depends upon the nature of the event. Affected entities are not included. To retrieve the
entities, use the DescribeAffectedEntities operation. If a specified event can't be
retrieved, an error message is returned for that event.  This operation supports
resource-level permissions. You can use this operation to allow or deny access to specific
Health events. For more information, see Resource- and action-based conditions in the
Health User Guide.

# Arguments
- `event_arns`: A list of event ARNs (unique identifiers). For example:
  \"arn:aws:health:us-east-1::event/EC2/EC2_INSTANCE_RETIREMENT_SCHEDULED/EC2_INSTANCE_RETIREM
  ENT_SCHEDULED_ABC123-CDE456\",
  \"arn:aws:health:us-west-1::event/EBS/AWS_EBS_LOST_VOLUME/AWS_EBS_LOST_VOLUME_CHI789_JKL101
  \"

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
"""
function describe_event_details(
    eventArns; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventDetails",
        Dict{String,Any}("eventArns" => eventArns);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_event_details(
    eventArns,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return health(
        "DescribeEventDetails",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("eventArns" => eventArns), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns detailed information about one or more specified events for one or more Amazon Web
Services accounts in your organization. This information includes standard event data (such
as the Amazon Web Services Region and service), an event description, and (depending on the
event) possible metadata. This operation doesn't return affected entities, such as the
resources related to the event. To return affected entities, use the
DescribeAffectedEntitiesForOrganization operation.  Before you can call this operation, you
must first enable Health to work with Organizations. To do this, call the
EnableHealthServiceAccessForOrganization operation from your organization's management
account.  When you call the DescribeEventDetailsForOrganization operation, specify the
organizationEventDetailFilters object in the request. Depending on the Health event type,
note the following differences:   To return event details for a public event, you must
specify a null value for the awsAccountId parameter. If you specify an account ID for a
public event, Health returns an error message because public events aren't specific to an
account.   To return event details for an event that is specific to an account in your
organization, you must specify the awsAccountId parameter in the request. If you don't
specify an account ID, Health returns an error message because the event is specific to an
account in your organization.    For more information, see Event.  This operation doesn't
support resource-level permissions. You can't use this operation to allow or deny access to
specific Health events. For more information, see Resource- and action-based conditions in
the Health User Guide.

# Arguments
- `organization_event_detail_filters`: A set of JSON elements that includes the
  awsAccountId and the eventArn.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
"""
function describe_event_details_for_organization(
    organizationEventDetailFilters; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventDetailsForOrganization",
        Dict{String,Any}(
            "organizationEventDetailFilters" => organizationEventDetailFilters
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_event_details_for_organization(
    organizationEventDetailFilters,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return health(
        "DescribeEventDetailsForOrganization",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "organizationEventDetailFilters" => organizationEventDetailFilters
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Returns the event types that meet the specified filter criteria. You can use this API
operation to find information about the Health event, such as the category, Amazon Web
Services service, and event code. The metadata for each event appears in the EventType
object.  If you don't specify a filter criteria, the API operation returns all event types,
in no particular order.   This API operation uses pagination. Specify the nextToken
parameter in the next request to return more results.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: Values to narrow the results returned.
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_event_types(; aws_config::AbstractAWSConfig=global_aws_config())
    return health(
        "DescribeEventTypes"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_event_types(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventTypes", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

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

 Returns information about events that meet the specified filter criteria. Events are
returned in a summary form and do not include the detailed description, any additional
metadata that depends on the event type, or any affected resources. To retrieve that
information, use the DescribeEventDetails and DescribeAffectedEntities operations. If no
filter criteria are specified, all events are returned. Results are sorted by
lastModifiedTime, starting with the most recent event.    When you call the DescribeEvents
operation and specify an entity for the entityValues parameter, Health might return public
events that aren't specific to that resource. For example, if you call DescribeEvents and
specify an ID for an Amazon Elastic Compute Cloud (Amazon EC2) instance, Health might
return events that aren't specific to that resource or service. To get events that are
specific to a service, use the services parameter in the filter object. For more
information, see Event.   This API operation uses pagination. Specify the nextToken
parameter in the next request to return more results.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: Values to narrow the results returned.
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_events(; aws_config::AbstractAWSConfig=global_aws_config())
    return health("DescribeEvents"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function describe_events(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEvents", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

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

Returns information about events across your organization in Organizations. You can use
thefilters parameter to specify the events that you want to return. Events are returned in
a summary form and don't include the affected accounts, detailed description, any
additional metadata that depends on the event type, or any affected resources. To retrieve
that information, use the following operations:    DescribeAffectedAccountsForOrganization
   DescribeEventDetailsForOrganization     DescribeAffectedEntitiesForOrganization    If
you don't specify a filter, the DescribeEventsForOrganizations returns all events across
your organization. Results are sorted by lastModifiedTime, starting with the most recent
event.  For more information about the different types of Health events, see Event. Before
you can call this operation, you must first enable Health to work with Organizations. To do
this, call the EnableHealthServiceAccessForOrganization operation from your organization's
management account.  This API operation uses pagination. Specify the nextToken parameter in
the next request to return more results.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filter"`: Values to narrow the results returned.
- `"locale"`: The locale (language) to return information in. English (en) is the default
  and the only supported value at this time.
- `"maxResults"`: The maximum number of items to return in one batch, between 10 and 100,
  inclusive.
- `"nextToken"`: If the results of a search are large, only a portion of the results are
  returned, and a nextToken pagination token is returned in the response. To retrieve the
  next batch of results, reissue the search request and include the returned token. When all
  results have been returned, the response does not contain a pagination token value.
"""
function describe_events_for_organization(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventsForOrganization";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_events_for_organization(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeEventsForOrganization",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation provides status information on enabling or disabling Health to work with
your organization. To call this operation, you must sign in as an IAM user, assume an IAM
role, or sign in as the root user (not recommended) in the organization's management
account.

"""
function describe_health_service_status_for_organization(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeHealthServiceStatusForOrganization";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_health_service_status_for_organization(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DescribeHealthServiceStatusForOrganization",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Disables Health from working with Organizations. To call this operation, you must sign in
as an Identity and Access Management (IAM) user, assume an IAM role, or sign in as the root
user (not recommended) in the organization's management account. For more information, see
Aggregating Health events in the Health User Guide. This operation doesn't remove the
service-linked role from the management account in your organization. You must use the IAM
console, API, or Command Line Interface (CLI) to remove the service-linked role. For more
information, see Deleting a Service-Linked Role in the IAM User Guide.  You can also
disable the organizational feature by using the Organizations DisableAWSServiceAccess API
operation. After you call this operation, Health stops aggregating events for all other
Amazon Web Services accounts in your organization. If you call the Health API operations
for organizational view, Health returns an error. Health continues to aggregate health
events for your Amazon Web Services account.

"""
function disable_health_service_access_for_organization(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DisableHealthServiceAccessForOrganization";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function disable_health_service_access_for_organization(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "DisableHealthServiceAccessForOrganization",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Enables Health to work with Organizations. You can use the organizational view feature to
aggregate events from all Amazon Web Services accounts in your organization in a
centralized location.  This operation also creates a service-linked role for the management
account in the organization.   To call this operation, you must meet the following
requirements:   You must have a Business, Enterprise On-Ramp, or Enterprise Support plan
from Amazon Web Services Support to use the Health API. If you call the Health API from an
Amazon Web Services account that doesn't have a Business, Enterprise On-Ramp, or Enterprise
Support plan, you receive a SubscriptionRequiredException error.   You must have permission
to call this operation from the organization's management account. For example IAM
policies, see Health identity-based policy examples.    If you don't have the required
support plan, you can instead use the Health console to enable the organizational view
feature. For more information, see Aggregating Health events in the Health User Guide.

"""
function enable_health_service_access_for_organization(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "EnableHealthServiceAccessForOrganization";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function enable_health_service_access_for_organization(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return health(
        "EnableHealthServiceAccessForOrganization",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
abstract BilinearForm
abstract LinearForm

"""
Represents a weak form of a differential equation.

To provide flexibility, different operators can be used.
This type is used to store operators of the following form:

$$
\sum_{i=1}^N A_i(u,v) = \sum_{k=1}^M b_k(v)
$$

"""
type WeakForm{Dof<:DofMap}
    A::Vector{BilinearForm}
    Ab::Dict{Int,Vector{BilinearForm}}
    F::Vector{LinearForm}
    Fb::Dict{Int,Vector{LinearForm}}
    dlift::Dict{Int,DirichiletLift}
end

ext: jl


text:
@testset "Tree" begin #Tree, Ensamble, Classification module test
    x = [5.1 3.5; 4.9 3.0; 4.7 3.2; 4.6 3.1; 5.0 3.6; 5.4 3.9; 4.6 3.4; 5.0 3.4; 4.4 2.9; 4.9 3.1; 5.4 3.7; 4.8 3.4; 4.8 3.0; 4.3 3.0; 5.8 4.0; 5.7 4.4; 5.4 3.9; 5.1 3.5; 5.7 3.8; 5.1 3.8; 5.4 3.4; 5.1 3.7; 4.6 3.6; 5.1 3.3; 4.8 3.4; 5.0 3.0; 5.0 3.4; 5.2 3.5; 5.2 3.4; 4.7 3.2; 4.8 3.1; 5.4 3.4; 5.2 4.1; 5.5 4.2; 4.9 3.1; 5.0 3.2; 5.5 3.5; 4.9 3.6; 4.4 3.0; 5.1 3.4; 5.0 3.5; 4.5 2.3; 4.4 3.2; 5.0 3.5; 5.1 3.8; 4.8 3.0; 5.1 3.8; 4.6 3.2; 5.3 3.7; 5.0 3.3; 7.0 3.2; 6.4 3.2; 6.9 3.1; 5.5 2.3; 6.5 2.8; 5.7 2.8; 6.3 3.3; 4.9 2.4; 6.6 2.9; 5.2 2.7; 5.0 2.0; 5.9 3.0; 6.0 2.2; 6.1 2.9; 5.6 2.9; 6.7 3.1; 5.6 3.0; 5.8 2.7; 6.2 2.2; 5.6 2.5; 5.9 3.2; 6.1 2.8; 6.3 2.5; 6.1 2.8; 6.4 2.9; 6.6 3.0; 6.8 2.8; 6.7 3.0; 6.0 2.9; 5.7 2.6; 5.5 2.4; 5.5 2.4; 5.8 2.7; 6.0 2.7; 5.4 3.0; 6.0 3.4; 6.7 3.1; 6.3 2.3; 5.6 3.0; 5.5 2.5; 5.5 2.6; 6.1 3.0; 5.8 2.6; 5.0 2.3; 5.6 2.7; 5.7 3.0; 5.7 2.9; 6.2 2.9; 5.1 2.5; 5.7 2.8; 6.3 3.3; 5.8 2.7; 7.1 3.0; 6.3 2.9; 6.5 3.0; 7.6 3.0; 4.9 2.5; 7.3 2.9; 6.7 2.5; 7.2 3.6; 6.5 3.2; 6.4 2.7; 6.8 3.0; 5.7 2.5; 5.8 2.8; 6.4 3.2; 6.5 3.0; 7.7 3.8; 7.7 2.6; 6.0 2.2; 6.9 3.2; 5.6 2.8; 7.7 2.8; 6.3 2.7; 6.7 3.3; 7.2 3.2; 6.2 2.8; 6.1 3.0; 6.4 2.8; 7.2 3.0; 7.4 2.8; 7.9 3.8; 6.4 2.8; 6.3 2.8; 6.1 2.6; 7.7 3.0; 6.3 3.4; 6.4 3.1; 6.0 3.0; 6.9 3.1; 6.7 3.1; 6.9 3.1; 5.8 2.7; 6.8 3.2; 6.7 3.3; 6.7 3.0; 6.3 2.5; 6.5 3.0; 6.2 3.4; 5.9 3.0]
    t = [0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2]
    tree, classes = Tree.DT.fit(x, t)
    p = Tree.DT.predict(x, tree, classes)
    println(p)
    Tree.MV("./test.dot", tree)
    forest, classes, using_feature = Ensemble.RF.fit(x, t)
    p = Ensemble.RF.predict(x, forest, classes, using_feature)
    println(p)
    label = ["Apple", "Apple", "Pear", "Pear", "Lemon", "Apple", "Pear", "Lemon"]
    l = Classification.LE(label)
    println(l)
    t = Classification.CTOH(t)
    w = Classification.MS.fit(x, t, alpha = 0.1)
    println(w)
    p = Classification.MS.predict(x, w)
    println(p)
    cee = LossFunction.CEE(p, t, t_f = true)
    println(cee)
    w = Classification.OVR.fit(x, t; alpha = 0.1)
    println("w:", w)
    p = Classification.OVR.predict(x, w)
    println("predict:", p)
    cee = LossFunction.CEE(p, t)
    println("cee:", cee)
end


text:
# One action, zero observation, one timestep long, +1 reward every timestep: This isolates
# the value network. If my agent can't learn that the value of the only observation it ever
# sees it 1, there's a problem with the value loss calculation or the optimizer.
import AlphaZero.GI

struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  terminated :: Bool
end

const INITSTATE = (; terminated=false)

GI.init(::GameSpec, state=INITSTATE) = GameEnv(state.terminated)

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = false

const ACTIONS = collect(1:1)

GI.actions(::GameSpec) = ACTIONS

GI.actions_mask(g::GameEnv) = [true]

GI.current_state(g::GameEnv) = (; terminated=g.terminated)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.terminated

GI.white_reward(g::GameEnv) = 1.

function GI.play!(g::GameEnv, action)
  g.terminated = true
end

GI.vectorize_state(::GameSpec, state) = Float32[!state.terminated]


text:
# Univariate

const VectorOfUnivariate = Distributions.Product

function arraydist(dists::AbstractVector{<:UnivariateDistribution})
    return Product(dists)
end

function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractMatrix{<:Real})
    size(x, 1) == length(dist) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    # `eachcol` breaks Zygote, so we use `view` directly
    return map(i -> sum(map(logpdf, dist.v, view(x, :, i))), axes(x, 2))
end

struct MatrixOfUnivariate{
    S <: ValueSupport,
    Tdist <: UnivariateDistribution{S},
    Tdists <: AbstractMatrix{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::MatrixOfUnivariate) = size(dist.dists)
function arraydist(dists::AbstractMatrix{<:UnivariateDistribution})
    return MatrixOfUnivariate(dists)
end
function Distributions._logpdf(dist::MatrixOfUnivariate, x::AbstractMatrix{<:Real})
    # return sum(((d, xi),) -> logpdf(d, xi), zip(dist.dists, x))
    # Broadcasting here breaks Tracker for some reason
    return sum(map(logpdf, dist.dists, x))
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end

function Distributions.rand(rng::Random.AbstractRNG, dist::MatrixOfUnivariate)
    return rand.(Ref(rng), dist.dists)
end

# Multivariate

struct VectorOfMultivariate{
    S <: ValueSupport,
    Tdist <: MultivariateDistribution{S},
    Tdists <: AbstractVector{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::VectorOfMultivariate) = (length(dist.dists[1]), length(dist))
Base.length(dist::VectorOfMultivariate) = length(dist.dists)
function arraydist(dists::AbstractVector{<:MultivariateDistribution})
    return VectorOfMultivariate(dists)
end

function Distributions._logpdf(dist::VectorOfMultivariate, x::AbstractMatrix{<:Real})
    # `eachcol` breaks Zygote, so we use `view` directly
    return sum(i -> logpdf(dist.dists[i], view(x, :, i)), axes(x, 2))
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end

function Distributions.rand(rng::Random.AbstractRNG, dist::VectorOfMultivariate)
    init = reshape(rand(rng, dist.dists[1]), :, 1)
    return mapreduce(i -> rand(rng, dist.dists[i]), hcat, 2:length(dist); init = init)
end


text:
module QuantumNumbers

using Base.Iterators: Reverse, flatten, product, reverse
using Printf: @printf, @sprintf
using DataStructures: OrderedDict
using LinearAlgebra: norm
using ...Prerequisites: Float
using ...Prerequisites.Combinatorics: Combinations
using ...Prerequisites.NamedVectors: HomoNamedVector
using ...Prerequisites.VectorSpaces: EnumerativeVectorSpace

import Random
import ...Prerequisites.Traits: contentnames, getcontent
import ...Interfaces: ⊕, ⊗, dimension, expand, permute, decompose

export AbelianNumber, AbelianNumbers, periods, regularize, regularize!, @abeliannumber
export SQN, PQN, SPQN, SQNS, PQNS, SPQNS, Momentum, Momentum₁, Momentum₂, Momentum₃

"""
    positives(inputs::NTuple{N, Any}) where N -> NTuple{N, Int}

Return a tuple of all positive signs.
"""
@inline positives(inputs::NTuple{N, Any}) where N = ntuple(i->1, Val(N))

"Abstract type for all concrete quantum numbers for a single basis."
abstract type AbelianNumber <: HomoNamedVector{Float} end

"""
    periods

The generic function that used to define the period of each component of a concrete `AbelianNumber`.
"""
function periods end

"""
    dimension(::Type{<:AbelianNumber}) -> Int
    dimension(::AbelianNumber) -> Int

The dimension of the Hilbert space an `AbelianNumber` represents. Apparently, this is always 1.
"""
@inline dimension(::Type{<:AbelianNumber}) = 1
@inline dimension(::AbelianNumber) = 1

"""
    +(qn::AbelianNumber) -> typeof(qn)
    +(qn::QN, qns::QN...) where {QN<:AbelianNumber} -> QN

Overloaded `+` operator for `AbelianNumber`.

!!! note
    To ensure type stability, two `AbelianNumber` can be added together if and only if they are of the same type.
"""
@inline Base.:+(qn::AbelianNumber) = qn
@inline Base.:+(qn::QN, qns::QN...) where {QN<:AbelianNumber} = map(+, qn, qns...)

"""
    -(qn::AbelianNumber) -> typeof(qn)
    -(qn1::QN, qn2::QN) where {QN<:AbelianNumber} -> QN

Overloaded `-` operator for `AbelianNumber` and `AbelianNumbers`.

!!! note
    To ensure type stability, an `AbelianNumber` can be subtracted by another `AbelianNumber` if and only if they are of the same type.
"""
@inline Base.:-(qn::AbelianNumber) = map(-, qn)
@inline Base.:-(qn1::QN, qn2::QN) where {QN<:AbelianNumber} = map(-, qn1, qn2)

"""
    *(factor::Integer, qn::AbelianNumber) -> typeof(qn)
    *(qn::AbelianNumber, factor::Integer) -> typeof(qn)

Overloaded `*` operator for the multiplication between an integer and an `AbelianNumber`.
"""
@inline Base.:*(factor::Integer, qn::AbelianNumber) = qn * factor
@inline @generated function Base.:*(qn::AbelianNumber, factor::Integer)
    exprs = Expr[:(getfield(qn, $i)*factor) for i = 1:(qn|>fieldnames|>length)]
    return :(typeof(qn)($(exprs...)))
end

"""
    ^(qn::AbelianNumber, factor::Integer) -> typeof(qn)

Overloaded `^` operator for `AbelianNumber`.

This operation translates into the direct product of `factor` copies of `qn`.
"""
@inline Base.:^(qn::AbelianNumber, factor::Integer) = kron(NTuple{factor, typeof(qn)}(qn for i = 1:factor)...)

"""
    ⊗(qns::AbelianNumber...) -> eltype(qns)

Get the direct product of some `AbelianNumber`s.
"""
⊗(qns::AbelianNumber...) = kron(qns...)

"""
    kron(qns::Vararg{AbelianNumber}; signs=positives(qns))-> eltype(qns)

Get the direct product of some `AbelianNumber`s.
!!! note
    Physically, the direct product of a couple of `AbelianNumber`s are defined through the direct product of the bases of the Hilbert spaces they represent. Apparently, the result is still an `AbelianNumber` whose dimension is 1. At the same time, each component of the result is obtained by a summation of the corresponding components of the inputs with the correct signs. This is a direct consequence of the Abelian nature of our quantum numbers.
"""
Base.kron(qns::Vararg{AbelianNumber}; signs=positives(qns)) = sum(sign==1 ? qn : -qn for (sign, qn) in zip(signs, qns))

"""
    regularize!(::Type{QN}, array::AbstractVector{<:Real}) where {QN<:AbelianNumber} -> typeof(array)
    regularize!(::Type{QN}, array::AbstractMatrix{<:Real}) where {QN<:AbelianNumber} -> typeof(array)

Regularize the elements of an array in place so that it can represent quantum numbers.
"""
function regularize!(::Type{QN}, array::AbstractVector{<:Real}) where {QN<:AbelianNumber}
    @assert array|>length == QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i, period) in enumerate(QN|>periods)
        period==Inf || @inbounds begin
            remainder = array[i] % period
            array[i] = remainder<0 ? remainder+period : remainder
        end
    end
    return array
end
function regularize!(::Type{QN}, array::AbstractMatrix{<:Real}) where {QN<:AbelianNumber}
    @assert size(array, 1) == QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i, period) in enumerate(QN|>periods)
        period==Inf || @inbounds for j = 1:size(array, 2)
            remainder = array[i, j] % period
            array[i, j] = remainder<0 ? remainder+period : remainder
        end
    end
    return array
end

"""
    regularize(::Type{QN}, array::Union{AbstractVector{<:Real}, AbstractMatrix{<:Real}}) where {QN<:AbelianNumber} -> typeof(array)

Regularize the elements of an array and return a copy that can represent quantum numbers.
"""
regularize(::Type{QN}, array::Union{AbstractVector{<:Real}, AbstractMatrix{<:Real}}) where {QN<:AbelianNumber} = regularize!(QN, deepcopy(array))

"""
    @abeliannumber typename fields periods

Construct a concrete `AbelianNumber` with the type name being `typename`, fieldnames specified by `fields` and periods specified by `periods`.
"""
macro abeliannumber(typename, fields, periods)
    typename = Symbol(typename)
    fields = tuple(eval(fields)...)
    periods = tuple(eval(periods)...)
    arguments = ntuple(i->Symbol(:v, i), length(fields))
    @assert length(fields) == length(periods) "@abeliannumber error: number of fields($(length(fields))) and periods($(length(periods))) not equal."
    @assert all(isa(name, Symbol) for name in fields) "@abeliannumber error: all field names should be Symbol."
    @assert all(periods .> 0) "@abeliannumber error: all field periods should be positive."
    if all(periods .== Inf)
        title = Expr(:call, :($(esc(typename))), arguments...)
        body = Expr(:call, :new, arguments...)
    else
        title = Expr(:call, :($(esc(typename))), arguments..., Expr(:kw, :(regularize::Bool), :true))
        body = Expr(:call, :new, (p==Inf ? arg : :(regularize ? (r=$arg%$p; r<0 ? r+$p : r) : $arg) for (p, arg) in zip(periods, arguments))...)
    end
    newtype = Expr(:struct, false, :($(esc(typename))<:AbelianNumber), Expr(:block, (:($field::Float) for field in fields)..., Expr(:(=), title, body)))
    functions = Expr(:block, :(periods(::Type{<:$(esc(typename))}) = $periods))
    return Expr(:block, :(global periods), :(Base.@__doc__($newtype)), functions)
end

"""
    AbelianNumbers{QN<:AbelianNumber} <: EnumerativeVectorSpace{QN}

The whole quantum numbers of the total bases of a Hilbert space.
"""
struct AbelianNumbers{QN<:AbelianNumber} <: EnumerativeVectorSpace{QN}
    form::Char
    contents::Vector{QN}
    indptr::Vector{Int}
    function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int})
        @assert uppercase(form)∈('G', 'U', 'C') "AbelianNumbers error: 'form'($form) is not 'G', 'U' or 'C'."
        @assert length(contents)+1==length(indptr) "AbelianNumbers error: mismatch shapes of contents and indptr ($(length(contents))+1!=$length(indptr))."
        return new{contents|>eltype}(form|>uppercase, contents, indptr)
    end
end
@inline contentnames(::Type{<:AbelianNumbers}) = (:form, :table, :indptr)
@inline getcontent(qns::AbelianNumbers, ::Val{:table}) = qns.contents
@inline dimension(qns::AbelianNumbers) = @inbounds(qns.indptr[end])
@inline Base.size(qns::AbelianNumbers) = (length(qns.contents),)
@inline Base.issorted(qns::AbelianNumbers) = qns.form=='C'
Base.string(qns::AbelianNumbers) = @sprintf "QNS(%s, %s)" qns|>length qns|>dimension
Base.show(io::IO, qns::AbelianNumbers) = @printf io "QNS(%s)" join((@sprintf("%s=>%s", qn, slice) for (qn, slice) in pairs(qns, :indptr)), ", ")

"""
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, choice::Symbol)
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, counts::Vector{Int}, ::Val{:counts})
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, ::Val{:indptr})

Construct an `AbelianNumbers` from a vector of concrete quantum numbers and an vector containing their counts or indptr.
"""
@inline function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, choice::Symbol)
    return AbelianNumbers(form, contents, indptr, choice|>Val)
end
@inline function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, counts::Vector{Int}, ::Val{:counts})
    return AbelianNumbers(form|>uppercase, contents, [0, cumsum(counts)...])
end
@inline function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, ::Val{:indptr})
    return AbelianNumbers(form|>uppercase, contents, indptr)
end

"""
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber})

Construct an `AbelianNumbers` with a set of quantum numbers whose counts are all one.
"""
@inline function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber})
    return AbelianNumbers(form, contents, collect(0:length(contents)), Val(:indptr))
end

"""
    AbelianNumbers(qn::AbelianNumber, count::Int=1)

Construct an `AbelianNumbers` with one unique quantum number which occurs `count` times.
"""
@inline AbelianNumbers(qn::AbelianNumber, count::Int=1) = AbelianNumbers('C', [qn], [0, count], :indptr)

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:Integer})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their counts.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:Integer})
    contents = Vector{od|>keytype}(undef, length(od))
    indptr = zeros(Int, length(od)+1)
    for (i, (qn, count)) in enumerate(od)
        @inbounds contents[i] = qn
        @inbounds indptr[i+1] = indptr[i] + count
    end
    return AbelianNumbers('U', contents, indptr, :indptr)
end

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:UnitRange{<:Integer}})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their slices.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:UnitRange{<:Integer}})
    contents = Vector{od|>keytype}(undef, length(od))
    indptr = zeros(Int, length(od)+1)
    for (i, (qn, slice)) in enumerate(od)
        @inbounds contents[i] = qn
        @inbounds @assert indptr[i]+1==slice.start "AbelianNumbers error: slice not consistent."
        @inbounds indptr[i+1] = slice.stop
    end
    return AbelianNumbers('U', contents, indptr, :indptr)
end

"""
    getindex(qns::AbelianNumbers, slice::UnitRange{Int}) -> AbelianNumbers
    getindex(qns::AbelianNumbers, indices::Vector{Int}) -> AbelianNumbers

Overloaded `[]` operator.
!!! note
    For an `AbelianNumbers`, all the `getindex` functions act on its `contents`, i.e. its compressed data, but not on its expansion, i.e. the uncompressed data. This definition is consistent with the size of an `AbelianNumbers`.
"""
function Base.getindex(qns::AbelianNumbers, slice::UnitRange{Int})
    contents = qns.contents[slice]
    indptr = qns.indptr[slice.start:slice.stop+1]
    indptr .= indptr .- qns.indptr[slice.start]
    return AbelianNumbers(qns.form, contents, indptr, :indptr)
end
function Base.getindex(qns::AbelianNumbers, indices::Vector{Int})
    contents = Vector{qns|>eltype}(undef, length(indices))
    indptr = zeros(Int, length(indices)+1)
    for (i, index) in enumerate(indices)
        contents[i] = qns.contents[index]
        indptr[i+1] = indptr[i] + qns.indptr[index+1] - qns.indptr[index]
    end
    return AbelianNumbers('G', contents, indptr, :indptr)
end

"""
    keys(qns::AbelianNumbers) -> Vector{qns|>eltype}

Iterate over the concrete `AbelianNumber`s contained in an `AbelianNumbers`.
"""
@inline Base.keys(qns::AbelianNumbers) = qns.contents

"""
    values(qns::AbelianNumbers, choice::Symbol)
    values(qns::AbelianNumbers, ::Val{:indptr})
    values(qns::AbelianNumbers, ::Val{:counts})

Iterate over the slices/counts of the `AbelianNumbers`.
"""
@inline Base.values(qns::AbelianNumbers, choice::Symbol) = values(qns, choice|>Val)
@inline @views Base.values(qns::AbelianNumbers, ::Val{:indptr}) = ((start+1):stop for (start, stop) in zip(qns.indptr[1:end-1], qns.indptr[2:end]))
@inline @views Base.values(qns::AbelianNumbers, ::Val{:counts}) = (stop-start for (start, stop) in zip(qns.indptr[1:end-1], qns.indptr[2:end]))

"""
    pairs(qns::AbelianNumbers, choice::Symbol)
    pairs(qns::AbelianNumbers, ::Val{:indptr})
    pairs(qns::AbelianNumbers, ::Val{:counts})

Iterate over the `AbelianNumber=>slice` or `AbelianNumber=>count` pairs.
"""
@inline Base.pairs(qns::AbelianNumbers, choice::Symbol) = pairs(qns, choice|>Val)
@inline Base.pairs(qns::AbelianNumbers, choice::Union{Val{:indptr}, Val{:counts}}) = Base.Generator(=>, keys(qns), values(qns, choice))

"""
    OrderedDict(qns::AbelianNumbers, choice::Symbol) -> OrderedDict
    OrderedDict(qns::AbelianNumbers, ::Val{:indptr}) -> OrderedDict{qns|>eltype, UnitRange{Int}}
    OrderedDict(qns::AbelianNumbers, ::Val{:counts}) -> OrderedDict{qns|>eltype, Int}

Convert an `AbelianNumbers` to an ordered dict.
"""
@inline OrderedDict(qns::AbelianNumbers, choice::Symbol) = OrderedDict(qns, choice|>Val)
function OrderedDict(qns::AbelianNumbers, ::Val{:indptr})
    @assert qns.form≠'G' "OrderedDict error: input `AbelianNumbers` cannot be `G` formed."
    result = OrderedDict{qns|>eltype, UnitRange{Int}}()
    @inbounds for i = 1:length(qns)
        result[qns.contents[i]] = qns.indptr[i]+1:qns.indptr[i+1]
    end
    return result
end
function OrderedDict(qns::AbelianNumbers, ::Val{:counts})
    @assert qns.form≠'G' "OrderedDict error: input `AbelianNumbers` cannot be `G` formed."
    result = OrderedDict{qns|>eltype, Int}()
    @inbounds for i = 1:length(qns)
        result[qns.contents[i]] = qns.indptr[i+1] - qns.indptr[i]
    end
    return result
end

"""
    sort(qns::AbelianNumbers) -> Tuple{AbelianNumbers, Vector{Int}}

Sort the quantum numbers of an `AbelianNumbers`, return the sorted `AbelianNumbers` and the permutation array that sorts the expansion of the original `AbelianNumbers`.
"""
function Base.sort(qns::AbelianNumbers)
    ctpts = sortperm(qns.contents, alg=Base.Sort.QuickSort)
    masks = Vector{Bool}(undef, length(qns.contents))
    masks[1] = true
    nonduplicate = 1
    @inbounds for i = 2:length(ctpts)
        masks[i] = qns.contents[ctpts[i]]≠qns.contents[ctpts[i-1]]
        masks[i] && (nonduplicate += 1)
    end
    contents = Vector{qns|>eltype}(undef, nonduplicate)
    indptr = zeros(Int, nonduplicate+1)
    permutation = Vector{Int}(undef, dimension(qns))
    qncount, ptcount = 0, 0
    @inbounds for (mask, index) in zip(masks, ctpts)
        mask && (qncount += 1; contents[qncount] = qns.contents[index]; indptr[qncount+1] = indptr[qncount])
        indptr[qncount+1] += qns.indptr[index+1] - qns.indptr[index]
        for (i, p) in enumerate(qns.indptr[index]+1:qns.indptr[index+1])
            permutation[ptcount+i] = p
        end
        ptcount += qns.indptr[index+1] - qns.indptr[index]
    end
    return AbelianNumbers('C', contents, indptr, :indptr), permutation
end

"""
    findall(target::Union{QN, Tuple{Vararg{QN}}}, qns::AbelianNumbers{QN}, choice::Symbol) where {QN<:AbelianNumber} -> Vector{Int}
    findall(target::Union{QN, Tuple{Vararg{QN}}}, qns::AbelianNumbers{QN}, ::Val{:compression}) where {QN<:AbelianNumber} -> Vector{Int})
    findall(target::Union{QN, Tuple{Vararg{QN}}}, qns::AbelianNumbers{QN}, ::Val{:expansion}) where {QN<:AbelianNumber} -> Vector{Int}

Find all the indices of the target quantum numbers in the contents or the expansion of an `AbelianNumbers`.
"""
@inline Base.findall(target::Union{QN, Tuple{Vararg{QN}}}, qns::AbelianNumbers{QN}, choice::Symbol) where {QN<:AbelianNumber} = findall(target, qns, choice|>Val)
@inline Base.findall(target::QN, qns::AbelianNumbers{QN}, ::Val{:compression}) where {QN<:AbelianNumber} = findall((target,), qns, :compression|>Val)
@inline Base.findall(target::QN, qns::AbelianNumbers{QN}, ::Val{:expansion}) where {QN<:AbelianNumber} = findall((target,), qns, :expansion|>Val)
function Base.findall(targets::Tuple{Vararg{QN}}, qns::AbelianNumbers{QN}, ::Val{:compression}) where {QN<:AbelianNumber}
    result = Int[]
    if issorted(qns)
        for qn in targets
            range = searchsorted(qns.contents, qn)
            range.start<=range.stop && push!(result, range.start)
        end
    else
        for (i, qn) in enumerate(qns)
            qn∈targets && push!(result, i)
        end
    end
    return result
end
function Base.findall(targets::Tuple{Vararg{QN}}, qns::AbelianNumbers{QN}, ::Val{:expansion}) where {QN<:AbelianNumber}
    result = Int[]
    @inbounds for index in findall(targets, qns, :compression)
        append!(result, qns.indptr[index]+1:qns.indptr[index+1])
    end
    return result
end

"""
    filter(target::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} -> AbelianNumbers{QN}
    filter(targets::Tuple{Vararg{QN}}, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} -> AbelianNumbers{QN}

Find a subset of an `AbelianNumbers` by picking out the target quantum numbers.
"""
@inline Base.filter(target::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = qns[findall(target, qns, :compression)]
@inline Base.filter(targets::Tuple{Vararg{QN}}, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = qns[findall(targets, qns, :compression)]

"""
    permute(qns::AbelianNumbers, permutation::Vector{Int}, choice::Symbol) -> AbelianNumbers
    permute(qns::AbelianNumbers, permutation::Vector{Int}, ::Val{:compression}) -> AbelianNumbers
    permute(qns::AbelianNumbers, permutation::Vector{Int}, ::Val{:expansion}) -> AbelianNumbers

Reorder the quantum numbers contained in an `AbelianNumbers` with a permutation and return the new one.

For `:compression` case, the permutation is for the compressed contents of the original `AbelianNumbers` while for `:expansion` case, the permutation is for the expanded contents of the original `AbelianNumbers`.
"""
@inline permute(qns::AbelianNumbers, permutation::Vector{Int}, choice::Symbol) = permute(qns, permutation, choice|>Val)
@inline permute(qns::AbelianNumbers, permutation::Vector{Int}, ::Val{:compression}) = qns[permutation]
function permute(qns::AbelianNumbers, permutation::Vector{Int}, ::Val{:expansion})
    contents = Vector{qns|>eltype}(undef, length(permutation))
    indptr = zeros(Int, length(permutation)+1)
    expansion = expand(qns, :indices)
    @inbounds for (i, p) in enumerate(permutation)
        contents[i] = qns.contents[expansion[p]]
        indptr[i+1] = indptr[i] + 1
    end
    return AbelianNumbers('G', contents, indptr, :indptr)
end

"""
    expand(qns::AbelianNumbers, choice::Symbol)
    expand(qns::AbelianNumbers, ::Val{:contents}) -> Vector{eltype(qns)}
    expand(qns::AbelianNumbers, ::Val{:indices}) -> Vector{Int}

Expand the contents or indices of an `AbelianNumbers` to the uncompressed form.
"""
@inline expand(qns::AbelianNumbers, choice::Symbol) = expand(qns, choice|>Val)
function expand(qns::AbelianNumbers, ::Val{:contents})
    result = Vector{qns|>eltype}(undef, dimension(qns))
    @inbounds for i = 1:length(qns)
        for j = qns.indptr[i]+1:qns.indptr[i+1]
            result[j] = qns.contents[i]
        end
    end
    return result
end
function expand(qns::AbelianNumbers, ::Val{:indices})
    result = Vector{Int}(undef, dimension(qns))
    @inbounds for i = 1:length(qns)
        for j = qns.indptr[i]+1:qns.indptr[i+1]
            result[j] = i
        end
    end
    return result
end

"""
    +(qns::AbelianNumbers) -> AbelianNumbers
    +(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} -> AbelianNumbers{QN}
    +(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber} -> AbelianNumbers{QN}

Overloaded `+` operator for `AbelianNumber`/`AbelianNumbers`.
!!! note
    1. The addition between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `+` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` for direct sum and `⊗` for direct product.
    3. To ensure type stability, an `AbelianNumber` and an `AbelianNumbers` can be added together if and only if the former's type is the same with the latter's eltype.
"""
@inline Base.:+(qns::AbelianNumbers) = qns
@inline Base.:+(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = qns + qn
@inline function Base.:+(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber}
    return AbelianNumbers(qns.form=='G' ? 'G' : 'U', [iqn+qn for iqn in qns], qns.indptr, :indptr)
end

"""
    -(qns::AbelianNumbers) -> AbelianNumbers
    -(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} -> AbelianNumbers{QN}
    -(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber} -> AbelianNumbers{QN}

Overloaded `-` operator for `AbelianNumber`/`AbelianNumbers`.
!!! note
    1. The subtraction between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `-` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` with signs for direct sum and `⊗` with signs for direct product.
    3. To ensure type stability, an `AbelianNumber` can be subtracted by an `AbelianNumbers` or vice versa if and only if the former's type is the same with the latter's eltype.
"""
@inline function Base.:-(qns::AbelianNumbers)
    return AbelianNumbers(qns.form=='G' ? 'G' : 'U', -qns.contents, qns.indptr, :indptr)
end
@inline function Base.:-(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber}
    return AbelianNumbers(qns.form=='G' ? 'G' : 'U', [qn-iqn for iqn in qns], qns.indptr, :indptr)
end
@inline function Base.:-(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber}
    return AbelianNumbers(qns.form=='G' ? 'G' : 'U', [iqn-qn for iqn in qns], qns.indptr, :indptr)
end

"""
    *(factor::Integer, qns::AbelianNumbers) -> AbelianNumbers
    *(qns::AbelianNumbers, factor::Integer) -> AbelianNumbers

Overloaded `*` operator for the multiplication between an integer and an `AbelianNumbers`.
"""
@inline Base.:*(factor::Integer, qns::AbelianNumbers) = qns * factor
@inline function Base.:*(qns::AbelianNumbers, factor::Integer)
    return AbelianNumbers(qns.form=='G' ? 'G' : 'U', [qn*factor for qn in qns.contents], qns.indptr, :indptr)
end

"""
    ^(qns::AbelianNumbers, factor::Integer) -> AbelianNumbers

Overloaded `^` operator for `AbelianNumbers`.

This operation translates into the direct product of `factor` copies of `qns`.
"""
@inline Base.:^(qns::AbelianNumbers, factor::Integer) = kron(NTuple{factor, qns|>typeof}(qns for i = 1:factor)...)

"""
    ⊕(qns::AbelianNumber...) -> AbelianNumbers{qns|>eltype}
    ⊕(qnses::AbelianNumbers...) -> qnses|>eltype

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
"""
@inline ⊕(qns::AbelianNumber...) = union(qns...)
@inline ⊕(qnses::AbelianNumbers...) = union(qnses...)

"""
    ⊗(qnses::AbelianNumbers...) -> eltype(qnses)

Get the direct product of some `AbelianNumbers`es.
"""
⊗(qnses::AbelianNumbers...) = kron(qnses...)

"""
    union(qns::AbelianNumber...; signs=positives(qns)) -> AbelianNumbers
    union(qnses::AbelianNumbers{QN}...; signs=positives(qnses)) where {QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
!!! note
    1. Physically, the direct sum of a couple of `AbelianNumber`s or `AbelianNumbers`es is defined by the direct sum of the bases of the Hilbert spaces they represent. Therefore, the dimension of the result equals the summation of those of the inputs. As a consequence, even for `AbelianNumber`s, the result will be an `AbelianNumbers` because the dimension of the result is greater than 1.
    2. Signs of `AbelianNumber`s or `AbelianNumbers`es can be provided when getting their direct sums.
"""
@inline function Base.union(qns::AbelianNumber...; signs=positives(qns))
    contents = [sign==1 ? qn : -qn for (sign, qn) in zip(signs, qns)]
    indptr = collect(0:length(qns))
    return AbelianNumbers('G', contents, indptr, :indptr)
end
function Base.union(qnses::AbelianNumbers{QN}...; signs=positives(qnses)) where {QN<:AbelianNumber}
    lengths = NTuple{fieldcount(typeof(qnses)), Int}(length(qns) for qns in qnses)
    contents = Vector{QN}(undef, sum(lengths))
    indptr = zeros(Int, sum(lengths)+1)
    count = 0
    @inbounds for (sign, qns, length) in zip(signs, qnses, lengths)
        for i = 1:length
            contents[count+i] = (sign == 1) ? qns.contents[i] : -qns.contents[i]
            indptr[count+i+1] = indptr[count+i] + qns.indptr[i+1] - qns.indptr[i]
        end
        count += length
    end
    return AbelianNumbers('G', contents, indptr, :indptr)
end

"""
    kron(qnses::AbelianNumbers{QN}...; signs=positives(qnses)) where {QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct product of some `AbelianNumbers`es.
!!! note
    Physically, the direct product of a couple of `AbelianNumbers`es are defined by the direct product of the bases of the Hilbert spaces they represent. Therefore, the dimension of the result equals the product of those of the inputs. Meanwhile, each quantum number in the contents of the result is obtained by a summation of the corresponding quantum numbers of the inputs with the correct signs. This is a direct consequence of the Abelian nature of our quantum numbers.
"""
function Base.kron(qnses::AbelianNumbers{QN}...; signs=positives(qnses)) where {QN<:AbelianNumber}
    N = fieldcount(typeof(qnses))
    lengths = NTuple{N, Int}(i<N ? dimension(qns) : length(qns) for (i, qns) in enumerate(qnses))
    contents = Vector{QN}(undef, prod(lengths))
    indptr = zeros(Int, prod(lengths)+1)
    cache = Vector{QN}(undef, N)
    @inbounds expansions = NTuple{N-1, Vector{Int}}(expand(qnses[i], :indices) for i = 1:N-1)
    @inbounds for (i, indices) in enumerate(product(NTuple{N, UnitRange{Int64}}(1:length for length in reverse(lengths))...))
        for (j, (sign, qns, index)) in enumerate(zip(signs, qnses, reverse(indices)))
            pos = j<N ? expansions[j][index] : index
            cache[j] = sign==1 ? qns.contents[pos] : -qns.contents[pos]
        end
        indptr[i+1] = indptr[i] + qnses[end].indptr[indices[1]+1] - qnses[end].indptr[indices[1]]
        contents[i] = sum(cache)
    end
    return AbelianNumbers('G', contents, indptr, :indptr)
end

"""
    prod(qnses::AbelianNumbers{QN}...;
        signs=positives(qnses)
        ) where {QN<:AbelianNumber} -> AbelianNumbers{QN}, Dict{QN, Dict{NTuple{length(qnses), QN}, UnitRange{Int}}}

Unitary Kronecker product of several `AbelianNumbers`es. The product result as well as the records of the product will be returned.
!!! note
    1. All input `AbelianNumbers` must be 'U' formed or 'C' formed.
    2. Since duplicate quantum number are not allowed in 'U' formed and 'C' formed `AbelianNumbers`es, in general, there exists a merge process of duplicate quantum numbers in the result. Therefore, records are needed to keep track of this process, which will be returned along with the product result. The records are stored in a `Dict{QN, Dict{NTuple{NTuple{length(qnses), QN}, UnitRange{Int}}}` typed dict, in which, for each nonduplicate quantum number `qn` in the result, there exist a record `Dict((qn₁, qn₂, ...)=>start:stop, ...)` telling what quantum numbers `(qn₁, qn₂, ...)` a merged duplicate `qn` comes from and what slice `start:stop` this merged duplicate corresponds in the result.
"""
function Base.prod(qnses::AbelianNumbers{QN}...; signs=positives(qnses)) where {QN<:AbelianNumber}
    @assert all((qns.form == 'U') || (qns.form == 'C') for qns in qnses) "prod error: all input qnses should be 'U' formed or 'C' formed."
    N = fieldcount(typeof(qnses))
    lengths = NTuple{N, Int}(length(qns) for qns in qnses)
    cache = Vector{QN}(undef, N)
    container = OrderedDict{QN, Int}()
    records = Dict{QN, Dict{NTuple{N, QN}, UnitRange{Int}}}()
    for indices in product(NTuple{N, UnitRange{Int64}}(1:length for length in reverse(lengths))...)
        qn, count = QN|>zero, 1
        @inbounds for (j, (sign, qns, index)) in enumerate(zip(signs, qnses, reverse(indices)))
            cache[j] = qns.contents[index]
            qn = sign==1 ? (qn+cache[j]) : (qn-cache[j])
            count = count * (qns.indptr[index+1]-qns.indptr[index])
        end
        container[qn] = get(container, qn, 0) + count
        !haskey(records, qn) && (records[qn] = Dict{NTuple{N, QN}, UnitRange{Int}}())
        records[qn][NTuple{N, QN}(cache)] = (container[qn]-count+1):container[qn]
    end
    contents = QN[qn for qn in keys(container)]
    sort!(contents)
    indptr = zeros(Int, length(container)+1)
    @inbounds for (i, qn) in enumerate(contents)
        indptr[i+1] = indptr[i] + container[qn]
    end
    return AbelianNumbers('C', contents, indptr, :indptr), records
end

"""
    decompose(target::QN, qnses::AbelianNumbers{QN}...;
        signs=positives(qnses),
        method=:montecarlo,
        nmax=20
        ) where {QN<:AbelianNumber} -> Vector{NTuple{length(qnses), Int}}

Find a couple of decompositions of `target` with respect to `qnses`.
!!! note
    A tuple of integers `(i₁, i₂, ...)` is called a decomposition of a given `target` with respect to the given `qnses` if and only if they satisfy the "decomposition rule":
    ```math
    \\sum_\\text{j} \\text{signs}[\\text{j}]\\times\\text{qnses}[\\text{j}][\\text{i}_{\\text{j}}]==\\text{target}
    ```
    This equation is in fact a set of restricted [linear Diophantine equations](https://en.wikipedia.org/wiki/Diophantine_equation#Linear_Diophantine_equations). Indeed, our quantum numbers are always discrete Abelian ones and all instances of a concrete `AbelianNumber` forms a [module](https://en.wikipedia.org/wiki/Module_(mathematics)) over the [ring](https://en.wikipedia.org/wiki/Ring_(mathematics)) of integers. Therefore, each quantum number can be represented as a integral multiple of the unit element of the Abelian module, which results in the final reduction of the above equation to a set of linear Diophantine equations. Then finding a decomposition is equivalent to find a solution of the reduced linear Diophantine equations, with the restriction that the quantum numbers constructed from the solution should be in the corresponding `qnses`. Here we provide two methods to find such decompositions, one is by brute force, and the other is by Monte Carlo simulations.
"""
@inline function decompose(target::QN, qnses::AbelianNumbers{QN}...; signs=positives(qnses), method=:montecarlo, nmax=20) where {QN<:AbelianNumber}
    _decompose(Val(method), target, qnses...; signs=signs, nmax=nmax)
end
function _decompose(::Val{:bruteforce}, target::QN, qnses::AbelianNumbers{QN}...; signs, nmax) where {QN<:AbelianNumber}
    N = fieldcount(typeof(qnses))
    result = Set{NTuple{N, Int}}()
    cache = Vector{Int}(undef, N)
    dimensions = NTuple{N, Int}(dimension(qns) for qns in reverse(qnses))
    indices = findall(target, kron(qnses...; signs=signs), :expansion)
    nmax<length(indices) && (Random.shuffle!(Random.MersenneTwister(), indices); indices = @views indices[1:nmax])
    for index in indices
        @inbounds for (i, dimension) in enumerate(dimensions)
            cache[end+1-i] = (index-1)%dimension + 1
            index = (index-1)÷dimension + 1
        end
        push!(result, NTuple{N, Int}(cache))
    end
    return collect(NTuple{N, Int}, result)
end
function _decompose(::Val{:montecarlo}, target::QN, qnses::AbelianNumbers{QN}...; signs, nmax) where {QN<:AbelianNumber}
    Random.seed!()
    N = fieldcount(typeof(qnses))
    result = Set{NTuple{N, Int}}()
    expansions = Vector{QN}[expand(qns, :contents) for qns in qnses]
    dimensions = NTuple{N, Int}(dimension(qns) for qns in qnses)
    diff = indices->norm(sum(sign==1 ? expansion[index] : -expansion[index] for (sign, expansion, index) in zip(signs, expansions, indices))-target)
    count = 1
    while true
        oldindices, newindices = Int[rand(1:dimension) for dimension in dimensions], ones(Int, N)
        olddiff, newdiff = diff(oldindices), diff(newindices)
        @inbounds while newdiff > 0
            pos = rand(1:N)
            index = rand(1:dimensions[pos]-1)
            newindices[pos] = index<newindices[pos] ? index : index+1
            newdiff = diff(newindices)
            if newdiff<=olddiff || exp(olddiff-newdiff)>rand()
                oldindices[:] = newindices
                olddiff = newdiff
            end
            newindices[:] = oldindices
        end
        count = count + 1
        push!(result, NTuple{N, Int}(newindices))
        (length(result)>=nmax || count>nmax*5) && break
    end
    return collect(NTuple{N, Int}, result)
end

"""
    SQN(Sz::Real)

The concrete `AbelianNumber` of a quantum system with spin z-component `Sz` conserved.
"""
@abeliannumber "SQN" (:Sz,) (Inf,)

"""
    PQN(N::Real)

The concrete `AbelianNumber` of a quantum system with particle number `N` conserved.
"""
@abeliannumber "PQN" (:N,) (Inf,)

"""
    SPQN(N::Real, Sz::Real)

The concrete `AbelianNumber` of a quantum system with both particle number `N` and spin z-component `Sz` conserved.
"""
@abeliannumber "SPQN" (:N, :Sz) (Inf, Inf)

"""
    SQNS(S::Real) -> AbelianNumbers{SQN}

Construct the `AbelianNumbers` of the Hilbert space of a single spin `S`.
"""
SQNS(S::Real) = AbelianNumbers('C', [SQN(sz) for sz = -S:S], collect(0:Int(2*S+1)), :indptr)

"""
    PQNS(N::Real) -> AbelianNumbers{PQN}

Construct the `AbelianNumbers` of the Hilbert space of a single-particle state with at most `N` identical particles.
"""
PQNS(N::Real) = AbelianNumbers('C', [PQN(np) for np = 0:N], collect(0:Int(N)+1), :indptr)

"""
    SPQNS(S::Real) -> AbelianNumbers{SPQN}

Construct the `AbelianNumbers` of the Hilbert space of a single site with internal degrees of freedom that can be ascribed to a spin `S`.
"""
function SPQNS(S::Real)
    contents = [SPQN(0.0, 0.0)]
    for n = 1:Int(2*S+1)
        for szs in Combinations{n}(-S:S)
            push!(contents, SPQN(n, sum(szs)))
        end
    end
    return sort(AbelianNumbers('G', contents, collect(0:length(contents)), :indptr))[1]
end

"""
    Momentum <: AbelianNumber

Abstract type for momentum.
"""
abstract type Momentum <: AbelianNumber end

"""
    Momentum₁{N}(k::Real) where N

One dimensional momentum.
"""
struct Momentum₁{N} <: Momentum
    k::Float
    function Momentum₁{N}(k::Real) where N
        @assert isa(N, Integer) && N>0 "Momentum₁ error: wrong period ($N)."
        remainder = k % N
        remainder < 0 && (remainder = remainder + N)
        new{N}(remainder)
    end
end
periods(::Type{<:Momentum₁{N}}) where N = (N,)

"""
    Momentum₂{N}(k₁::Real, k₂::Real) where N
    Momentum₂{N₁, N₂}(k₁::Real, k₂::Real) where {N₁, N₂}

Two dimensional momentum.
"""
struct Momentum₂{N₁, N₂} <: Momentum
    k₁::Float
    k₂::Float
    Momentum₂{N}(k₁::Real, k₂::Real) where N =  Momentum₂{N, N}(k₁, k₂)
    function Momentum₂{N₁, N₂}(k₁::Real, k₂::Real) where {N₁, N₂}
        @assert isa(N₁, Integer) && N₁>0 "Momentum₂ error: wrong 1st period ($N₁)."
        @assert isa(N₂, Integer) && N₂>0 "Momentum₂ error: wrong 2nd period ($N₂)."
        r₁ = k₁ % N₁
        r₂ = k₂ % N₂
        r₁ < 0 && (r₁ = r₁ + N₁)
        r₂ < 0 && (r₂ = r₂ + N₂)
        new{N₁, N₂}(r₁, r₂)
    end
end
periods(::Type{<:Momentum₂{N₁, N₂}}) where {N₁, N₂} = (N₁, N₂)

"""
    Momentum₃{N}(k₁::Real, k₂::Real, k₃::Real) where N
    Momentum₃{N₁, N₂, N₃}(k₁::Real, k₂::Real, k₃::Real) where {N₁, N₂, N₃}

Three dimensional momentum.
"""
struct Momentum₃{N₁, N₂, N₃} <: Momentum
    k₁::Float
    k₂::Float
    k₃::Float
    Momentum₃{N}(k₁::Real, k₂::Real, k₃::Real) where N =  Momentum₃{N, N, N}(k₁, k₂, k₃)
    function Momentum₃{N₁, N₂, N₃}(k₁::Real, k₂::Real, k₃::Real) where {N₁, N₂, N₃}
        @assert isa(N₁, Integer) && N₁>0 "Momentum₃ error: wrong 1st period ($N₁)."
        @assert isa(N₂, Integer) && N₂>0 "Momentum₃ error: wrong 2nd period ($N₂)."
        @assert isa(N₃, Integer) && N₃>0 "Momentum₃ error: wrong 3rd period ($N₃)."
        r₁ = k₁ % N₁
        r₂ = k₂ % N₂
        r₃ = k₃ % N₃
        r₁ < 0 && (r₁ = r₁ + N₁)
        r₂ < 0 && (r₂ = r₂ + N₂)
        r₃ < 0 && (r₃ = r₃ + N₃)
        new{N₁, N₂, N₃}(r₁, r₂, r₃)
    end
end
periods(::Type{<:Momentum₃{N₁, N₂, N₃}}) where {N₁, N₂, N₃} = (N₁, N₂, N₃)

end #module


text:
struct ExtractString{T} <: AbstractExtractor
	datatype::Type{T}
	n::Int
	b::Int
	m::Int
end

ExtractString(::Type{T}) where {T<:String} = ExtractString(T, 3, 256, 2053)
(s::ExtractString)(v)   = ArrayNode(Mill.NGramMatrix([v], s.n, s.b, s.m))
(s::ExtractString)(v::S) where {S<:Nothing}= ArrayNode(Mill.NGramMatrix([""], s.n, s.b, s.m))
extractsmatrix(s::ExtractString) = false

function Base.show(io::IO, m::ExtractString;pad = [], key::String="")
	c = COLORS[(length(pad)%length(COLORS))+1]
	key *= isempty(key) ? "" : ": ";
	paddedprint(io,"$(key)$(m.datatype)\n", color = c)
end

"""
	extractscalar(Type{String}, n = 3, b = 256, m = 2053)

	represents strings as ngrams with
		--- `n` (the degree of ngram),
		--- `b` base of string,
		--- `m` modulo on index of the token to reduce dimension

	extractscalar(Type{Number}, m = 0, s = 1)
	extracts number subtracting `m` and multiplying by `s`
"""
extractscalar(::Type{T}, n = 3, b = 256, m = 2053) where {T<:AbstractString} = ExtractString(T, n, b, m)
extractscalar(::Type{T}, m = zero(T), s = one(T)) where {T<:Number} = ExtractScalar(T, m, s)

Base.hash(e::ExtractString, h::UInt) = hash((e.datatype, e.n, e.b, e.m), h)
Base.:(==)(e1::ExtractString, e2::ExtractString) = e1.datatype == e2.datatype && e1.n === e2.n && e1.b === e2.b && e1.m === e2.m


text:
include("discrete.jl")

mutable struct CliffWalkingEnv <: DiscreteEnv
    shape
    start_state_index

    _cliff
    discenv_obj
end

UP = 1
RIGHT = 2
DOWN = 3
LEFT = 4

@inline ind2sub(index, shape) = Int[((index-1) % shape[1])+1, ceil(index/shape[1])]

function CliffWalkingEnv()
    shape = [4, 12]
    start_state_index = 4

    nS = prod(shape)
    nA = 4

    _cliff = zeros(Bool, shape...)
    _cliff[4, 2:11] .= true

    P = Dict()

    """Determine the outcome for an action. Transition Prob is always 1.0
    returns a tuple of form `(1.0, new_state, reward, done)`"""
    function calculate_transition_prob(current, delta)
        new_position = current .+ delta
        limit_coordinates!(new_position)
        new_state = shape[1] * (new_position[2] - 1) + new_position[1]
        _cliff[new_position...] &&
            return [(1.0, start_state_index, -100, false)]
        terminal_state = shape
        is_done = all(new_position .== terminal_state)
        return [(1.0, new_state, -1, is_done)]
    end

    """Don't fall off..."""
    function limit_coordinates!(coord)
        coord[1] = clamp(coord[1], 1, shape[1])
        coord[2] = clamp(coord[2], 1, shape[2])
    end

    P = Dict()
    for s=1:nS
        position = ind2sub(s, shape)
        P[s] = Dict(a => [] for a=1:nA)
        P[s][UP]    = calculate_transition_prob(position, [-1, 0])
        P[s][RIGHT] = calculate_transition_prob(position, [0, 1] )
        P[s][DOWN]  = calculate_transition_prob(position, [1, 0] )
        P[s][LEFT]  = calculate_transition_prob(position, [0, -1])
    end

    isd = zeros(Float32, nS)
    isd[start_state_index] = 1.0

    discenv_obj = DiscreteEnvObj(nS, nA, P, isd)
    CliffWalkingEnv(shape, start_state_index, _cliff, discenv_obj)
end

function render(env::CliffWalkingEnv)
    output = ""
    for y ∈ 1:env.shape[1]
        output_line = ""
        for x ∈ 1:env.shape[2]
            s = env.shape[1]*(x-1) + y
            if env.s == s
                output_line *= " x "
            elseif (y, x) == (4, 12)
                output_line *= " T "
            elseif env._cliff[y, x]
                output_line *= " C "
            else
                output_line *= " o "
            end

            if x == 1
                output_line = lstrip(output_line)
            end
            if x == 12
                output_line = rstrip(output_line)
                output *= '\n'
            end
        end
        output *= output_line
    end
    println(output)
end


text:
using PowerModels, Ipopt, Gurobi
using CSV;
using DataFrames;
using JuMP
include("sol_min.jl");

list_file = readdir("./data/")
normal_data = readdir("./data/normal")

function in_normal_data(normal_data, itema)

    file_name = itema[1:end-2] * ".csv"
    if (file_name in normal_data) == false
        return false
    end

    return true

end

for itema in list_file
    if occursin(".m", itema) && (!in_normal_data(normal_data, itema))

        print(itema)

        case_name = AbstractString[]
        obj_val = Float64[]
        time_limit = Float64[]
        solve_time = Float64[]
        Optimal = []
        dual_feasibility = []
        primal_feasibility = []

        out_file_name = "./data/normal/" * itema[1:end-2] * ".csv"

        network_data = parse_file("./data/" * itema)
        gen_keys = [i for i in keys(network_data["gen"])]

        network_data_new = deepcopy(network_data)
        #relax the inf term
        for i in 1:length(gen_keys)
            if network_data_new["gen"][gen_keys[i]]["qmax"] == Inf
                network_data_new["gen"][gen_keys[i]]["qmax"] = 10^20
            end

            if network_data_new["gen"][gen_keys[i]]["qmax"] == -Inf
                network_data_new["gen"][gen_keys[i]]["qmax"] = -10^20
            end

            if network_data_new["gen"][gen_keys[i]]["qmin"] == Inf
                network_data_new["gen"][gen_keys[i]]["qmin"] = 10^20
            end

            if network_data_new["gen"][gen_keys[i]]["qmin"] == -Inf
                network_data_new["gen"][gen_keys[i]]["qmin"] = -10^20
            end
        end

        if haskey(network_data_new, "dcline")
            dcline_keys = [i for i in keys(network_data["dcline"])]

            for i in 1:length(dcline_keys)
                if network_data_new["dcline"][dcline_keys[i]]["qmaxf"] == Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmaxf"] = 10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qmaxf"] == -Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmaxf"] = -10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qmint"] == Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmint"] = 10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qmint"] == -Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmint"] = -10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qminf"] == Inf
                    network_data_new["dcline"][dcline_keys[i]]["qminf"] = 10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qminf"] == -Inf
                    network_data_new["dcline"][dcline_keys[i]]["qminf"] = -10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qmaxt"] == Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmaxt"] = 10^20
                end

                if network_data_new["dcline"][dcline_keys[i]]["qmaxt"] == -Inf
                    network_data_new["dcline"][dcline_keys[i]]["qmaxt"] = -10^20
                end
            end
        end

        result = run_ac_opf(network_data_new,
            optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3, "max_cpu_time" => 3000.0))

        push!(case_name, itema[1:end-2])
        push!(solve_time, result["solve_time"])
        push!(Optimal, result["termination_status"])
        push!(primal_feasibility, result["primal_status"])
        push!(dual_feasibility, result["dual_status"])
        push!(obj_val, result["objective"])
        push!(time_limit, 3000.0)

        df_result = DataFrame()

        df_result[!, "case_name"] = case_name
        df_result[!, "solve_time"] = solve_time
        df_result[!, "Optimal"] = Optimal
        df_result[!, "primal_feasibility"] = primal_feasibility
        df_result[!, "dual_feasibility"] = dual_feasibility
        df_result[!, "obj_val"] = obj_val
        df_result[!, "time_limit"] = time_limit

        sol_min(result, df_result)

        CSV.write(out_file_name, df_result, bufsize = 2^29)
    end
end

max_stars_repo_name: zhouliang-yu/Grid-ML

max_issues_repo_name: zhouliang-yu/Grid-ML

max_forks_repo_name: zhouliang-yu/Grid-ML


text:
#find out the test statistic when m0=1
import LatentGaussianMixtureModel
import Distributions, StatsBase
import Yeppp

@everywhere using LatentGaussianMixtureModel
@everywhere using Distributions, Yeppp, StatsBase, Random

#Brun calculate the statistic for one data set;
#b is the the random number seed, from 1 to 100
@everywhere function Brun(b::Integer, Ctrue::Integer, Calternative::Integer; debuginfo::Bool=false, ntrials::Int=5, showpower=false)
    nF = 282
    Random.seed!(100)
    n_ij = round.(Int64, rand(Poisson(55), 282).+rand(Exponential(95), 282))
    N = sum(n_ij)

    groupindex = inverse_rle(1:nF, n_ij)
    J=2  #42 #beta dimension
    Random.seed!(200*b)
    X = rand(Distributions.Normal(), N, 2)
    betas_true=ones(J)
    if Ctrue == 1
        mu_true = [log(1/0.779-1)]
        wi_true = [1.0]
        sigmas_true = [0.5]
    elseif Ctrue == 2
        mu_all = log(1/0.779-1)
        mu_true = [mu_all - 2.0, mu_all + 2.0]
        wi_true = [0.5, 0.5]
        sigmas_true = [1.2, 0.8]
    elseif Ctrue == 3
        mu_true = [log(1/0.779 - 1) - 4.0, log(1/0.779 - 1) + 1.0, log(1/0.779 - 1) + 4.0]
        wi_true = [.3, .4, .3]
        sigmas_true = [1.2, .8, .9]
    end
    if showpower
        mu_all = log(1/0.779 - 1)
        if Ctrue == 2
            mu_true = [mu_all - 1.0, mu_all + 0.8]
            wi_true =  [.6, .4]
            sigmas_true = [1.2, .8]
        elseif Ctrue == 3
            mu_true = [mu_all - 2.0, mu_all + 1.0, mu_all + 3.5]
            wi_true = [.3, .4, .3]
            sigmas_true = [1.2, .8, .9]
        end
    end

    #Randomly data generation based on the setting on datagen.jl
    Random.seed!(b * 1000)
    m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
    gamma_true = rand(m, nF)

    prob = exp.(gamma_true[groupindex] .+ X*betas_true)
    prob= prob ./ (1 .+ prob)
    Y = Bool[rand(Binomial(1, prob[i])) == 1 for i in 1:N];
    X = X .- mean(X, dims=1);

    lr=LatentGaussianMixtureModel.EMtest(X, Y, groupindex, Calternative-1, ntrials=ntrials, debuginfo=debuginfo, ctauparallel=false, ngh=100)
    println("Mission $b completed with lr=$lr")
    lr
end


text:
using Gadfly, DataFrames, CategoricalArrays

set_default_plot_size(6inch, 3inch)

a = DataFrame(diff = CategoricalArray([1,2,3,3,3,4,3,2]))
plot(a, x="diff", Geom.histogram)


text:
export Moment, print, ==

function Moment(expo::Exponent, clique::String)
    α, β = Exponent(), Exponent()

    for (var, deg) in expo
        product!(α, Exponent(SortedDict(var=>Degree(0, deg.conjvar))))
        product!(β, Exponent(SortedDict(var=>Degree(deg.explvar, 0))))
    end
    return Moment(α, β, clique)
end

hash(mom::Moment, h::UInt) = hash(mom.conj_part, hash(mom.expl_part, hash(mom.clique, h)))

isless(mom1::Moment, mom2::Moment) = isless((mom1.conj_part, mom1.expl_part, mom1.clique),
                                            (mom2.conj_part, mom2.expl_part, mom2.clique))

function print(io::IO, mom::Moment)
    print(io, "$(mom.expl_part * mom.conj_part) ($(mom.clique))")
end

==(mom1::Moment, mom2::Moment) = ((mom1.conj_part, mom1.expl_part, mom1.clique) == (mom2.conj_part, mom2.expl_part, mom2.clique))


text:
using Optim
using Luxor
using JuMP
using Clp
using Cbc
using GLPK
using Ipopt
using Glob
using Base.Filesystem

include("OptMaps.jl")
using .OptMaps

function solve_model(solver, model_log, solver_options=Dict())
    model = Model(solver)

    for (k, v) in solver_options
        set_optimizer_attribute(model, k, v)
    end

    #sets
    foods = [:brownie, :icecream, :cola, :cheesecake]
    resources = [:calories, :chocolate, :sugar, :fat]

    #parameters
    cost_per_food = [50 20 30 80]
    rsc_per_food = [2000 200 150 500;
                    3 2 0 0;
                    2 2 4 4;
                    2 4 1 5]

    demand_per_rsc = [645 8 10 8]

    n_foods = length(foods)
    n_resources =length(resources)

    n_variables = n_foods
    n_constraints = n_resources

    ###############################################################
    #decision variables
    @variable(model, amt_per_food[1:n_foods] >= 0)

    starting_point = zeros(Float64, n_foods)
    solution_point = zeros(Float64, n_foods)

    for i in 1:n_foods
        v = start_value(amt_per_food[i])
        starting_point[i] = v === nothing ? 0.0 : v
    end

    ###############################################################
    #objective

    #regular function
    total_cost_f(amt_per_food) = sum([cost_per_food[f] * amt_per_food[f] for f
                                      in 1:n_foods])

    #jump expression
    total_cost = @expression(model, total_cost, sum(cost_per_food[f] *
                                                    amt_per_food[f] for f in
                                                        1:n_foods))
    @objective(model, Min, total_cost)

    ###############################################################
    #constraints: fulfill min nutrition demand (demand_per_rsc)

    #regular functions
    constraints_f = []

    for r in 1:n_resources
        push!(constraints_f, (amt_per_food) ->
              sum([amt_per_food[f]*rsc_per_food[r,f]  for f in 1:n_foods]) >=
              demand_per_rsc[r])
    end

    #jump expressions
    @constraint(model, dmd[r = 1:n_resources],
                sum((amt_per_food[f]*rsc_per_food[r,f]  for f in 1:n_foods)) >=
                demand_per_rsc[r])

    ###############################################################
    #save model textual description
    open(model_log, "w") do f
        print(f, model)

        ###############################################################
        #solve model
        elapsed_time = @elapsed optimize!(model)

        println(f, "-----------------------------------------------")
        println(f, "Elapsed time: $elapsed_time")
        println(f, "-----------------------------------------------")
        println(f, "Status:    $(termination_status(model))")
        println(f, "Objective: $(objective_value(model))")
        println(f, "-----------------------------------------------")
        println(f, "Variables:")

        for fd in 1:n_foods
            println(f, "$(foods[fd]): $(value(amt_per_food[fd]))")
            solution_point[fd] = value(amt_per_food[fd])
        end
        # println(f, "-----------------------------------------------")
        # println(f, "Shadow prices:")
        # for r in 1:n_resources
        #     println(f, "$r - $(resources[r]):  $(shadow_price(dmd[r]))")
        # end
    end

    variable_names = ["amt_per_food[$(String(f))]" for f in foods]
    constraint_names = ["dmd[$(String(r))]" for r in resources]

    obj_val = objective_value(model)

    return total_cost_f, constraints_f, variable_names, constraint_names,
    starting_point, solution_point, obj_val
end

#The datatype used in the tuples is used to determine if the variable is "real"
#or integer
variable_domains = [(0.0,4.0), (0.0,4.0), (0.0,4.0), (0.0,4.0)]

optimizers = [GLPK.Optimizer, Ipopt.Optimizer, Cbc.Optimizer]
opt_names = ["glpk", "ipopt", "cbc"]

for (opt, opt_name) in zip(optimizers, opt_names)
    objective_f, constraints_f, variable_names, constraint_names,
    starting_point, solution_point, obj_val = solve_model(opt,
                                       "test_optmap_$(opt_name).log")
end

opt = Clp.Optimizer
opt_name = "clp"

objective_f, constraints_f, variable_names, constraint_names, starting_point,
solution_point, obj_val = solve_model(opt, "test_optmap_$(opt_name).log")

path_to_solution = zeros(Float64, length(variable_domains), 2)
path_to_solution[:,1] = starting_point
path_to_solution[:,2] = solution_point

O = OptMap(variable_names, variable_domains, objective_f;
           constraint_names=constraint_names, constraints_f=constraints_f,
           verbose=true)
create_map!(O)

draw_solution!(O, path_to_solution)
save_all(O, "test_optmap_$(opt_name)", ".")

test_files = glob("test_optmap*.*", ".")

#Precompile Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
initial_x = [0.0, 0.0]

solvers = [NelderMead, LBFGS, GradientDescent, SimulatedAnnealing,
           ConjugateGradient, Newton]

for solver in solvers
    result = optimize(f, initial_x, solver(), Optim.Options(store_trace = true,
                                                            extended_trace=true,
                                                            iterations=50000))

    if solver == NelderMead
        tr = Optim.trace(result)
        points = [state.metadata["centroid"] for state in tr]
    else
        points = Optim.x_trace(result)
    end
end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
using Zygote
using Statistics

# We modify (the implementation of) batchnorm to be more ammenable to TPUs
struct TPUBatchNorm{F,V,W}
   λ::F  # activation function
   β::V  # bias
   γ::V  # scale
   μ::W  # moving mean
   σ::W  # moving std
   ϵ::XRTArray{Float32, (), 0}
   momentum::XRTArray{Float32, (), 0}
   #active::Bool
end

Flux.children(BN::TPUBatchNorm) =
  (BN.λ, BN.β, BN.γ, BN.μ, BN.σ, BN.ϵ, BN.momentum)

Flux.mapchildren(f, BN::TPUBatchNorm) =  # e.g. mapchildren(cu, BN)
  TPUBatchNorm(BN.λ, f(BN.β), f(BN.γ), f(BN.μ), f(BN.σ), BN.ϵ, BN.momentum #=, BN.active=#)

function compute_affine_shape(x)
  dims = length(size(x))
  channels = size(x, dims-1)
  affine_shape = let dims=dims, channels=channels
   ntuple(i->i == dims - 1 ? channels : 1, dims)
  end
  #=
  m = let sz = size(x)
    prod(ntuple(i->sz[i], dims-2)) * sz[end]
  end
  =#
  affine_shape
end
@Zygote.adjoint compute_affine_shape(x) = compute_affine_shape(x), Δ->(nothing,)

# This is a bit of a trick. We use Zygote's backward pass infrastructure
# to backpropagate the batchnorm statistics to the parameter update
function tpu_batchnorm_statistics(active::Bool, bn_μ, bn_σ, bn_ε, x)
  #@assert !BN.active
  affine_shape = compute_affine_shape(x)
  μ = reshape(bn_μ, affine_shape...)
  σ = reshape(bn_σ, affine_shape...)
  return (x .- μ)./σ
end

@Zygote.adjoint function tpu_batchnorm_statistics(active::Bool, bn_μ, bn_σ, bn_ϵ, x)
  ϵ = convert(eltype(x), bn_ϵ)
  axes = (1, 2, 4)
  m = Int32(size(x, 1) * size(x, 2) * size(x, 4))

  # Calculate μ and σ for the "forward pass"
  μ = mean(x, dims = axes)
  meansub = (x .- μ)
  sq = meansub .* meansub
  σ = sqrt.(mean(sq, dims = axes) .+ ϵ)

  x̂ = (x .- μ)./σ
  # Return "forward pass" calculation and closure that returns our "gradients"
  x̂, let μ_dropped=dropdims(μ, dims = axes),
         σ_dropped=dropdims(σ, dims = axes) .* convert(eltype(x), m) / convert(eltype(x), m - 1)
    function(Δ)
      let axes = (1, 2, 4), x̂ = x̂
        # https://kratzert.github.io/2016/02/12/understanding-the-gradient-flow-through-the-batch-normalization-layer.html
        Δ_x = (Δ .- mean(Δ, dims=axes) .- x̂ .* mean(Δ .* x̂, dims=axes))./σ
        return (nothing, μ_dropped, σ_dropped, nothing, Δ_x)
      end
    end
  end
end

function (BN::TPUBatchNorm)(x)
  γ = BN.γ
  β = BN.β
  affine_shape = compute_affine_shape(x)

  x̂ = tpu_batchnorm_statistics(true, BN.μ, BN.σ, BN.ϵ, x)

  let λ = BN.λ
    λ.(reshape(γ, affine_shape...) .* x̂ .+ reshape(β, affine_shape...))
  end
end

function update_params(BN::TPUBatchNorm{F,V,W}, updates, state, η, β) where {F,V,W}
    mtm = BN.momentum
    new_λ, new_λ_state = update_params(BN.λ, updates.λ, state[1], η, β)
    new_β, new_β_state = update_params(BN.β, updates.β, state[2], η, β)
    new_γ, new_γ_state = update_params(BN.γ, updates.γ, state[3], η, β)

    new_BN = TPUBatchNorm{F,V,W}(
      new_λ,
      new_β,
      new_γ,
      (XRTArray(1f0) - mtm) .* BN.μ + mtm .* updates.μ,
      (XRTArray(1f0) - mtm) .* BN.σ + mtm .* updates.σ,
      BN.ϵ, mtm #, BN.active
    )
    new_state = tuple(
      new_λ_state,
      new_β_state,
      new_γ_state,
      state[4],
      state[5],
      state[6],
      state[7],
    )
    return (new_BN, new_state)
end


text:
using WindMixing
using OceanParameterizations

PATH = joinpath(pwd(), "extracted_training_output")
FILE_PATH = joinpath(pwd(), "Output")

train_files = [
                "wind_-5e-4_diurnal_5e-8"
                "wind_-5e-4_diurnal_3.5e-8"
                "wind_-5e-4_diurnal_3e-8"
                "wind_-5e-4_diurnal_2e-8"
                "wind_-5e-4_diurnal_1e-8"
                "wind_-3.5e-4_diurnal_5e-8"
                "wind_-3.5e-4_diurnal_3.5e-8"
                "wind_-3.5e-4_diurnal_3e-8"
                "wind_-3.5e-4_diurnal_2e-8"
                "wind_-3.5e-4_diurnal_1e-8"

                "wind_-2e-4_diurnal_5e-8"
                "wind_-2e-4_diurnal_3.5e-8"
                "wind_-2e-4_diurnal_3e-8"
                "wind_-2e-4_diurnal_2e-8"
                "wind_-2e-4_diurnal_1e-8"
               ]

# train_files = ["wind_-1e-3_cooling_4e-8",
#                "wind_-2e-4_cooling_1e-8"]
# train_files = ["wind_-2e-4_cooling_1e-8",]
𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

VIDEO_NAME = "LES_simulations_test_2"

animate_training_data_profiles_fluxes(train_files, joinpath(FILE_PATH, VIDEO_NAME))


text:
int_rules_1_2_3_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.3.1*(a+b*x^n+c*x^(2*n))^p =#
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (2 * ~n * ~p) * (~c + ~b * (~x) ^ -(~n) + ~a * (~x) ^ (-2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (LtQ(~n, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * (~x) ^ (k * ~n) + ~c * (~x) ^ (2 * k * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && FractionQ(~n))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -(Subst(Antiderivative((~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p)) * Antiderivative((~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IGtQ(~p, 0)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (-(~x) * (((~b) ^ 2 - 2 * ~a * ~c) + ~b * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((((~b) ^ 2 - 2 * ~a * ~c) + ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) + ~b * ~c * (~n * (2 * ~p + 3) + 1) * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && ILtQ(~p, -1)))
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt(~a / ~c, 2)], With([r = Rt(2q - ~b / ~c, 2)], (1 / (2 * ~c * q * r)) * Antiderivative((r - (~x) ^ (~n / 2)) / ((q - r * (~x) ^ (~n / 2)) + (~x) ^ ~n), ~x) + (1 / (2 * ~c * q * r)) * Antiderivative((r + (~x) ^ (~n / 2)) / (q + r * (~x) ^ (~n / 2) + (~x) ^ ~n), ~x))) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n / 2, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))))
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~c / q) * Antiderivative(1 / ((~b / 2 - q / 2) + ~c * (~x) ^ ~n), ~x) - (~c / q) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) / ((1 + (2 * ~c * (~x) ^ ~n) / (~b + Rt((~b) ^ 2 - 4 * ~a * ~c, 2))) ^ FracPart(~p) * (1 + (2 * ~c * (~x) ^ ~n) / (~b - Rt((~b) ^ 2 - 4 * ~a * ~c, 2))) ^ FracPart(~p))) * Antiderivative((1 + (2 * ~c * (~x) ^ ~n) / (~b + sqrt((~b) ^ 2 - 4 * ~a * ~c))) ^ ~p * (1 + (2 * ~c * (~x) ^ ~n) / (~b - sqrt((~b) ^ 2 - 4 * ~a * ~c))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~a + ~(b') * (~u) ^ ~n + ~(c') * (~u) ^ ~(n2')) ^ ~p, ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (LinearQ(~u, ~x) && NeQ(~u, ~x)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative((~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (EqQ(~mn, -(~n)) && (IntegerQ(~p) && PosQ(~n)))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~p, ~x) => (((~x) ^ (~n * FracPart(~p)) * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ ~n) ^ FracPart(~p)) / (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~mn, -(~n)) && (Not(IntegerQ(~p)) && PosQ(~n)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.3 General/1.2.3.1 (a+b x^n+c x^(2 n))^p.jl