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

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

Associates your quota request template with your organization. When a new account is
created in your organization, the quota increase requests in the template are automatically
applied to the account. You can add a quota increase request for any adjustable quota to
your template.

"""
function associate_service_quota_template(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "AssociateServiceQuotaTemplate";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function associate_service_quota_template(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "AssociateServiceQuotaTemplate",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_service_quota_increase_request_from_template(aws_region, quota_code, service_code)
    delete_service_quota_increase_request_from_template(aws_region, quota_code, service_code, params::Dict{String,<:Any})

Deletes the quota increase request for the specified quota from your quota request template.

# Arguments
- `aws_region`: The AWS Region.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function delete_service_quota_increase_request_from_template(
    AwsRegion, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "DeleteServiceQuotaIncreaseRequestFromTemplate",
        Dict{String,Any}(
            "AwsRegion" => AwsRegion, "QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_service_quota_increase_request_from_template(
    AwsRegion,
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "DeleteServiceQuotaIncreaseRequestFromTemplate",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AwsRegion" => AwsRegion,
                    "QuotaCode" => QuotaCode,
                    "ServiceCode" => ServiceCode,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Disables your quota request template. After a template is disabled, the quota increase
requests in the template are not applied to new accounts in your organization. Disabling a
quota request template does not apply its quota increase requests.

"""
function disassociate_service_quota_template(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "DisassociateServiceQuotaTemplate";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function disassociate_service_quota_template(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "DisassociateServiceQuotaTemplate",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the status of the association for the quota request template.

"""
function get_association_for_service_quota_template(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetAssociationForServiceQuotaTemplate";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_association_for_service_quota_template(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetAssociationForServiceQuotaTemplate",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_awsdefault_service_quota(quota_code, service_code)
    get_awsdefault_service_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the default value for the specified quota. The default value does not reflect any
quota increases.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function get_awsdefault_service_quota(
    QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetAWSDefaultServiceQuota",
        Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_awsdefault_service_quota(
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "GetAWSDefaultServiceQuota",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves information about the specified quota increase request.

# Arguments
- `request_id`: The ID of the quota increase request.

"""
function get_requested_service_quota_change(
    RequestId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetRequestedServiceQuotaChange",
        Dict{String,Any}("RequestId" => RequestId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_requested_service_quota_change(
    RequestId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "GetRequestedServiceQuotaChange",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("RequestId" => RequestId), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_service_quota(quota_code, service_code)
    get_service_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the applied quota value for the specified quota. For some quotas, only the
default values are available. If the applied quota value is not available for a quota, the
quota is not retrieved.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function get_service_quota(
    QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetServiceQuota",
        Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_service_quota(
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "GetServiceQuota",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_service_quota_increase_request_from_template(aws_region, quota_code, service_code)
    get_service_quota_increase_request_from_template(aws_region, quota_code, service_code, params::Dict{String,<:Any})

Retrieves information about the specified quota increase request in your quota request
template.

# Arguments
- `aws_region`: The AWS Region.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function get_service_quota_increase_request_from_template(
    AwsRegion, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "GetServiceQuotaIncreaseRequestFromTemplate",
        Dict{String,Any}(
            "AwsRegion" => AwsRegion, "QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_service_quota_increase_request_from_template(
    AwsRegion,
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "GetServiceQuotaIncreaseRequestFromTemplate",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AwsRegion" => AwsRegion,
                    "QuotaCode" => QuotaCode,
                    "ServiceCode" => ServiceCode,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the default values for the quotas for the specified AWS service. A default value does
not reflect any quota increases.

# Arguments
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
function list_awsdefault_service_quotas(
    ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListAWSDefaultServiceQuotas",
        Dict{String,Any}("ServiceCode" => ServiceCode);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_awsdefault_service_quotas(
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "ListAWSDefaultServiceQuotas",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ServiceCode" => ServiceCode), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the quota increase requests for the specified service.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"ServiceCode"`: The service identifier.
- `"Status"`: The status of the quota increase request.
"""
function list_requested_service_quota_change_history(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListRequestedServiceQuotaChangeHistory";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_requested_service_quota_change_history(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListRequestedServiceQuotaChangeHistory",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_requested_service_quota_change_history_by_quota(quota_code, service_code)
    list_requested_service_quota_change_history_by_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the quota increase requests for the specified quota.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"Status"`: The status value of the quota increase request.
"""
function list_requested_service_quota_change_history_by_quota(
    QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListRequestedServiceQuotaChangeHistoryByQuota",
        Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_requested_service_quota_change_history_by_quota(
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "ListRequestedServiceQuotaChangeHistoryByQuota",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("QuotaCode" => QuotaCode, "ServiceCode" => ServiceCode),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the quota increase requests in the specified quota request template.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AwsRegion"`: The AWS Region.
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"ServiceCode"`: The service identifier.
"""
function list_service_quota_increase_requests_in_template(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListServiceQuotaIncreaseRequestsInTemplate";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_service_quota_increase_requests_in_template(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListServiceQuotaIncreaseRequestsInTemplate",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the applied quota values for the specified AWS service. For some quotas, only the
default values are available. If the applied quota value is not available for a quota, the
quota is not retrieved.

# Arguments
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
function list_service_quotas(ServiceCode; aws_config::AbstractAWSConfig=global_aws_config())
    return service_quotas(
        "ListServiceQuotas",
        Dict{String,Any}("ServiceCode" => ServiceCode);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_service_quotas(
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "ListServiceQuotas",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ServiceCode" => ServiceCode), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the names and codes for the services integrated with Service Quotas.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
function list_services(; aws_config::AbstractAWSConfig=global_aws_config())
    return service_quotas(
        "ListServices"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_services(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListServices", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

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

Returns a list of the tags assigned to the specified applied quota.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota for which you want
  to list tags. You can get this information by using the Service Quotas console, or by
  listing the quotas using the list-service-quotas AWS CLI command or the ListServiceQuotas
  AWS API operation.

"""
function list_tags_for_resource(
    ResourceARN; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "ListTagsForResource",
        Dict{String,Any}("ResourceARN" => ResourceARN);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    ResourceARN,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "ListTagsForResource",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ResourceARN" => ResourceARN), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_service_quota_increase_request_into_template(aws_region, desired_value, quota_code, service_code)
    put_service_quota_increase_request_into_template(aws_region, desired_value, quota_code, service_code, params::Dict{String,<:Any})

Adds a quota increase request to your quota request template.

# Arguments
- `aws_region`: The AWS Region.
- `desired_value`: The new, increased value for the quota.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function put_service_quota_increase_request_into_template(
    AwsRegion,
    DesiredValue,
    QuotaCode,
    ServiceCode;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "PutServiceQuotaIncreaseRequestIntoTemplate",
        Dict{String,Any}(
            "AwsRegion" => AwsRegion,
            "DesiredValue" => DesiredValue,
            "QuotaCode" => QuotaCode,
            "ServiceCode" => ServiceCode,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_service_quota_increase_request_into_template(
    AwsRegion,
    DesiredValue,
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "PutServiceQuotaIncreaseRequestIntoTemplate",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AwsRegion" => AwsRegion,
                    "DesiredValue" => DesiredValue,
                    "QuotaCode" => QuotaCode,
                    "ServiceCode" => ServiceCode,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    request_service_quota_increase(desired_value, quota_code, service_code)
    request_service_quota_increase(desired_value, quota_code, service_code, params::Dict{String,<:Any})

Submits a quota increase request for the specified quota.

# Arguments
- `desired_value`: The new, increased value for the quota.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
function request_service_quota_increase(
    DesiredValue, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "RequestServiceQuotaIncrease",
        Dict{String,Any}(
            "DesiredValue" => DesiredValue,
            "QuotaCode" => QuotaCode,
            "ServiceCode" => ServiceCode,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function request_service_quota_increase(
    DesiredValue,
    QuotaCode,
    ServiceCode,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "RequestServiceQuotaIncrease",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "DesiredValue" => DesiredValue,
                    "QuotaCode" => QuotaCode,
                    "ServiceCode" => ServiceCode,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Adds tags to the specified applied quota. You can include one or more tags to add to the
quota.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota. You can get this
  information by using the Service Quotas console, or by listing the quotas using the
  list-service-quotas AWS CLI command or the ListServiceQuotas AWS API operation.
- `tags`: The tags that you want to add to the resource.

"""
function tag_resource(ResourceARN, Tags; aws_config::AbstractAWSConfig=global_aws_config())
    return service_quotas(
        "TagResource",
        Dict{String,Any}("ResourceARN" => ResourceARN, "Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    ResourceARN,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "TagResource",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("ResourceARN" => ResourceARN, "Tags" => Tags),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Removes tags from the specified applied quota. You can specify one or more tags to remove.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota that you want to
  untag. You can get this information by using the Service Quotas console, or by listing the
  quotas using the list-service-quotas AWS CLI command or the ListServiceQuotas AWS API
  operation.
- `tag_keys`: The keys of the tags that you want to remove from the resource.

"""
function untag_resource(
    ResourceARN, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return service_quotas(
        "UntagResource",
        Dict{String,Any}("ResourceARN" => ResourceARN, "TagKeys" => TagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    ResourceARN,
    TagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return service_quotas(
        "UntagResource",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("ResourceARN" => ResourceARN, "TagKeys" => TagKeys),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
## Defines all hashing functionality
## As usual, check out http://www.lysator.liu.se/~nisse/nettle/nettle.html#Hash-functions

import Base: show
export Hasher, update!, digest, digest!, hexdigest!, hexdigest

immutable Hasher
    hash_type::HashType
    state::Vector{UInt8}
end

# Constructor for Hasher
function Hasher(name::AbstractString)
    hash_types = get_hash_types()
    name = uppercase(name)
    if !haskey(hash_types, name)
        throw(ArgumentError("Invalid hash type $name: call Nettle.get_hash_types() to see available list"))
    end

    # Construct Hasher object for this type and initialize using Nettle's init functions
    hash_type = hash_types[name]
    state = Vector{UInt8}(hash_type.context_size)
    ccall(hash_type.init, Void, (Ptr{Void},), state)
    return Hasher(hash_type, state)
end

# Update hash state with new data
function update!(state::Hasher, data)
    ccall(state.hash_type.update, Void, (Ptr{Void},Csize_t,Ptr{UInt8}), state.state, sizeof(data), pointer(data))
    return state
end

# Spit out a digest of the current hash state and reset it
function digest!(state::Hasher)
    digest = Vector{UInt8}(state.hash_type.digest_size)
    ccall(state.hash_type.digest, Void, (Ptr{Void},UInt32,Ptr{UInt8}), state.state, sizeof(digest), pointer(digest))
    return digest
end

# Take a digest, and convert it to a printable hex representation
hexdigest!(state::Hasher) = bytes2hex(digest!(state))

# The one-shot functions that makes this whole thing so easy.
digest(hash_name::AbstractString, data) = digest!(update!(Hasher(hash_name), data))
digest(hash_name::AbstractString, io::IO) = digest(hash_name, readall(io))
hexdigest(hash_name::AbstractString, data) = hexdigest!(update!(Hasher(hash_name), data))
hexdigest(hash_name::AbstractString, io::IO) = hexdigest(hash_name, readall(io))

# Custom show overrides make this package have a little more pizzaz!
function show(io::IO, x::HashType)
    write(io, "$(x.name) Hash\n")
    write(io, "  Context size: $(x.context_size) bytes\n")
    write(io, "  Digest size: $(x.digest_size) bytes\n")
    write(io, "  Block size: $(x.block_size) bytes")
end
show(io::IO, x::Hasher) = write(io, "$(x.hash_type.name) Hash state")

ext: jl


text:
using BenchmarkTools
using Arrow
using CategoricalArrays

data = convert(Vector{UInt8}, reinterpret(UInt8, [2,3,5,7]))
data = vcat(0x0f, data)

# v = NullablePrimitive([2,missing,5,7])

v = NullablePrimitive{Int64}(data, 1, 2, 4)


text:
# shortestpath.jl
# ===============
#
# A DS_DN_Counter for use with the DS_DN_enumerator in DS_DN_enumerator.jl
# Only keeps and returns the DS and DN values for
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/NaturalSelection.jl/blob/master/LICENSE.md

struct ShortestPath <: DS_DN_Counter end
@inline DS_DN_init(::Type{ShortestPath}) = 0.0, 3.0
@inline weighting(::Type{ShortestPath}, ::Val{T}) where T = 1.0
@inline function accumulate(::Type{ShortestPath}, DS::Float64, DN::Float64, DS_i::Float64, DN_i::Float64)
    switchout = DN_i < DN
    DN = ifelse(switchout, DN_i, DN)
    DS = ifelse(switchout, DS_i, DS)
    return DS, DN
end


text:
## %PREDICTMD_GENERATED_BY%

import PredictMDExtra
PredictMDExtra.import_all()

import PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "boston_housing", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin linear regression code

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

categorical_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "categorical_feature_names.jld2",
    )
continuous_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "continuous_feature_names.jld2",
    )
categorical_feature_names = FileIO.load(
    categorical_feature_names_filename,
    "categorical_feature_names",
    )
continuous_feature_names = FileIO.load(
    continuous_feature_names_filename,
    "continuous_feature_names",
    )
feature_names = vcat(categorical_feature_names, continuous_feature_names)

single_label_name = :MedV
continuous_label_names = Symbol[single_label_name]
categorical_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

# PREDICTMD IF INCLUDE TEST STATEMENTS
Test.@test(
    PredictMD.columns_are_linearly_independent(training_features_df)
    )
Test.@test(
    PredictMD.columns_are_linearly_independent(
        training_features_df,
        feature_names,
        )
    )
Test.@test(
    length(PredictMD.linearly_dependent_columns(df)) == 0
    )
Test.@test(
    length(
        PredictMD.linearly_dependent_columns(
            training_features_df,
            feature_names,
            ),
        ) == 0
    )
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

show(
    logger_stream, PredictMD.linearly_dependent_columns(df)
    )

show(
    logger_stream, PredictMD.linearly_dependent_columns(
        training_features_df,
        feature_names,
        )
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
Test.@test_throws ErrorException PredictMD.single_labeldataframelinearregression(
    feature_names,
    single_label_name;
    package = :thispackagedoesnotexist,
    intercept = true,
    interactions = 1,
    name = "Linear regression",
    )
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

linear_regression = PredictMD.single_labeldataframelinearregression(
    feature_names,
    single_label_name;
    package = :GLM,
    intercept = true,
    interactions = 1,
    name = "Linear regression",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
DEEPCOPY_training_features_df = deepcopy(training_features_df)
DEEPCOPY_linear_regression = deepcopy(linear_regression)
DEEPCOPY_training_labels_df = deepcopy(training_labels_df)
double_DEEPCOPY_training_features_df = hcat(DEEPCOPY_training_features_df, DEEPCOPY_training_features_df; makeunique=true)
PredictMD.fit!(DEEPCOPY_linear_regression,double_DEEPCOPY_training_features_df,training_labels_df)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.predict(linear_regression, training_features_df)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

PredictMD.fit!(linear_regression,training_features_df,training_labels_df)

PredictMD.get_underlying(linear_regression) # TODO: fix this error

linear_regression_plot_training =
    PredictMD.plotsinglelabelregressiontrueversuspredicted(
        linear_regression,
        training_features_df,
        training_labels_df,
        single_label_name,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "linear_regression_plot_training",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, linear_regression_plot_training)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(linear_regression_plot_training)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "linear_regression_plot_training.pdf",
        ),
    linear_regression_plot_training,
    )

linear_regression_plot_testing =
    PredictMD.plotsinglelabelregressiontrueversuspredicted(
        linear_regression,
        testing_features_df,
        testing_labels_df,
        single_label_name
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "linear_regression_plot_testing",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, linear_regression_plot_testing)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(linear_regression_plot_testing)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "linear_regression_plot_testing.pdf",
        ),
    linear_regression_plot_testing,
    )

show(
    logger_stream, PredictMD.singlelabelregressionmetrics(
        linear_regression,
        training_features_df,
        training_labels_df,
        single_label_name,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

show(
    logger_stream, PredictMD.singlelabelregressionmetrics(
        linear_regression,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

linear_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "linear_regression.jld2",
    )

PredictMD.save_model(linear_regression_filename, linear_regression)

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.predict(linear_regression, training_features_df)
Test.@test_throws ErrorException PredictMD.predict_proba(linear_regression, training_features_df)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
linear_regression_filename_bson = joinpath(
    PredictMD.maketempdir(),
    "linear_regression.bson",
    )
PredictMD.save_model(
    linear_regression_filename_bson,
    linear_regression,
    )
test_load_bson = PredictMD.load_model(
    linear_regression_filename_bson,
    )
Test.@test_throws(
    ErrorException,
    PredictMD.save_model("test.nonexistentextension", linear_regression)
    )
Test.@test_throws(
    ErrorException,
    PredictMD.save_model_jld2("test.nonexistentextension", linear_regression)
    )
Test.@test_throws(
    ErrorException,
    PredictMD.save_model_bson("test.nonexistentextension", linear_regression)
    )
Test.@test_throws(
    ErrorException,
    PredictMD.load_model("test.nonexistentextension")
    )
Test.@test_throws(
    ErrorException,
    PredictMD.load_model_jld2("test.nonexistentextension")
    )
Test.@test_throws(
    ErrorException,
    PredictMD.load_model_bson("test.nonexistentextension")
    )
linear_regression = nothing
Test.@test isnothing(linear_regression)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End linear regression code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "boston_housing", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
module libcsfml
	using CBinding

	# NOTE: Libglvnd_jll and GLU_jll are only available for Linux
	@static if Sys.islinux()
		@eval using Libglvnd_jll
		@eval using GLU_jll
		@eval using CSFML_jll

		@eval c`$([
			Main.SYSROOT...,

			"-I$(Libglvnd_jll.artifact_dir)/include",
			"-L$(dirname(Libglvnd_jll.libGL_path))", "-lGL",

			"-I$(GLU_jll.artifact_dir)/include",
			"-L$(dirname(GLU_jll.libGLU_path))", "-lGLU",

			"-I$(CSFML_jll.artifact_dir)/include",
			"-L$(dirname(CSFML_jll.libcsfml_system_path))", "-lcsfml-system",
			"-L$(dirname(CSFML_jll.libcsfml_audio_path))", "-lcsfml-audio",
			"-L$(dirname(CSFML_jll.libcsfml_graphics_path))", "-lcsfml-graphics",
			"-L$(dirname(CSFML_jll.libcsfml_network_path))", "-lcsfml-network",
			"-L$(dirname(CSFML_jll.libcsfml_window_path))", "-lcsfml-window",
		])`

		@eval const c"size_t" = Csize_t
		@eval const c"wchar_t" = Cwchar_t

		@eval c"""
			#include <SFML/System.h>
			#include <SFML/Config.h>
			#include <SFML/GPUPreference.h>
			#include <SFML/Graphics.h>
			#include <SFML/Network.h>
			#include <SFML/OpenGL.h>
			#include <SFML/System.h>
			#include <SFML/Window.h>
		"""ji
	end
end


text:
module SemidefiniteProgramming

import
    Base.print,
    Base.show,
    Base.getindex,
    Base.setindex!,
    Base.size,
    Base.copy,
    Base.(*),
    Base.delete!,
    Base.(==),
    Base.isapprox
    # Base.start,# This is deprciated
    # Base.next, # This is deprciated and is now called Base.nextind but not the same
    # Base.done, # This is deprciated

# export
    # SparseSDP,
    # ismaximizationproblem,
    # setmaximizationproblem!,
    # setobj!,
    # obj,
    # setrhs!,
    # rhs,
    # setcon!,
    # cons,
    # ncons,
    # blocksizes,
    #
    # writesdpasparse,
    # readsdpasparse,
    #
    # SparseSDPSolution,
    # dualobj,
    # primalmatrix,
    # dualvector,
    # dualmatrix,
    #
    # SDPSolver,
    # SDPA,
    # SDPAQD,
    # SDPAGMP,
    # CSDP,
    # solve

export
 # sparsesymmetricmatrix.jl
    sanityCheck,
    SparseSymmetricMatrix,
    setindex!,
    getindex,
    size,
    *,
    delete!,
    ==,
    isapprox, # This one still gives errors.
# sparsesymmetricblockmatrix.jl
    SparseSymmetricBlockMatrix,
    blocks,
    setindex!,
    getindex,
    copy,
    *,
    show,
    ==,
    isapprox,
# sparsesdp.jl
    SparseSDP,
    ismaximizationproblem,
    setmaximizationproblem!,
    normalized,
    copy,
    setobj!,
    obj,
    freeobj,
    setcon!,
    setfreecon!,
    cons,
    getcon,
    setrhs!,
    rhs,
    rhsdense,
    ncons,
    blocksizes
# normalization.jl
    sanitycheck
    # IndexMap,
    # newindex!
    # normalize

include("sparsesymmetricmatrix.jl")
include("sparsesymmetricblockmatrix.jl")
include("sparsesdp.jl")
include("normalization.jl")
# include("sdpasparse.jl")
# include("solvers.jl")
# include("sdpsolution.jl")

end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
# driver for benchmarking different load balancing algorithms

using DataFrames
using CSV

P = 2
addprocs(P)

@everywhere using LoadBalancers

header = [:strategy, :nprocs, :nwork, :walltime, :overhead, :speedup, :efficiency]

df = ["strategy" 1 1 1 1 1 1;]

W = 20
Uᵧ = 3
Tᵧ = W*Uᵧ

@printf("ARR Strategy\n")
work = WorkUnit(W, Uᵧ)
Tₚ = @elapsed parallel_lb_arr(work)

Tₒ = P*Tₚ - Tᵧ
S = Tᵧ/Tₚ
E = S/P

df = [df; "arr" P W Tₚ Tₒ S E]
#@printf("Walltime: %4.2f\nOverhead/P: %4.2f\nSpeedup: %4.2f\nEfficiency: %4.2f\n", Tₚ, Tₒ/P, S, E)

@printf("RP Strategy\n")
work = WorkUnit(W, Uᵧ)
Tₚ = @elapsed parallel_lb_rp(work)

Tₒ = P*Tₚ - Tᵧ
S = Tᵧ/Tₚ
E = S/P

#@printf("Walltime: %4.2f\nOverhead/P: %4.2f\nSpeedup: %4.2f\nEfficiency: %4.2f\n", Tₚ, Tₒ/P, S, E)
df = [df; "rp" P W Tₚ Tₒ S E]

@printf("SB Strategy\n")
work = WorkUnit(W, Uᵧ)
Tₚ = @elapsed parallel_lb_sb(work)

Tₒ = P*Tₚ - Tᵧ
S = Tᵧ/Tₚ
E = S/P

#@printf("Walltime: %4.2f\nOverhead/P: %4.2f\nSpeedup: %4.2f\nEfficiency: %4.2f\n", Tₚ, Tₒ/P, S, E)
df = [df; "sb" P W Tₚ Tₒ S E]

df = DataFrame(df[2:end,1:end], header)

display(df)

CSV.write("out.csv", df)

max_stars_repo_name: dillondaudert/jlance

max_issues_repo_name: dillondaudert/jlance

max_forks_repo_name: dillondaudert/jlance


text:
#---------------------------------------------------------
# # [ImageGeoms overview](@id 1-overview)
#---------------------------------------------------------

# This page explains the Julia package
# [`ImageGeoms`](https://github.com/JuliaImageRecon/ImageGeoms.jl).

# ### Setup

# Packages needed here.

using ImageGeoms
using MIRTjim: jim, prompt
using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# When performing tomographic image reconstruction,
# one must specify the geometry of the grid of image pixels.
# (In contrast, for image denoising and image deblurring problems,
# one works with the given discrete image
# and no physical coordinates are needed.)

# The key parameters of a grid of image pixels are
# * the size (dimensions) of the grid, e.g., `128 × 128`,
# * the spacing of the pixels, e.g., `1mm × 1mm`,
# * the offset of the pixels relative to the origin, e.g., `(0,0)`

# The data type `ImageGeom` describes such a geometry,
# for arbitrary dimensions (2D, 3D, etc.).

# There are several ways to construct this structure.
# The default is a `128 × 128` grid
# with pixel size ``\Delta_X = \Delta_Y = 1`` (unitless) and zero offset:

ig = ImageGeom()

# Here is a 3D example with non-cubic voxel size:

ig = ImageGeom( (512,512,128), (1,1,2), (0,0,0) )

# To avoid remembering the order of the arguments,
# named keyword pairs are also supported:

ig = ImageGeom( dims=(512,512,128), deltas=(1,1,2), offsets=(0,0,0) )

# ### Units

# The pixel dimensions `deltas` can (and should!) be values with units.

# Here is an example for a video (2D+time) with 12 frames per second:
using UnitfulRecipes
using Unitful: mm, s

ig = ImageGeom( dims=(640,480,1000), deltas=(1mm,1mm,(1//12)s) )

# ### Methods

# An ImageGeom object has quite a few methods;
# `axes` and `axis` are especially useful:

ig = ImageGeom( dims=(7,8), deltas=(3,2), offsets=(0,0.5) )
axis(ig, 2)

# For an axis of length `n` with spacing `Δ` (possibly with units)
# and (always unitless but possibly non-integer) `offset` the axis
# is a subtype of `AbstractRange` of the form
# `( (0:n-1) .- ((n - 1)/2 + offset) ) * Δ`

# These axes are useful for plotting:
ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm), offsets=(0.5,0.5) )

#
showgrid = (ig) -> begin # x,y grid locations of pixel centers
    x = axis(ig, 1)
    y = axis(ig, 2)
    (xg, yg) = grids(ig)
    scatter(xg, yg, label="", xlabel="x", ylabel="y",
        xlims = maximum(abs, x) * 1.2 .* (-1,1),
        xticks = [x[1], zero(eltype(x)), x[end]],
        ylims = maximum(abs, y) * 1.2 .* (-1,1),
        yticks = [y[1], zero(eltype(y)), y[end]],
        aspect_ratio=1, title="offsets $(ig.offsets)")
end
showgrid(ig)

# Unit labels on the axes due to `UnitfulRecipes.jl`
prompt();

# ### Offsets (unitless translation of grid)

# The default `offsets` are zeros,
# corresponding to symmetric sampling around origin:

ig = ImageGeom( dims=(12,10), deltas=(1mm,1mm) )
p = showgrid(ig)

#
prompt();

# That default for `offsets` is natural for tomography
# when considering finite pixel size:

square = (x,y,Δ) -> plot!(p, label="", color=:black,
    x .+ Δ[1] * ([0,1,1,0,0] .- 0.5),
    y .+ Δ[2] * ([0,0,1,1,0] .- 0.5),
)
square2 = (x,y) -> square(x, y, ig.deltas)
square2.(grids(ig)...)
plot!(p)

#
prompt();

# In that default geometry, the center `(0,0)` of the image
# is at a corner of the middle 4 pixels (for even image sizes).
# That default is typical for tomographic imaging (e.g., CT, PET, SPECT).
# One must be careful when using operations like `imrotate` or `fft`.

# ### AxisArrays

# There is a natural connection between `ImageGeom` and `AxisArrays`.
# Note the automatic labeling of units (when relevant) on all axes by
# [MIRTjim.jim](https://github.com/JeffFessler/MIRTjim.jl).

using AxisArrays
using Unitful: mm
ig = ImageGeom( dims=(60,48), deltas=(1.5mm,1mm) )
za = AxisArray( ellipse(ig) * 10/mm ; x=axis(ig,1), y=axis(ig,2) )
jim(za, "AxisArray example")

#
prompt();

# ### Resizing

# Often we have a target grid in mind but want coarser sampling for debugging.
# The `downsample` method is useful for this.

ig = ImageGeom( dims = (512,512), deltas = (500mm,500mm) ./ 512 )
ig_down = downsample(ig, 4)

# Other times we want to avoid an "inverse crime" by using finer sampling
# to simulate data; use `oversample` for this.

ig_over = oversample(ig, (2,2))

# Just for debugging:
using Pkg
Pkg.status("ImageGeoms")


text:
# infrastructure.jl

using DataFrames
#using SQLite
#using Dates
#using SQLite

using Serialization

include("../api/api.jl")
include("./db.jl")
include("./doc.jl")

const VAT_ACCOUNT = 4000

process(entries::Array{JournalEntry, 1}; path_journal="./test_journal.txt", path_ledger="./test_ledger.txt") = begin
    # archive journal entries
    add_to_file(path_journal, entries)

    # store records
    for entry in entries
        debit = create_record_from(entry)
        credit = create_record_to(entry)
        vat = create_record_vat(entry, VAT_ACCOUNT)

        records = entry.vat != 0 ? [debit, credit, vat] : [debit, credit]

        add_to_file(path_ledger, records)
    end

end # process


text:
using Makie

s=Scene()
xy1, xz1, xy2, xz2 = -1., -1., 1., 1.
lines!(s,[0.,0.,0.,0.,0.],[xy1,xy1,xy2,xy2,xy1],[xz1,xz2,xz2,xz1,xz1])
b=vec(Point3f0.([0.,0.,0.,0.],[xy1,xy1,xy2,xy2],[xz1,xz2,xz2,xz1]))
u=vec(Point3f0.([0.,0.,0.,0.],[0,xy2-xy1,0.,xy1-xy2],[xz2-xz1,0.,xz1-xz2,0]))
arrows!(s,b,u)

yz1, yx1, yz2, yx2 = -1., -1., 1., 1.
lines!(s,[yx1,yx1,yx2,yx2,yx1],[0.,0.,0.,0.,0.],[yz1,yz2,yz2,yz1,yz1])
b=vec(Point3f0.([yx1,yx2,yx2,yx1],[0.,0.,0.,0.],[yz1,yz1,yz2,yz2]))
u=vec(Point3f0.([yx2-yx1,0,yx1-yx2,0],[0.,0.,0.,0.],[0,yz2-yz1,0,yz1-yz2]))
arrows!(s,b,u)

zx1, zy1, zx2, zy2 = -1., -1., 1., 1.
lines!(s,[zx1,zx1,zx2,zx2,zx1],[zy1,zy2,zy2,zy1,zy1],[0.,0.,0.,0.,0.])
b=vec(Point3f0.([zx1,zx2,zx2,zx1],[zy1,zy1,zy2,zy2],[0.,0.,0.,0.]))
u=vec(Point3f0.([zx2-zx1,0,zx1-zx2,0],[0,zy2-zy1,0,zy1-zy2],[0.,0.,0.,0.]))
arrows!(s,b,u)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
```
mutation(loglikelihood::Function, parameters::ParameterVector{U},
         data::Matrix{S}, p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
         blocks_free::Vector{Vector{Int}}, blocks_all::Vector{Vector{Int}},
         ϕ_n::S, ϕ_n1::S; c::S = 1., α::S = 1., n_mh_steps::Int = 1,
         old_data::T = T(undef, size(data, 1), 0)) where {S<:AbstractFloat,
                                                          T<:AbstractMatrix, U<:Number}
```

Execute one proposed move of the Metropolis-Hastings algorithm for a given parameter

### Arguments:
- `loglikelihood::Function`: Likelihood function of model being estimated.
- `parameters::ParameterVector{U}`: Model parameter vector, which stores parameter
    values, prior dists, and bounds
- `data::Matrix{Float64}`: Matrix of data
- `p::Vector{Float64}`: Initial particle value
- `d::Distribution`: A distribution with μ = the weighted mean, and Σ = the weighted
    variance/covariance matrix
- `blocks_free::Vector{Vector{Int64}}`: A vector of index blocks, where the indices
    in each block corresponds to the ordering of free parameters only (e.g. all the
    indices will be ∈ 1:n_free_parameters)
- `blocks_all::Vector{Vector{Int64}}`: A vector of index blocks, where the indices in
    each block corresponds to the ordering of all parameters (e.g. all the indices will
    be in ∈ 1:n_para, free and fixed)
- `ϕ_n::Float64`: The current tempering factor
- `ϕ_n1::Float64`: The previous tempering factor

### Keyword Arguments:
- `c::Float64`: The scaling parameter for the proposed covariance matrix
- `α::Float64`: The mixing proportion
- `n_mh_steps::Int`: Number of Metropolis Hastings steps to attempt
- `old_data::Matrix{Float64}`: The matrix of old data to be used in calculating the
    old_loglh, old_logpost in time tempering

### Outputs:
- `p::Vector{Float64}`: An updated particle containing updated parameter values,
    log-likelihood, prior, and acceptance indicator.

"""
function mutation(loglikelihood::Function, parameters::ParameterVector{U},
                  data::Matrix{S}, p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                  blocks_free::Vector{Vector{Int}}, blocks_all::Vector{Vector{Int}},
                  ϕ_n::S, ϕ_n1::S; c::S = 1., α::S = 1., n_mh_steps::Int = 1,
                  old_data::T = T(undef, size(data, 1), 0)) where {S<:AbstractFloat,
                                                                   T<:AbstractMatrix, U<:Number}

    n_free_para = length([!θ.fixed for θ in parameters])
    step_prob   = rand() # Draw initial step probability

    N         = length(p)
    para      = p[1:ind_para_end(N)]
    like      = p[ind_loglh(N)]
    logprior  = p[ind_logprior(N)]
    like_prev = p[ind_old_loglh(N)] # Likelihood evaluated at the old data (for time tempering)
    accept    = 0.0

    d = MvNormal(d_μ, d_Σ)

    for step in 1:n_mh_steps
        for (block_f, block_a) in zip(blocks_free, blocks_all)

            # Index out parameters corresponding to given random block, create distribution
            # centered at weighted mean, with Σ corresponding to the same random block
            para_subset = para[block_a]
            d_subset    = MvNormal(d.μ[block_f], d.Σ.mat[block_f, block_f])
            para_draw   = mvnormal_mixture_draw(para_subset, d_subset; c = c, α = α)

            q0, q1 = compute_proposal_densities(para_draw, para_subset,
                                                d_subset, c = c, α = α)

            para_new          = deepcopy(para)
            para_new[block_a] = para_draw

            like_init, prior_init = like, logprior
            prior_new = like_new = like_old_data = -Inf
            try
                update!(parameters, para_new)
                para_new  = [θ.value for θ in parameters]
                prior_new = prior(parameters)
                like_new  = loglikelihood(parameters, data)

                if like_new == -Inf
                    prior_new = like_old_data = -Inf
                end

                like_old_data = isempty(old_data) ? 0. : loglikelihood(parameters, old_data)

            catch err
                if isa(err, ParamBoundsError) || isa(err, LinearAlgebra.LAPACKException) ||
                   isa(err, PosDefException)  || isa(err, SingularException)
                    prior_new = like_new = like_old_data = -Inf
                else
                    throw(err)
                end
            end

            η = exp(ϕ_n * (like_new - like_init) + (1 - ϕ_n) * (like_old_data - like_prev) +
                    (prior_new - prior_init) + (q0 - q1))

            if step_prob < η
                para      = para_new
                like      = like_new
                logprior  = prior_new
                like_prev = like_old_data
                accept   += length(block_a)
            end
            step_prob = rand() # Draw again for next step
        end
    end
    update_mutation!(p, para, like, logprior, like_prev, accept / n_free_para)
    return p
end


text:
module MyExample2

using ForwardDiff
# Write your package code here.

greet() = print("Hello World!")
include("demo.jl")
df(x,y) = ForwardDiff.derivative(x->f(x, y),x)

end


text:
include("C:/Users/stefano/Documents/Zhulab/julia_code/model_fit.jl")
include("C:/Users/stefano/Documents/Zhulab/julia_code/models.jl")
include("C:/Users/stefano/Documents/Zhulab/julia_code/helper_functions.jl")
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using BlackBoxOptim
using BenchmarkTools
using RecursiveArrayTools
# --------------------------create dummy data-----------------------------------
f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens,cons) # enclose constants

# u0 = [0.0;0.0]
u₀ = [0.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
p = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
dens = [100., 10., 1.]
cons = []

# Setup the ODE problem, then solve
prob = ODEProblem(f,u₀,tspan,p)
sol = solve(prob,Tsit5())
# sol = solve(prob,Tsit5())
plot(sol)
dataset1 = generate_data(sol)

model_opt2 = (1.0e-6,10.0)
# pinit = [2.0e-4,0.6]
pinit = [-4.0,0.0]
t = collect(range(0,stop=10,length=101))
rates = optimization(bimolecular!,dataset1,dens',[],model_opt2,t)

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

# u0 = [0.0;0.0]
u₀ = [0.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
p = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
dens_new = [100., 50., 1.]
cons = []

f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens_new,cons) # enclose constants
# Setup the ODE problem, then solve
prob = ODEProblem(f,u₀,tspan,p)
sol = solve(prob,Tsit5())
# sol = solve(prob,Tsit5())
plot!(sol)
dataset3 = generate_data(sol)
dataset_comb = vcat(dataset1,dataset3)

dens_comb = reduce(vcat, [dens', dens_new'])
rates = optimization(bimolecular!,dataset_comb,dens_comb,cons,model_opt2,t)

#----------------------------------------------------------------------------
# try TCR-pMHC only data from Muaz

using CSV
tcr_data = CSV.read("C:/Users/stefano/Documents/Zhulab/julia_code/tcr_only_data.csv")

t, Pa, n_se_tcr = clean_data(tcr_data)
n = Pa2n(Pa)

scatter(t,Pa)
scatter!(t,n)
scatter!(n_se_tcr[:,1],n_se_tcr[:,2],yerror=n_se_tcr[:,3],markercolor=:black)

densities_tcr = [25.0, 38.0, 1.0]
tspan = (0,maximum(t)+1)
u₀ = [0.0]
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_tcr,[]) # enclose constants

pinit = [-4.5,-0.7]

rates_tcr = optimization(bimolecular!,n',densities_tcr',[],model_opt2,t)

prob = ODEProblem(f,u₀,tspan,10 .^ rates_tcr)
# prob = remake(prob, p=10 .^ rates)
sol = solve(prob,Rosenbrock23())
plot!(sol,linecolor=:black)
summ = reduce(vcat,sum(sol,dims=1))
h = plot!(sol.t,summ,linecolor=:black)

# prob = ODEProblem(f,u₀,tspan,[1e-4,0.5])

#------------------------------------------------------------------------------
# try TCR-pMHC-CD4 fit with 2-path, 2-step

tcr_cd4_1 = CSV.read("C:/Users/stefano/Documents/Zhulab/julia_code/tcr_cd4_1.csv")
tcr_cd4_2 = CSV.read("C:/Users/stefano/Documents/Zhulab/julia_code/tcr_cd4_2.csv")
tcr_cd4_3 = CSV.read("C:/Users/stefano/Documents/Zhulab/julia_code/tcr_cd4_3.csv")
t_tc1, Pa_tc1, n_se_1 = clean_data(tcr_cd4_1)
t_tc2, Pa_tc2,  n_se_2  = clean_data(tcr_cd4_2)
t_tc3, Pa_tc3,  n_se_3 = clean_data(tcr_cd4_3)

n_tc1 = Pa2n(Pa_tc1)
n_tc2 = Pa2n(Pa_tc2)
n_tc3 = Pa2n(Pa_tc3)

# scatter(t,Pa)
scatter(t_tc1,n_tc1)
scatter!(t_tc2,n_tc2)
scatter!(t_tc3,n_tc3)

densities_tc1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
tspan = (0,maximum(t)+1)
u₀ = [0.0, 0.0]

cons_tcr = 10 .^ rates_tcr
# f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc1,cons_tcr)
# rates_tc1 = optimization(trimolecular!,n_tc1',densities_tc1',cons_tcr,model_opt2,t_tc1)

u₀ = [0.0, 0.0, 0.0]
# CD4 rates from Muaz, need to get full CD4 curve
kon2 = 7.01979340915551e-06 / 6
koff2 = 0.4000
cons_tc = vcat(cons_tcr, kon2, koff2)

# cons_tcr_from_mat = [1.4e-4, 0.5]
# cons_tc = vcat(cons_tcr_from_mat, kon2, koff2)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc1,cons_tc) # enclose constants
# p = [9.32e-7,1.81e-1,1.7e-2,5.23e-5]
p = [0.0, 0.0, 0.0, 0.0]
prob = ODEProblem(f,u₀,tspan,p)
# prob = remake(prob, p=10 .^ rates)
sol = solve(prob,Rosenbrock23())
summ = reduce(vcat,sum(sol,dims=1))
plot!(sol.t,summ)

pinit = [1e-2,0.5,1e-2,0.5]
rates_tc1_2paths = optimization(two_path!,n_tc1',densities_tc1',cons_tc,model_opt2,t_tc1)

densities_tc2 = [8.0, 38.0, 17.0, 1.0]
rates_tc2_2paths = optimization(two_path!,n_tc2',densities_tc2',cons_tc,model_opt2,t_tc2)

densities_tc3 = [3.0, 38.0, 22.0, 1.0]
rates_tc3_2paths = optimization(two_path!,n_tc3',densities_tc3',cons_tc,model_opt2,t_tc3)

scatter(n_se_1[:,1],n_se_1[:,2],yerror=n_se_1[:,3],markercolor=:red)
scatter!(n_se_2[:,1],n_se_2[:,2],yerror=n_se_2[:,3],markercolor=:blue)
scatter!(n_se_3[:,1],n_se_3[:,2],yerror=n_se_3[:,3],markercolor=:cyan)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc1,cons_tc)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc1_2paths)
sol = solve(prob,Rosenbrock23())
summ1 = reduce(vcat,sum(sol,dims=1))
time = sol.t
plot!(time,summ1,linecolor=:red)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc2,cons_tc)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc2_2paths)
sol = solve(prob,Rosenbrock23())
summ = reduce(vcat,sum(sol,dims=1))
time = sol.t
plot!(sol.t,summ,linecolor=:blue)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc3,cons_tc)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc3_2paths)
sol = solve(prob,Rosenbrock23())
summ = reduce(vcat,sum(sol,dims=1))
plot!(sol.t,summ,linecolor=:cyan)

function plot_sum(rates)
    prob = ODEProblem(f,u₀,tspan,rates)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Rosenbrock23())
    summ = reduce(vcat,sum(sol,dims=1))
    plot!(sol.t,summ)
end

# dprob = DiscreteProblem(u₀, tspan, 10 .^ rates)
#
# # now we create a JumpProblem, and specify Gillespie's Direct Method as the solver:\n",
# jprob = JumpProblem(dprob, Direct(), f)
#
# # now let's solve and plot the jump process:
# sol = solve(jprob, SSAStepper())
# plot(sol, fmt=:svg, ylims = (0, 2))
#
# prob = DiscreteProblem(bimol, u₀, tspan, 10 .^ rates)
# jump_prob = JumpProblem(bimol, prob, Direct())
# sol = solve(jump_prob, SSAStepper())
#
# plot(sol)
#
# function bimol(du,u,p,t)
#     # unpack rates and constants
#     nᵣ = u[1]
#     k₁,k₋₁  = p
#     mᵣ,mₗ,A = [100,10,1]
#     # model
#     du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ
#
# end


text:
"""
```
sendto(p::Int; args...)
```
Function to send data from master process to particular worker, p. Code from ChrisRackauckas, avavailable at: https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, Core.eval(Main, Expr(:(=), nm, val)))
    end
end

"""
```
sendto(ps::AbstractVector{Int}; args...)
```
Function to send data from master process to list of workers. Code from ChrisRackauckas, available at: https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(ps::AbstractVector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end

"""
```
abbrev_symbol(s::Symbol, n::Int=4)
```

Abbreviate the symbol `s` to an `n`-character string.
"""
function abbrev_symbol(s::Symbol, n::Int=4)
    str = string(s)
    if length(str) ≤ n
        return str
    else
        return str[1:n]
    end
end

function sorted_list_insert!(v::Vector{T}, x::T) where T
    insert_index = 1
    for val in v
        if x<val
            break
        else
            insert_index += 1
        end
    end
    insert!(v,insert_index,x)
end

"""
```
quarter_range(t0::Date, t1::Date)
```

Returns a vector of `Dates`, consisting of the last days of each quarter between
`t0` and `t1`, inclusive.
"""
function quarter_range(t0::Date, t1::Date)
    dr = t0:Day(1):t1
    return Base.filter(d -> Dates.lastdayofquarter(d) == d, dr)
end

"""
```
detexify(s::String)

detexify(s::Symbol)
```

Remove Unicode characters from the string `s`, replacing them with ASCII
equivalents. For example, `detexify(\"π\")` returns `\"pi\"`.
"""
function detexify(s::String)
    s = replace(s, "α" => "alpha")
    s = replace(s, "β" => "beta")
    s = replace(s, "γ" => "gamma")
    s = replace(s, "δ" => "delta")
    s = replace(s, "ϵ" => "epsilon")
    s = replace(s, "ε" => "epsilon")
    s = replace(s, "ζ" => "zeta")
    s = replace(s, "η" => "eta")
    s = replace(s, "θ" => "theta")
    s = replace(s, "ι" => "iota")
    s = replace(s, "κ" => "kappa")
    s = replace(s, "λ" => "lambda")
    s = replace(s, "μ" => "mu")
    s = replace(s, "ν" => "nu")
    s = replace(s, "ξ" => "xi")
    s = replace(s, "π" => "pi")
    s = replace(s, "ρ" => "rho")
    s = replace(s, "ϱ" => "rho")
    s = replace(s, "σ" => "sigma")
    s = replace(s, "ς" => "sigma")
    s = replace(s, "τ" => "tau")
    s = replace(s, "υ" => "upsilon")
    s = replace(s, "ϕ" => "phi")
    s = replace(s, "φ" => "phi")
    s = replace(s, "χ" => "chi")
    s = replace(s, "ψ" => "psi")
    s = replace(s, "ω" => "omega")

    s = replace(s, "Α" => "Alpha")
    s = replace(s, "Β" => "Beta")
    s = replace(s, "Γ" => "Gamma")
    s = replace(s, "Δ" => "Delta")
    s = replace(s, "Ε" => "Epsilon")
    s = replace(s, "Ζ" => "Zeta")
    s = replace(s, "Η" => "Eta")
    s = replace(s, "Θ" => "Theta")
    s = replace(s, "Ι" => "Iota")
    s = replace(s, "Κ" => "Kappa")
    s = replace(s, "Λ" => "Lambda")
    s = replace(s, "Μ" => "Mu")
    s = replace(s, "Ν" => "Nu")
    s = replace(s, "Ξ" => "Xi")
    s = replace(s, "Π" => "Pi")
    s = replace(s, "Ρ" => "Rho")
    s = replace(s, "Σ" => "Sigma")
    s = replace(s, "Τ" => "Tau")
    s = replace(s, "Υ" => "Upsilon")
    s = replace(s, "Φ" => "Phi")
    s = replace(s, "Χ" => "Chi")
    s = replace(s, "Ψ" => "Psi")
    s = replace(s, "Ω" => "Omega")

    s = replace(s, "′" => "'")

    return s
end

function detexify(s::Symbol)
    Symbol(detexify(string(s)))
end

dispfns = [:print, :println]
for disp in dispfns
    @eval begin
        function Base.$disp(verbose::Symbol, min::Symbol, xs...)
            if VERBOSITY[verbose] >= VERBOSITY[min]
                $disp(xs...)
            end
        end
    end
end

function info_print(verbose::Symbol, min::Symbol, x::String)
    if VERBOSITY[verbose] >= VERBOSITY[min]
        @info(x)
    end
end

function warn_print(verbose::Symbol, min::Symbol, x::String)
    if VERBOSITY[verbose] >= VERBOSITY[min]
        @warn(x)
    end
end

## Testing functions

function test_matrix_eq2(expect::AbstractArray,
                         actual::AbstractArray,
                         expstr::String,
                         actstr::String,
                         ϵ_abs::Float64 = 1e-6,
                         ϵ_rel::Float64 = 1e-2) where {T<:AbstractFloat}
    if length(expect) ≠ length(actual)
        error("") #="lengths of ", expstr, " and ", actstr, " do not match: ",
              "\n  ", expstr, " (length $(length(expect))) = ", expect,
              "\n  ", actstr, " (length $(length(actual))) = ", actual)=#
    end

    # Absolute difference filter
    abs_diff   = abs.(actual .- expect) .> ϵ_abs
    n_abs_diff = sum(skipmissing(abs_diff))

    # Relative difference filter
    rel_diff   = 100*abs.((actual .- expect) ./ expect) .> ϵ_rel
    n_rel_diff = sum(skipmissing(rel_diff))

    # Element is only problematic if it fails *both* tests.
    mixed_diff   = abs_diff .& rel_diff
    n_mixed_diff = sum(skipmissing(mixed_diff))

    if n_mixed_diff ≠ 0
        sdiff = string("|a - b| <= ", ϵ_abs,
                   " or |a - b|/|b| <= ", ϵ_rel, "%,",
                   " ∀ a ∈ ", actstr, ",",
                   " ∀ b ∈ ",expstr)
        @warn "assertion failed:\n",
             "    ", sdiff,
             "\n$(n_abs_diff) entries fail absolute equality filter",
             "\n$(n_rel_diff) entries fail relative equality filter",
             "\n$(n_mixed_diff) entries fail both equality filters\n"
        return false
    end
    return true
end

"""
    @test_matrix_approx_eq(a, b)

Test two matrices of floating point numbers `a` and `b` for approximate equality.
"""
macro test_matrix_approx_eq(a,b)
    :(test_matrix_eq2($(esc(a)),$(esc(b)),$(string(a)),$(string(b))))
end

"""
    @test_matrix_approx_eq_eps(a, b, ϵ_abs, ϵ_rel)

Test two matrices of floating point numbers `a` and `b` for equality taking in account a
margin of absolute tolerance given by `ϵ_abs` and a margin of relative tolerance given by
`ϵ_rel` (in comparison to `b`).
"""
macro test_matrix_approx_eq_eps(a,b,c,d)
    :(test_matrix_eq2($(esc(a)),$(esc(b)),$(string(a)),$(string(b)),$(esc(c)),$(esc(d))))
end

"""
Sparse identity matrix - since deprecated in 0.7
"""
function speye(n::Integer)
    return SparseMatrixCSC{Float64}(I, n, n)
end

"""
Sparse identity matrix - since deprecated in 0.7
"""
function speye(T::Type, n::Integer)
    return SparseMatrixCSC{T}(I, n, n)
end

"""
    <(a::Complex, b::Complex)

Compare real values of complex numbers.
"""
function <(a::Complex, b::Complex)
    return a.re < b.re
end

"""
    <(a::Real, b::Complex)

Compare real values of complex numbers.
"""
function <(a::Real, b::Complex)
    return a < b.re
end

"""
    <(a::Complex, b::Real)

Compare real values of complex numbers.
"""
function <(a::Complex, b::Real)
    return a.re < b
end

function min(a::Complex, b::Real)
    return min(a.re, b)
end

function min(a::Complex, b::Complex)
    return min(a.re, b.re)
end

function min(a::Real, b::Complex)
    return min(a, b.re)
end

function max(a::Complex, b::Real)
    return max(a.re, b)
end

function max(a::Complex, b::Complex)
    return max(a.re, b.re)
end

function max(a::Real, b::Complex)
    return max(a, b.re)
end


text:
"""
    DiskTrav(A::AbstractMatrix)

    takes coefficients as provided by the Zernike polynomial layout of FastTransforms.jl and
    makes them accessible sorted such that in each block the m=0 entries are always in first place,
    followed by alternating sin and cos terms of increasing |m|.
"""
struct DiskTrav{T, AA<:AbstractMatrix{T}} <: AbstractBlockVector{T}
    matrix::AA
    function DiskTrav{T, AA}(matrix::AA) where {T,AA<:AbstractMatrix{T}}
        m,n = size(matrix)
        isodd(n) && m == n ÷ 4 + 1 || throw(ArgumentError("size must match"))
        new{T,AA}(matrix)
    end
end

DiskTrav{T}(matrix::AbstractMatrix{T}) where T = DiskTrav{T,typeof(matrix)}(matrix)
DiskTrav(matrix::AbstractMatrix{T}) where T = DiskTrav{T}(matrix)

function DiskTrav(v::AbstractVector{T}) where T
    N = blocksize(v,1)
    m = N ÷ 2 + 1
    n = 4(m-1) + 1
    mat = zeros(T, m, n)
    for K in blockaxes(v,1)
        K̃ = Int(K)
        w = v[K]
        if isodd(K̃)
            mat[K̃÷2 + 1,1] = w[1]
            for j = 2:2:K̃-1
                mat[K̃÷2-j÷2+1,2(j-1)+2] = w[j]
                mat[K̃÷2-j÷2+1,2(j-1)+3] = w[j+1]
            end
        else
            for j = 1:2:K̃
                mat[K̃÷2-j÷2,2(j-1)+2] = w[j]
                mat[K̃÷2-j÷2,2(j-1)+3] = w[j+1]
            end
        end
    end
    DiskTrav(mat)
end

axes(A::DiskTrav) = (blockedrange(oneto(div(size(A.matrix,2),2,RoundUp))),)

function getindex(A::DiskTrav, K::Block{1})
    k = Int(K)
    k == 1 && return A.matrix[1:1]
    k == 2 && return A.matrix[1,2:3]
    st = stride(A.matrix,2)
    if isodd(k)
        # nonnegative terms
        p = A.matrix[range(k÷2+1, step=4*st-1, length=k÷2+1)]
        # negative terms
        n = A.matrix[range(k÷2+3*st, step=4*st-1, length=k÷2)]
        interlace(p,n)
    else
        # negative terms
        n = A.matrix[range(st+k÷2, step=4*st-1, length=k÷2)]
        # positive terms
        p = A.matrix[range(2st+k÷2, step=4*st-1, length=k÷2)]
        interlace(n,p)
    end
end

getindex(A::DiskTrav, k::Int) = A[findblockindex(axes(A,1), k)]

ClassicalOrthogonalPolynomials.checkpoints(d::UnitDisk{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]

"""
    ZernikeWeight(a, b)

is a quasi-vector representing `r^(2a) * (1-r^2)^b`
"""
struct ZernikeWeight{T} <: Weight{T}
    a::T
    b::T
end

"""
    ZernikeWeight(b)

is a quasi-vector representing `(1-r^2)^b`
"""

ZernikeWeight(b) = ZernikeWeight(zero(b), b)
ZernikeWeight{T}(b) where T = ZernikeWeight{T}(zero(T), b)
ZernikeWeight{T}() where T = ZernikeWeight{T}(zero(T))
ZernikeWeight() = ZernikeWeight{Float64}()

copy(w::ZernikeWeight) = w

axes(::ZernikeWeight{T}) where T = (Inclusion(UnitDisk{T}()),)

==(w::ZernikeWeight, v::ZernikeWeight) = w.a == v.a && w.b == v.b

function getindex(w::ZernikeWeight, xy::StaticVector{2})
    r = norm(xy)
    r^(2w.a) * (1-r^2)^w.b
end

"""
    Zernike(a, b)

is a quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`
"""
struct Zernike{T} <: BivariateOrthogonalPolynomial{T}
    a::T
    b::T
    Zernike{T}(a::T, b::T) where T = new{T}(a, b)
end
Zernike{T}(a, b) where T = Zernike{T}(convert(T,a), convert(T,b))
Zernike(a::T, b::V) where {T,V} = Zernike{float(promote_type(T,V))}(a, b)
Zernike{T}(b) where T = Zernike{T}(zero(b), b)
Zernike{T}() where T = Zernike{T}(zero(T))

"""
    Zernike(b)

is a quasi-matrix orthogonal `(1-r^2)^b`
"""
Zernike(b) = Zernike(zero(b), b)
Zernike() = Zernike{Float64}()

axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))

==(w::Zernike, v::Zernike) = w.a == v.a && w.b == v.b

copy(A::Zernike) = A

orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)

zerniker(ℓ, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
zerniker(ℓ, m, r) = zerniker(ℓ, m, zero(r), r)

function zernikez(ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
    r,θ = rθ.r,rθ.θ
    m = abs(ms)
    zerniker(ℓ, m, a, b, r) * (signbit(ms) ? sin(m*θ) : cos(m*θ))
end

zernikez(ℓ, ms, a, b, xy::StaticVector{2}) = zernikez(ℓ, ms, a, b, RadialCoordinate(xy))
zernikez(ℓ, ms, b, xy::StaticVector{2}) = zernikez(ℓ, ms, zero(b), b, xy)
zernikez(ℓ, ms, xy::StaticVector{2,T}) where T = zernikez(ℓ, ms, zero(T), xy)

function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
    ℓ = Int(block(B))-1
    k = blockindex(B)
    m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
    zernikez(ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
end

getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate(xy), B]
getindex(Z::Zernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)]
getindex(Z::Zernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])])

###
# Transforms
###

const FiniteZernike{T} = SubQuasiArray{T,2,Zernike{T},<:Tuple{<:Inclusion,<:BlockSlice{BlockRange1{OneTo{Int}}}}}

function grid(S::FiniteZernike{T}) where T
    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
    M = 4N-3

    r = sinpi.((N .-(0:N-1) .- one(T)/2) ./ (2N))

    # The angular grid:
    θ = (0:M-1)*convert(T,2)/M
    RadialCoordinate.(r, π*θ')
end

_angle(rθ::RadialCoordinate) = rθ.θ

function plotgrid(S::FiniteZernike{T}) where T
    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
    g = grid(parent(S)[:,Block.(OneTo(2N))]) # double sampling
    θ = [map(_angle,g[1,:]); 0]
    [permutedims(RadialCoordinate.(1,θ)); g g[:,1]; permutedims(RadialCoordinate.(0,θ))]
end

function plotgrid(S::SubQuasiArray{<:Any,2,<:Zernike})
    kr,jr = parentindices(S)
    Z = parent(S)
    plotgrid(Z[kr,Block.(OneTo(Int(findblock(axes(Z,2),maximum(jr)))))])
end

struct ZernikeTransform{T} <: Plan{T}
    N::Int
    disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
    analysis::FastTransforms.FTPlan{T,2,FastTransforms.DISKANALYSIS}
end

struct ZernikeITransform{T} <: Plan{T}
    N::Int
    disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
    synthesis::FastTransforms.FTPlan{T,2,FastTransforms.DISKSYNTHESIS}
end

function ZernikeTransform{T}(N::Int, a::Number, b::Number) where T<:Real
    Ñ = N ÷ 2 + 1
    ZernikeTransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_analysis(T, Ñ, 4Ñ-3))
end
function ZernikeITransform{T}(N::Int, a::Number, b::Number) where T<:Real
    Ñ = N ÷ 2 + 1
    ZernikeITransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_synthesis(T, Ñ, 4Ñ-3))
end

*(P::ZernikeTransform{T}, f::AbstractArray) where T = P * convert(Matrix{T}, f)
*(P::ZernikeTransform{T}, f::Matrix{T}) where T = DiskTrav(P.disk2cxf \ (P.analysis * f))[Block.(1:P.N)]
*(P::ZernikeITransform, f::AbstractVector) = P.synthesis * (P.disk2cxf * DiskTrav(f).matrix)

factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))

# gives the entries for the Laplacian times (1-r^2) * Zernike(1)
struct WeightedZernikeLaplacianDiag{T} <: AbstractBlockVector{T} end

axes(::WeightedZernikeLaplacianDiag) = (blockedrange(oneto(∞)),)
copy(R::WeightedZernikeLaplacianDiag) = R

MemoryLayout(::Type{<:WeightedZernikeLaplacianDiag}) = LazyLayout()
Base.BroadcastStyle(::Type{<:Diagonal{<:Any,<:WeightedZernikeLaplacianDiag}}) = LazyArrayStyle{2}()

function Base.view(W::WeightedZernikeLaplacianDiag{T}, K::Block{1}) where T
    K̃ = Int(K)
    d = K̃÷2 + 1
    if isodd(K̃)
        v = (d:K̃) .* (d:(-1):1)
        convert(AbstractVector{T}, -4*interlace(v, v[2:end]))
    else
        v = (d:K̃) .* (d-1:(-1):1)
        convert(AbstractVector{T}, -4*interlace(v, v))
    end
end

getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,1),k)]

@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike})
    @assert WZ.P.a == 0 && WZ.P.b == 1
    WZ.P * Diagonal(WeightedZernikeLaplacianDiag{eltype(eltype(WZ))}())
end

###
# Fractional Laplacian
###

function *(L::AbsLaplacianPower, WZ::Weighted{<:Any,<:Zernike{<:Any}})
    @assert axes(L,1) == axes(WZ,1) && WZ.P.a == 0 && WZ.P.b == L.α
    WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(L.α)}(L.α))
end

# gives the entries for the (negative!) fractional Laplacian (-Δ)^(α) times (1-r^2)^α * Zernike(α)
struct WeightedZernikeFractionalLaplacianDiag{T} <: AbstractBlockVector{T}
    α::T
end

axes(::WeightedZernikeFractionalLaplacianDiag) = (blockedrange(oneto(∞)),)
copy(R::WeightedZernikeFractionalLaplacianDiag) = R

MemoryLayout(::Type{<:WeightedZernikeFractionalLaplacianDiag}) = LazyLayout()
Base.BroadcastStyle(::Type{<:Diagonal{<:Any,<:WeightedZernikeFractionalLaplacianDiag}}) = LazyArrayStyle{2}()

getindex(W::WeightedZernikeFractionalLaplacianDiag, k::Integer) = W[findblockindex(axes(W,1),k)]

function Base.view(W::WeightedZernikeFractionalLaplacianDiag{T}, K::Block{1}) where T
    l = Int(K)
    if isodd(l)
        m = Vcat(0,interlace(Array(2:2:l),Array(2:2:l)))
    else #if iseven(l)
        m = Vcat(interlace(Array(1:2:l),Array(1:2:l)))
    end
    return convert(AbstractVector{T}, 2^(2*W.α)*fractionalcfs2d.(l-1,m,W.α))
end

# generic d-dimensional ball fractional coefficients without the 2^(2*β) factor. m is assumed to be entered as abs(m)
function fractionalcfs(l::Integer, m::Integer, α::T, d::Integer) where T
    n = (l-m)÷2
    return exp(loggamma(α+n+1)+loggamma((2*α+2*n+d+2*m)/2)-loggamma(one(T)+n)-loggamma((2*one(T)*m+2*n+d)/2))
end
# 2 dimensional special case, again without the 2^(2*β) factor
fractionalcfs2d(l::Integer, m::Integer, β) = fractionalcfs(l,m,β,2)

"""
    ModalInterlace
"""
struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
    ops
    MN::MMNN
    bandwidths::NTuple{2,Int}
end

ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T =
    ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)

# act like lazy array
MemoryLayout(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyBandedBlockBandedLayout()
Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()

axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))

blockbandwidths(R::ModalInterlace) = R.bandwidths
subblockbandwidths(::ModalInterlace) = (0,0)

function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
    K,J = KJ.n
    dat = Matrix{T}(undef,1,J)
    l,u = blockbandwidths(R)
    if iseven(J-K) && -l ≤ J - K ≤ u
        sh = (J-K)÷2
        if isodd(K)
            k = K÷2+1
            dat[1,1] = R.ops[1][k,k+sh]
        end
        for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
            k = K÷2-m÷2+isodd(K)
            dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
        end
    else
        fill!(dat, zero(T))
    end
    _BandedMatrix(dat, K, 0, 0)
end

getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]

function getindex(R::ModalInterlace{T}, KR::BlockOneTo, JR::BlockOneTo) where T
    M,N = Int(last(KR)), Int(last(JR))
    ModalInterlace{T}([R.ops[m][1:(M-m+2)÷2,1:(N-m+2)÷2] for m=1:min(N,M)], (M,N), R.bandwidths)
end

function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
    TV = promote_type(T,V)
    A.a == B.a && A.b == B.b && return Eye{TV}(∞)
    @assert A.a == 0 && A.b == 1
    @assert B.a == 0 && B.b == 0
    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(1,0:∞)) .\ Normalized.(Jacobi{TV}.(0,0:∞))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (0,2))
end

function \(A::Zernike{T}, B::Weighted{V,Zernike{V}}) where {T,V}
    TV = promote_type(T,V)
    A.a == B.P.a == A.b == B.P.b == 0 && return Eye{TV}(∞)
    if A.a == A.b == 0
        @assert B.P.a == 0 && B.P.b == 1
        ModalInterlace{TV}((Normalized.(Jacobi{TV}.(0, 0:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(1, 0:∞)))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (2,0))
    else
        Z = Zernike{TV}()
        (A \ Z) * (Z \ B)
    end
end


text:
"""
The DAQ supertype corresponds to the component of the
    DAQ and electronics setup that saves the waveform after
    it triggered
"""
abstract type DAQ end

"""
GenericDAQ is a dummy DAQ model that stores the waveform
    after it triggered, saving <baseline_length> samples before the
    trigger index, and in total <nsamples> samples from start to end
"""
@with_kw struct GenericDAQ <: DAQ
    "Number of samples i.e. length of the final stored waveform"
    nsamples::Int = 4000

    "Length of stored baseline in number of samples"
    baseline_length::Int = 400
end

"""
    GenericDAQ(sim_conf)

PropDict -> GenericDAQ

Construct a GenericDAQ struct based on given simulation configuration
"""
function GenericDAQ(sim_conf::PropDict)
    GenericDAQ(
        nsamples = sim_conf.setup.daq.nsamples,
        baseline_length = sim_conf.setup.daq.baseline_length
    )
end

"""
    DAQ(sim_conf)

PropDict -> <DAQ>

Construct a DAQ supertype struct based on given simulation configuration.
Type of returned instance depends on settings in <sim_conf>
Currently only one type of DAQ available (GenericDAQ),
    rendering this function temporarily redundant.
"""
function DAQ(sim_config::PropDict)
    GenericDAQ(sim_config)
end

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

"""
    simulate(wf, trigger_index, daq)

RDWaveform, Int, DAQ -> RDWaveform

Simulate how the DAQ stores the waveform after it receives a trigger.
"""
function simulate(wf::RDWaveform, trigger_index::Int, daq::GenericDAQ)
    # the stored waveform will start from zero
    # and go until the number of samples DAQ is configured to save
    ts = range(T(0)u"ns", step = step(wf.time), length = daq.nsamples)

    # in case the waveform didn't trigger (which is what 0 index stands for)...
    if(trigger_index == 0)
        # ...return all zeros
        return RDWaveform(ts, similar(wf.value[1:daq.nsamples])) # zeros
    else
        # ...otherwise store the waveform
        iStart = trigger_index - daq.baseline_length
        return RDWaveform(ts, wf.value[iStart : iStart + daq.nsamples-1])
    end
end


text:
using XGPaint
using Healpix
using HDF5
import ThreadsX
using JLD2

# ENV["SCRATCH"] = "/home/zequnl/scratch/"
@time halo_pos, halo_mass = read_halo_catalog_hdf5(
    joinpath(ENV["SCRATCH"],"websky_halos-light.hdf5"));

## Load halos from HDF5 files, establish a CIB model and cosmology
cosmo = get_cosmology(Float32, h = 0.7, OmegaM = 0.25)
radio_model = Radio_Sehgal2009{Float32}(a_0 = -1)
@time begin
    sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass);
end;

@save joinpath(scratch_dir, "before_match/catalog_$(freq).h5") sources

freqs = [
     "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", "63.9", "67.8",
     "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", "153", "164", "189", "210",
     "217", "232", "256", "275", "294", "306", "314", "340", "353", "375", "409", "467",
     "525", "545", "584", "643", "729", "817", "857", "906", "994", "1080"
]

# freqs = ["143.0"]

m = Map{Float64,RingOrder}(radio_model.nside)

function generate_maps()
    scratch_dir = ENV["SCRATCH"]
    println("SCRATCH: ", scratch_dir)

    for freq in freqs
        @time begin
#             flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II = paint!(
#                 m, parse(Float32, freq) * 1.0f9, radio_model, sources,
#                 return_fluxes=true)

#             flux = vcat(flux_I, flux_II)
#             redshift = vcat(redshift_I, redshift_II)
#             theta = vcat(θ_I, θ_II)
#             phi = vcat(ϕ_I, ϕ_II)
#             perm = sortperm(flux, rev=true, alg=ThreadsX.MergeSort)

#             h5open(joinpath(scratch_dir, "before_match/catalog_$(freq).h5"), "w") do file
#                 write(file, "flux", flux[perm])
#                 write(file, "redshift", redshift[perm])
#                 write(file, "theta", theta[perm])
#                 write(file, "phi", phi[perm])
#             end

            flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II = paint!(
                m, parse(Float32, freq) * 1.0f9, radio_model, sources,
                return_fluxes=true)

            flux = vcat(flux_I, flux_II)
            redshift = vcat(redshift_I, redshift_II)
            theta = vcat(θ_I, θ_II)
            phi = vcat(ϕ_I, ϕ_II)
            # perm = sortperm(flux, rev=true, alg=ThreadsX.MergeSort)

            h5open(joinpath(scratch_dir, "before_match/catalog_$(freq).h5"), "w") do file
                write(file, "flux", flux)
                write(file, "redshift", redshift)
                write(file, "theta", theta)
                write(file, "phi", phi)
            end

            map_path = joinpath(scratch_dir, "before_match/maps/radio$(freq).fits")
            Healpix.saveToFITS(m, "!$(map_path)")
        end
    end
end

generate_maps()
##

max_stars_repo_name: xzackli/WebskyRadioGalaxies

max_issues_repo_name: xzackli/WebskyRadioGalaxies

max_forks_repo_name: xzackli/WebskyRadioGalaxies


text:
module WaterIF97

    using Unitful
    using ThermoState
    using Roots

    using ThermoState.Types
    using ThermoState.QuickStates
    using ThermoState.StatePoints

    import ThermoState: pressure,temperature,mass,moles,molar_mass
    import ThermoState: mass_volume, mol_volume, total_volume
    import ThermoState: mass_density, mol_density
    import ThermoState: mass_enthalpy, mol_enthalpy,total_enthalpy
    import ThermoState: mass_gibbs, mol_gibbs, total_gibbs
    import ThermoState: mass_helmholtz, mol_helmholtz, total_helmholtz
    import ThermoState: mass_internal_energy, mol_internal_energy, total_internal_energy
    import ThermoState: mass_entropy, mol_entropy, total_entropy
    import ThermoState: mass_fraction, mol_fraction
    import ThermoState: mass_number, mol_number
    import ThermoState: ThermoModel
    import ThermoState: mol_cp,mass_cp,mol_cv,mol_cp,sound_speed

    struct IF97Region{T} <: ThermoModel end

    const IF97_R =  0.461526
    const IF97_Tc = 647.096      #K          Critical temperature of water
    const IF97_Pc = 22.064       #Pa         Critical pressure of water
    const IF97_ρc = 322.         #kg/m3      Critical density of water
    const IF97_T3 = 273.16       #K          Triple point temperature of water
    const IF97_P3 = 611.657E-6   #MPa        Triple point pressure of water
    const IF97_Mr = 18.01528     #kg/kmol    Molecular weight of water
    const IF97_Tc_inv = 1.0/647.096
    const IF97_ρc_inv = 1.0/322.
    struct IndustrialWater <: ThermoModel end

    temperature(::IndustrialWater,st::CriticalPoint,unit=u"K")  = convert_unit(u"K",unit,647.096)
    pressure(::IndustrialWater,st::CriticalPoint,unit=u"Pa")  = convert_unit(u"Pa",unit,2.2064e7)
    mol_volume(::IndustrialWater,st::CriticalPoint,unit=u"m^3/mol")  = convert_unit(u"m^3/mol",unit,5.594803726708074e-5)
    mol_density(::IndustrialWater,st::CriticalPoint,unit=u"mol/m^3")  = convert_unit(u"mol/m^3",unit,17873.72799560906)
    mass_density(::IndustrialWater,st::CriticalPoint,unit=u"kg/m^3")  = convert_unit(u"kg/m^3",unit,322.0)
    mass_volume(::IndustrialWater,st::CriticalPoint,unit=u"m^3/kg")  = convert_unit(u"m^3/kg",unit,0.003105590062111801)
    temperature(::IndustrialWater,st::TriplePoint,unit=u"K")  = convert_unit(u"K",unit,273.16)
    temperature(::IndustrialWater,st::NormalBoilingPoint,unit=u"K")  = convert_unit(u"K",unit,373.15)
    pressure(::IndustrialWater,st::TriplePoint,unit=u"Pa")  = convert_unit(u"Pa",unit,611.657)
    molecular_weight(::IndustrialWater) = 18.01528

    include("region1.jl")
    include("region2.jl")
    include("region3.jl")
    include("region4.jl")
    include("region5.jl")
    include("region_id.jl")
    include("industrialwater.jl")
    include("sat.jl")

    export IndustrialWater
    export mol_cp,mass_cp,mol_cv,mol_cp,sound_speed
end


text:
## PARTE 2
cd(@__DIR__)
input = map(line -> parse.(Bool, collect(line)), eachline("input.txt"));

##
function keep_only(pred, input)
    res = deepcopy(input)
    for j in axes(input[1], 1)
        length(res) == 1 && break
        pattern = pred(length(res))(foldl((s, row) -> s+row[j], res; init=0))
        filter!(row -> row[j] == pattern, res)
    end
    first(res)
end

##
oxygen_b = keep_only(d -> >=(cld(d, 2)), input)
CO2_b = keep_only(d -> <(cld(d, 2)), input)
n = length(oxygen_b)
oxygen, CO2 = foldl(0:n-1; init=(0,0)) do (o, c), i
    o |= oxygen_b[n-i] << i
    c |= CO2_b[n-i] << i
    o, c
end
oxygen * CO2


text:
suite["generators"] = BenchmarkGroup()

n=10; k=3; seed=17
suite["generators"]["rrg1"] = @benchmarkable random_regular_graph($n, $k, seed=$seed)
suite["generators"]["erdos1"] = @benchmarkable erdos_renyi($n, $k, seed=$seed)

n=100; k=3; seed=17
suite["generators"]["rrg2"] = @benchmarkable random_regular_graph($n, $k, seed=$seed)
suite["generators"]["erdos2"] = @benchmarkable erdos_renyi($n, $k, seed=$seed)

n=1000; k=3; seed=17
suite["generators"]["rrg3"] = @benchmarkable random_regular_graph($n, $k, seed=$seed)
suite["generators"]["erdos3"] = @benchmarkable erdos_renyi($n, $k, seed=$seed)