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

"""
    BatchCheckLayerAvailability()

Checks the availability of one or more image layers within a repository in a public registry. When an image is pushed to a repository, each image layer is checked to verify if it has been uploaded before. If it has been uploaded, then the image layer is skipped.  This operation is used by the Amazon ECR proxy and is not generally used by customers for pulling and pushing images. In most cases, you should use the docker CLI to pull, tag, and push images.

# Required Parameters
- `layerDigests`: The digests of the image layers to check.
- `repositoryName`: The name of the repository that is associated with the image layers to check.

# Optional Parameters
- `registryId`: The AWS account ID associated with the public registry that contains the image layers to check. If you do not specify a registry, the default public registry is assumed.
"""
batch_check_layer_availability(layerDigests, repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("BatchCheckLayerAvailability", Dict{String, Any}("layerDigests"=>layerDigests, "repositoryName"=>repositoryName); aws_config=aws_config)
batch_check_layer_availability(layerDigests, repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("BatchCheckLayerAvailability", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("layerDigests"=>layerDigests, "repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    BatchDeleteImage()

Deletes a list of specified images within a repository in a public registry. Images are specified with either an imageTag or imageDigest. You can remove a tag from an image by specifying the image's tag in your request. When you remove the last tag from an image, the image is deleted from your repository. You can completely delete an image (and all of its tags) by specifying the image's digest in your request.

# Required Parameters
- `imageIds`: A list of image ID references that correspond to images to delete. The format of the imageIds reference is imageTag=tag or imageDigest=digest.
- `repositoryName`: The repository in a public registry that contains the image to delete.

# Optional Parameters
- `registryId`: The AWS account ID associated with the registry that contains the image to delete. If you do not specify a registry, the default public registry is assumed.
"""
batch_delete_image(imageIds, repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("BatchDeleteImage", Dict{String, Any}("imageIds"=>imageIds, "repositoryName"=>repositoryName); aws_config=aws_config)
batch_delete_image(imageIds, repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("BatchDeleteImage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("imageIds"=>imageIds, "repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    CompleteLayerUpload()

Informs Amazon ECR that the image layer upload has completed for a specified public registry, repository name, and upload ID. You can optionally provide a sha256 digest of the image layer for data validation purposes. When an image is pushed, the CompleteLayerUpload API is called once per each new image layer to verify that the upload has completed.  This operation is used by the Amazon ECR proxy and is not generally used by customers for pulling and pushing images. In most cases, you should use the docker CLI to pull, tag, and push images.

# Required Parameters
- `layerDigests`: The sha256 digest of the image layer.
- `repositoryName`: The name of the repository in a public registry to associate with the image layer.
- `uploadId`: The upload ID from a previous InitiateLayerUpload operation to associate with the image layer.

# Optional Parameters
- `registryId`: The AWS account ID associated with the registry to which to upload layers. If you do not specify a registry, the default public registry is assumed.
"""
complete_layer_upload(layerDigests, repositoryName, uploadId; aws_config::AWSConfig=global_aws_config()) = ecr_public("CompleteLayerUpload", Dict{String, Any}("layerDigests"=>layerDigests, "repositoryName"=>repositoryName, "uploadId"=>uploadId); aws_config=aws_config)
complete_layer_upload(layerDigests, repositoryName, uploadId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("CompleteLayerUpload", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("layerDigests"=>layerDigests, "repositoryName"=>repositoryName, "uploadId"=>uploadId), args)); aws_config=aws_config)

"""
    CreateRepository()

Creates a repository in a public registry. For more information, see Amazon ECR repositories in the Amazon Elastic Container Registry User Guide.

# Required Parameters
- `repositoryName`: The name to use for the repository. This appears publicly in the Amazon ECR Public Gallery. The repository name may be specified on its own (such as nginx-web-app) or it can be prepended with a namespace to group the repository into a category (such as project-a/nginx-web-app).

# Optional Parameters
- `catalogData`: The details about the repository that are publicly visible in the Amazon ECR Public Gallery.
"""
create_repository(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("CreateRepository", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
create_repository(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("CreateRepository", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    DeleteRepository()

Deletes a repository in a public registry. If the repository contains images, you must either delete all images in the repository or use the force option which deletes all images on your behalf before deleting the repository.

# Required Parameters
- `repositoryName`: The name of the repository to delete.

# Optional Parameters
- `force`:  If a repository contains images, forces the deletion.
- `registryId`: The AWS account ID associated with the public registry that contains the repository to delete. If you do not specify a registry, the default public registry is assumed.
"""
delete_repository(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("DeleteRepository", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
delete_repository(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DeleteRepository", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    DeleteRepositoryPolicy()

Deletes the repository policy associated with the specified repository.

# Required Parameters
- `repositoryName`: The name of the repository that is associated with the repository policy to delete.

# Optional Parameters
- `registryId`: The AWS account ID associated with the public registry that contains the repository policy to delete. If you do not specify a registry, the default public registry is assumed.
"""
delete_repository_policy(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("DeleteRepositoryPolicy", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
delete_repository_policy(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DeleteRepositoryPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    DescribeImageTags()

Returns the image tag details for a repository in a public registry.

# Required Parameters
- `repositoryName`: The name of the repository that contains the image tag details to describe.

# Optional Parameters
- `maxResults`: The maximum number of repository results returned by DescribeImageTags in paginated output. When this parameter is used, DescribeImageTags only returns maxResults results in a single page along with a nextToken response element. The remaining results of the initial request can be seen by sending another DescribeImageTags request with the returned nextToken value. This value can be between 1 and 1000. If this parameter is not used, then DescribeImageTags returns up to 100 results and a nextToken value, if applicable. This option cannot be used when you specify images with imageIds.
- `nextToken`: The nextToken value returned from a previous paginated DescribeImageTags request where maxResults was used and the results exceeded the value of that parameter. Pagination continues from the end of the previous results that returned the nextToken value. This value is null when there are no more results to return. This option cannot be used when you specify images with imageIds.
- `registryId`: The AWS account ID associated with the public registry that contains the repository in which to describe images. If you do not specify a registry, the default public registry is assumed.
"""
describe_image_tags(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeImageTags", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
describe_image_tags(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeImageTags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    DescribeImages()

Returns metadata about the images in a repository in a public registry.  Beginning with Docker version 1.9, the Docker client compresses image layers before pushing them to a V2 Docker registry. The output of the docker images command shows the uncompressed image size, so it may return a larger image size than the image sizes returned by DescribeImages.

# Required Parameters
- `repositoryName`: The repository that contains the images to describe.

# Optional Parameters
- `imageIds`: The list of image IDs for the requested repository.
- `maxResults`: The maximum number of repository results returned by DescribeImages in paginated output. When this parameter is used, DescribeImages only returns maxResults results in a single page along with a nextToken response element. The remaining results of the initial request can be seen by sending another DescribeImages request with the returned nextToken value. This value can be between 1 and 1000. If this parameter is not used, then DescribeImages returns up to 100 results and a nextToken value, if applicable. This option cannot be used when you specify images with imageIds.
- `nextToken`: The nextToken value returned from a previous paginated DescribeImages request where maxResults was used and the results exceeded the value of that parameter. Pagination continues from the end of the previous results that returned the nextToken value. This value is null when there are no more results to return. This option cannot be used when you specify images with imageIds.
- `registryId`: The AWS account ID associated with the public registry that contains the repository in which to describe images. If you do not specify a registry, the default public registry is assumed.
"""
describe_images(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeImages", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
describe_images(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeImages", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    DescribeRegistries()

Returns details for a public registry.

# Optional Parameters
- `maxResults`: The maximum number of repository results returned by DescribeRegistries in paginated output. When this parameter is used, DescribeRegistries only returns maxResults results in a single page along with a nextToken response element. The remaining results of the initial request can be seen by sending another DescribeRegistries request with the returned nextToken value. This value can be between 1 and 1000. If this parameter is not used, then DescribeRegistries returns up to 100 results and a nextToken value, if applicable.
- `nextToken`: The nextToken value returned from a previous paginated DescribeRegistries request where maxResults was used and the results exceeded the value of that parameter. Pagination continues from the end of the previous results that returned the nextToken value. This value is null when there are no more results to return.  This token should be treated as an opaque identifier that is only used to retrieve the next items in a list and not for other programmatic purposes.
"""
describe_registries(; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeRegistries"; aws_config=aws_config)
describe_registries(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeRegistries", args; aws_config=aws_config)

"""
    DescribeRepositories()

Describes repositories in a public registry.

# Optional Parameters
- `maxResults`: The maximum number of repository results returned by DescribeRepositories in paginated output. When this parameter is used, DescribeRepositories only returns maxResults results in a single page along with a nextToken response element. The remaining results of the initial request can be seen by sending another DescribeRepositories request with the returned nextToken value. This value can be between 1 and 1000. If this parameter is not used, then DescribeRepositories returns up to 100 results and a nextToken value, if applicable. This option cannot be used when you specify repositories with repositoryNames.
- `nextToken`: The nextToken value returned from a previous paginated DescribeRepositories request where maxResults was used and the results exceeded the value of that parameter. Pagination continues from the end of the previous results that returned the nextToken value. This value is null when there are no more results to return. This option cannot be used when you specify repositories with repositoryNames.  This token should be treated as an opaque identifier that is only used to retrieve the next items in a list and not for other programmatic purposes.
- `registryId`: The AWS account ID associated with the registry that contains the repositories to be described. If you do not specify a registry, the default public registry is assumed.
- `repositoryNames`: A list of repositories to describe. If this parameter is omitted, then all repositories in a registry are described.
"""
describe_repositories(; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeRepositories"; aws_config=aws_config)
describe_repositories(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("DescribeRepositories", args; aws_config=aws_config)

"""
    GetAuthorizationToken()

Retrieves an authorization token. An authorization token represents your IAM authentication credentials and can be used to access any Amazon ECR registry that your IAM principal has access to. The authorization token is valid for 12 hours. This API requires the ecr-public:GetAuthorizationToken and sts:GetServiceBearerToken permissions.

"""
get_authorization_token(; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetAuthorizationToken"; aws_config=aws_config)
get_authorization_token(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetAuthorizationToken", args; aws_config=aws_config)

"""
    GetRegistryCatalogData()

Retrieves catalog metadata for a public registry.

"""
get_registry_catalog_data(; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRegistryCatalogData"; aws_config=aws_config)
get_registry_catalog_data(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRegistryCatalogData", args; aws_config=aws_config)

"""
    GetRepositoryCatalogData()

Retrieve catalog metadata for a repository in a public registry. This metadata is displayed publicly in the Amazon ECR Public Gallery.

# Required Parameters
- `repositoryName`: The name of the repository to retrieve the catalog metadata for.

# Optional Parameters
- `registryId`: The AWS account ID associated with the registry that contains the repositories to be described. If you do not specify a registry, the default public registry is assumed.
"""
get_repository_catalog_data(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRepositoryCatalogData", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
get_repository_catalog_data(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRepositoryCatalogData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    GetRepositoryPolicy()

Retrieves the repository policy for the specified repository.

# Required Parameters
- `repositoryName`: The name of the repository with the policy to retrieve.

# Optional Parameters
- `registryId`: The AWS account ID associated with the public registry that contains the repository. If you do not specify a registry, the default public registry is assumed.
"""
get_repository_policy(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRepositoryPolicy", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
get_repository_policy(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("GetRepositoryPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    InitiateLayerUpload()

Notifies Amazon ECR that you intend to upload an image layer. When an image is pushed, the InitiateLayerUpload API is called once per image layer that has not already been uploaded. Whether or not an image layer has been uploaded is determined by the BatchCheckLayerAvailability API action.  This operation is used by the Amazon ECR proxy and is not generally used by customers for pulling and pushing images. In most cases, you should use the docker CLI to pull, tag, and push images.

# Required Parameters
- `repositoryName`: The name of the repository to which you intend to upload layers.

# Optional Parameters
- `registryId`: The AWS account ID associated with the registry to which you intend to upload layers. If you do not specify a registry, the default public registry is assumed.
"""
initiate_layer_upload(repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("InitiateLayerUpload", Dict{String, Any}("repositoryName"=>repositoryName); aws_config=aws_config)
initiate_layer_upload(repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("InitiateLayerUpload", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    PutImage()

Creates or updates the image manifest and tags associated with an image. When an image is pushed and all new image layers have been uploaded, the PutImage API is called once to create or update the image manifest and the tags associated with the image.  This operation is used by the Amazon ECR proxy and is not generally used by customers for pulling and pushing images. In most cases, you should use the docker CLI to pull, tag, and push images.

# Required Parameters
- `imageManifest`: The image manifest corresponding to the image to be uploaded.
- `repositoryName`: The name of the repository in which to put the image.

# Optional Parameters
- `imageDigest`: The image digest of the image manifest corresponding to the image.
- `imageManifestMediaType`: The media type of the image manifest. If you push an image manifest that does not contain the mediaType field, you must specify the imageManifestMediaType in the request.
- `imageTag`: The tag to associate with the image. This parameter is required for images that use the Docker Image Manifest V2 Schema 2 or Open Container Initiative (OCI) formats.
- `registryId`: The AWS account ID associated with the public registry that contains the repository in which to put the image. If you do not specify a registry, the default public registry is assumed.
"""
put_image(imageManifest, repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutImage", Dict{String, Any}("imageManifest"=>imageManifest, "repositoryName"=>repositoryName); aws_config=aws_config)
put_image(imageManifest, repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutImage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("imageManifest"=>imageManifest, "repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    PutRegistryCatalogData()

Create or updates the catalog data for a public registry.

# Optional Parameters
- `displayName`: The display name for a public registry. The display name is shown as the repository author in the Amazon ECR Public Gallery.  The registry display name is only publicly visible in the Amazon ECR Public Gallery for verified accounts.
"""
put_registry_catalog_data(; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutRegistryCatalogData"; aws_config=aws_config)
put_registry_catalog_data(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutRegistryCatalogData", args; aws_config=aws_config)

"""
    PutRepositoryCatalogData()

Creates or updates the catalog data for a repository in a public registry.

# Required Parameters
- `catalogData`: An object containing the catalog data for a repository. This data is publicly visible in the Amazon ECR Public Gallery.
- `repositoryName`: The name of the repository to create or update the catalog data for.

# Optional Parameters
- `registryId`: The AWS account ID associated with the public registry the repository is in. If you do not specify a registry, the default public registry is assumed.
"""
put_repository_catalog_data(catalogData, repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutRepositoryCatalogData", Dict{String, Any}("catalogData"=>catalogData, "repositoryName"=>repositoryName); aws_config=aws_config)
put_repository_catalog_data(catalogData, repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("PutRepositoryCatalogData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("catalogData"=>catalogData, "repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    SetRepositoryPolicy()

Applies a repository policy to the specified public repository to control access permissions. For more information, see Amazon ECR Repository Policies in the Amazon Elastic Container Registry User Guide.

# Required Parameters
- `policyText`: The JSON repository policy text to apply to the repository. For more information, see Amazon ECR Repository Policies in the Amazon Elastic Container Registry User Guide.
- `repositoryName`: The name of the repository to receive the policy.

# Optional Parameters
- `force`: If the policy you are attempting to set on a repository policy would prevent you from setting another policy in the future, you must force the SetRepositoryPolicy operation. This is intended to prevent accidental repository lock outs.
- `registryId`: The AWS account ID associated with the registry that contains the repository. If you do not specify a registry, the default public registry is assumed.
"""
set_repository_policy(policyText, repositoryName; aws_config::AWSConfig=global_aws_config()) = ecr_public("SetRepositoryPolicy", Dict{String, Any}("policyText"=>policyText, "repositoryName"=>repositoryName); aws_config=aws_config)
set_repository_policy(policyText, repositoryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("SetRepositoryPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("policyText"=>policyText, "repositoryName"=>repositoryName), args)); aws_config=aws_config)

"""
    UploadLayerPart()

Uploads an image layer part to Amazon ECR. When an image is pushed, each new image layer is uploaded in parts. The maximum size of each image layer part can be 20971520 bytes (or about 20MB). The UploadLayerPart API is called once per each new image layer part.  This operation is used by the Amazon ECR proxy and is not generally used by customers for pulling and pushing images. In most cases, you should use the docker CLI to pull, tag, and push images.

# Required Parameters
- `layerPartBlob`: The base64-encoded layer part payload.
- `partFirstByte`: The position of the first byte of the layer part witin the overall image layer.
- `partLastByte`: The position of the last byte of the layer part within the overall image layer.
- `repositoryName`: The name of the repository to which you are uploading layer parts.
- `uploadId`: The upload ID from a previous InitiateLayerUpload operation to associate with the layer part upload.

# Optional Parameters
- `registryId`: The AWS account ID associated with the registry to which you are uploading layer parts. If you do not specify a registry, the default public registry is assumed.
"""
upload_layer_part(layerPartBlob, partFirstByte, partLastByte, repositoryName, uploadId; aws_config::AWSConfig=global_aws_config()) = ecr_public("UploadLayerPart", Dict{String, Any}("layerPartBlob"=>layerPartBlob, "partFirstByte"=>partFirstByte, "partLastByte"=>partLastByte, "repositoryName"=>repositoryName, "uploadId"=>uploadId); aws_config=aws_config)
upload_layer_part(layerPartBlob, partFirstByte, partLastByte, repositoryName, uploadId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ecr_public("UploadLayerPart", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("layerPartBlob"=>layerPartBlob, "partFirstByte"=>partFirstByte, "partLastByte"=>partLastByte, "repositoryName"=>repositoryName, "uploadId"=>uploadId), args)); aws_config=aws_config)


text:
module FeatureTest

    using FactCheck

    include(Pkg.dir("MethylUtils", "src","FeatureServices.jl"))
    path = Pkg.dir("MethylUtils", "testdata","moabs-format-1000lines.txt")

    sf = load_moabs_format_as_sequence_features(path)
    println(sf)

    rf = read_moabs_format_as_sequence_features(path)

    ensgene_path = Pkg.dir("MethylUtils","testdata","ensGene1000lines.txt.gz")
    ensgene_path = joinpath(homedir(), ".GenomeAnnotations.jl/annotation/mm10/ensGene.txt.gz")
    FT = make_ensgene_gene_feature_table( ensgene_path, "mm10")
    FR = FeatureRegions(FT)
    println( collect(keys(FR.regions)) )
    #local_store_path =  Pkg.dir("MethylUtils", "testdata","store", "annotation")
    #cgs = load_moabs(path)

    #SET_LOCAL_STORE_PATH(local_store_path)

    #facts("load_moabs") do
    #    @fact length(cgs) => 1000
    #end

    #transcripts = ensgene_dataframe("mm10","ensGene1000lines.txt.gz" )
    #println(transcripts)

    #facts("ensgene_dataframe") do
    #    @fact size(transcripts,1) => 1000
    #end

    #genes = ensgene_regions("mm10","ensGene1000lines.txt.gz" )

    # println( genes )
    # println("length: ", size(genes.dataframe, 1))
    # We want to test for ensgene_regions that it is making the kind of ensgene_regions
    # we want, i.e. the region covering all transcripts under a gene id

    # set config path back to nothing
    #SET_LOCAL_STORE_PATH("nothing")
end

ext: jl


text:
using Zygote
using ForwardDiff
using BenchmarkTools

function besselj(ν, z; atol=1e-8)
    k = 0
    s = (z/2)^ν / factorial(ν)
    out = s
    while abs(s) > atol
        k += 1
        s *= (-1) / k / (k+ν) * (z/2)^2
        out += s
    end
    out
end

function grad_besselj_manual(ν, z; atol=1e-8)
    (besselj(ν-1, z; atol=atol) - besselj(ν+1, z); atol=atol)/2
end

println("Benchmarking Julia")
display(@benchmark besselj(2, 1.0))
println("Benchmarking Manual")
display(@benchmark grad_besselj_manual(2, 1.0))
println("Benchmarking Zygote")
display(@benchmark Zygote.gradient(besselj, 2, 1.0))
println("Benchmarking ForwardDiff")
display(@benchmark ForwardDiff.derivative(x->besselj(2, x), 1.0))


text:
using PENBBControl
using Test

@testset "PENBBControl.jl" begin
    # Write your tests here.
end


text:
include("7.2.1 (a+b arccosh(c x))^n.jl")
include("7.2.2 (d x)^m (a+b arccosh(c x))^n.jl")
include("7.2.3 (d+e x^2)^p (a+b arccosh(c x))^n.jl")
include("7.2.4 (f x)^m (d+e x^2)^p (a+b arccosh(c x))^n.jl")
include("7.2.5 u (a+b arccosh(c x))^n.jl")
include("7.2.6 Miscellaneous inverse hyperbolic cosine.jl")

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.2 Inverse hyperbolic cosine/.jl


text:
export Spectra

module Spectra

using DelimitedFiles

abstract type AbstractSpectra{T} end

include("./opp.jl")

end

import .Spectra


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\src\ScoreDrivenERGM.jl\src\test\test_new_dyn_net_model copy.jl
#Created Date: Monday May 10th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Friday June 4th 2021 12:56:09 am
#Modified By:  Domenico Di Gangi
#-----
#Description: Test SDERGM_pml for different combinations of statistics
#-----
########

using Test
using ScoreDrivenERGM
using StatsBase
using Distributed

import ScoreDrivenERGM:StaticNets,DynNets, Utilities, ErgmRcall

# sample dgp
N = 50
T = 300

ergmTermsString = "edges + mutual + gwesp(decay = 0.25, fixed = TRUE, cutoff=10)"
model = DynNets.SdErgmPml(ergmTermsString, true)
nErgmPar = model.nErgmPar
@test nErgmPar == 3
dgpVals0 = DynNets.static_estimate(model, Int8.(rand(Bool, N, N)))
@test all(isfinite.(dgpVals0))

@testset " Mappings of static parameters beween restricted and unrestricted spaces" begin
    all_perm(xs, n) = vec(map(collect, Iterators.product(ntuple(_ -> xs, n)...)))

    listIndTvPar =  all_perm([true, false], nErgmPar)
    listIndTvPar = listIndTvPar[.!(all.( [.!indTvPar for indTvPar in listIndTvPar]))]

    for indTvPar in listIndTvPar

        vUnPar, _ = DynNets.starting_point_optim(model, dgpVals0)
        @test all(isfinite.(vUnPar))

        vResPar =  DynNets.restrict_all_par(model, vUnPar)
        @test all(isfinite.(vResPar))

        @test all(isapprox.(DynNets.unrestrict_all_par(model, vResPar), vUnPar, atol = 1e-8))
        @test all(isapprox.(DynNets.restrict_all_par(model, vUnPar), vResPar, atol = 1e-8))

        for i=1:10
            vResParRand = rand(length(vResPar))
            vUnParRand =  DynNets.unrestrict_all_par(model, vResParRand)
            @test all(isfinite.(vUnParRand))
            @test all(isapprox.(DynNets.restrict_all_par(model, vUnParRand), vResParRand, atol = 1e-8))
        end

    end

end

@testset " Sampling Change stats and estimates " begin
    T =100
    dgpParT = hcat(dgpVals0.*ones(nErgmPar, round(Int,T/2)),2 *dgpVals0.* ones(nErgmPar, round(Int, T/2)))
    # test reasonable sampling
    A_T = DynNets.sample_ergm_sequence(model, N, dgpParT, 1)
    eps = 0.01
    @test  eps < mean(A_T) < (1-eps)  # density
    @test !any(isnan.(A_T))
    @test all(isfinite.(A_T))

    obsT = [DynNets.stats_from_mat(model, A_T[:,:,t]) for t in 1:T ]
    @test all([all(isfinite.(obs)) for obs in obsT])

end

ENV["JULIA_DEBUG"] = nothing

# integrated version
ergmTermsString = "edges + mutual"
model = DynNets.SdErgmPml(ergmTermsString, true)
nErgmPar = model.nErgmPar
model.options["integrated"] = false
model.options["initMeth"] ="estimateFirstObs"# "uncMean"#

DynNets.get_option(model, "initMeth")

dgpVals0 = DynNets.static_estimate(model, Int8.(rand(Bool, N, N)))

dgpParT = hcat(dgpVals0.*ones(nErgmPar, round(Int,T/2)),2 *dgpVals0.* ones(nErgmPar, round(Int, T/2)))
# test reasonable sampling
A_T = DynNets.sample_ergm_sequence(model, N, dgpParT, 1)
obsT = [DynNets.stats_from_mat(model, A_T[:,:,t]) for t in 1:T ]

estSdResPar, conv_flag, UM_mple, ftot_0 = DynNets.estimate(model, N, obsT; indTvPar=model.indTvPar, show_trace = true )

vEstSdResParAll = DynNets.array_2_vec_all_par(model, estSdResPar, model.indTvPar)

vEstSdResPar, vConstPar = DynNets.divide_SD_par_from_const(model, vEstSdResParAll)

# vEstSdResPar[3:3:end] .= 0.001
fVecT_filt , target_fun_val_T, sVecT_filt = DynNets.score_driven_filter(model, N, obsT,  vEstSdResPar, model.indTvPar;ftot_0 = ftot_0, vConstPar=vConstPar)

DynNets.plot_filtered(model, N, fVecT_filt; parDgpTIn=dgpParT)

ss_filt = DynNets.estimate_single_snap_sequence(model, obsT)
DynNets.plot_filtered(model, N, ss_filt; parDgpTIn=dgpParT)

DynNets.seq_loglike_sd_filter(model, N, obsT,  vecUnPar, ftot_0Fun::Function)

indTvPar = model.indTvPar
vEstSdResParAll = DynNets.array_2_vec_all_par(model, estSdResPar, indTvPar)

vEstSdResPar, vConstPar = DynNets.divide_SD_par_from_const(model, vEstSdResParAll)

fVecT_filt , target_fun_val_T, sVecT_filt = DynNets.score_driven_filter(model, N, obsT,  vEstSdResPar, indTvPar;ftot_0 = ftot_0, vConstPar=vConstPar)

using PyPlot

# res_est = DynNets.estimate_and_filter(model, N, obsT; show_trace = true)
fig, ax = DynNets.plot_filtered(model, N, fVecT_filt)

est_SS = DynNets.estimate_single_snap_sequence(model, obsT)
fig, ax = DynNets.plot_filtered(model, N, est_SS, ax=ax, lineType = ".", lineColor="r")

# correctly specified filter

# estimate


text:
export l₂loss

function l₂loss(𝐲̂, 𝐲)
    feature_dims = 2:(ndims(𝐲)-1)

    loss = sum(.√(sum(abs2, 𝐲̂-𝐲, dims=feature_dims)))
    y_norm = sum(.√(sum(abs2, 𝐲, dims=feature_dims)))

    return loss / y_norm
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libinput"
version = v"1.18.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://gitlab.freedesktop.org/libinput/libinput", "2599b0d89935631176a2d6a3626e30c6bc1040db"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libinput/
mkdir build
cd build

meson -D libwacom=false -D tests=false -D documentation=false -D debug-gui=false ../ --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libinput", :libinput),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("eudev_jll"),
    Dependency("mtdev_jll"),
    Dependency("libevdev_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Dates, Documenter, TimeZones

makedocs(
    modules=[TimeZones],
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages=[
        "Introduction" => "index.md",
        "Types" => "types.md",
        "Converting" => "conversions.md",
        "Arithmetic" => "arithmetic.md",
        "Rounding" => "rounding.md",
        "Current Time" => "current.md",
        "Frequently Asked Questions" => "faq.md",
        "API – Public" => "api-public.md",
        "API – Private" => "api-private.md",
    ],
    sitename="TimeZones.jl",
    checkdocs=:exports,
    linkcheck=true,
    strict=true,
)


text:
# projection functions for the vectorized problem
module prj
	using ..prb
	using ..asm
	using LinearAlgebra

	# useful constants
	const nx = copy(prb.nx)
	const nu = copy(prb.nu)
	const nxnu = nx+nu
	const N = copy(prb.N)

	function project_Kdual!(w,κK)
		κK .= w[nx*(N-1)+1:end]
		w[nx*(N-1)+1:end] .= max.(κK,0.0)
	end
	function project_Kdual(w,κK)
		ww = deepcopy(w)
		project_Kdual!(ww,κK)
		return ww
	end
	function project_K!(w,κK)
		w[1:nx*(N-1)] .= 0.0
		κK .= w[nx*(N-1)+1:end]
		w[nx*(N-1)+1:end] .= max.(κK,0.0)
	end
	function project_K(w,κK)
		ww = deepcopy(w)
		project_K!(ww,κK)
		return ww
	end
	function project_Kpolar!(w,κK)
		w .= w .- project_K(w,κK)
	end
	function project_Kpolar(w,κK)
		ww = deepcopy(w)
		project_Kpolar!(ww,κK)
		return ww
	end
	function project_Z!(z,u,x)
		for j = 1:N-1
			u .= z[(j-1)*nxnu+1:(j-1)*nxnu+nu]
			asm.project_U!(u)
			z[(j-1)*nxnu+1:(j-1)*nxnu+nu] .= u

			x .= z[(j-1)*nxnu+nu+1:nu+(j-1)*nxnu+nx]
			asm.project_X!(x)
			z[(j-1)*nxnu+nu+1:nu+(j-1)*nxnu+nx] .= x
		end
	end
	function project_Z(z,u,x)
		zz = deepcopy(z)
		project_Z!(zz,u,x)
		return zz
	end

	const compute_relerr(z,zopt,nrm2_zopt) = (norm(z .- zopt)^2)/nrm2_zopt
	const compute_relfval(z,fopt) = abs( 0.5 * sum(z .* (asm.G*z)) + asm.h⋅z - fopt )/abs( fopt )
	const compute_rd2K(z,nrm_z_jump,κK) = (norm( ((asm.H*z .- asm.g) .* asm.scl_val_Hg) .- project_K(((asm.H*z .- asm.g) .* asm.scl_val_Hg),κK))^2)/nrm_z_jump

end

module jump

	using ..prb
	using ..asm
	using JuMP, MosekTools, ECOS, Gurobi, COSMO, SCS
	using LinearAlgebra

	# useful constants
	const nx = copy(prb.nx)
	const nu = copy(prb.nu)
	const nxnu = nx+nu
	const nxb2 = Int64(nx/2)
	const N = copy(prb.N)

	# initialize main variable
	z = randn(nxnu*(N-1))
	z_xu = randn(nxnu*(N-1))

	function solver_proxy!(z,slvr=:ecos,verbosity=true)
		if slvr == :ecos
			model= Model(ECOS.Optimizer)
			set_optimizer_attribute(model,"printlevel",0)
			set_optimizer_attribute(model,"feastol",prb.ϵ_pd)
			set_optimizer_attribute(model,"abstol",prb.ϵ_gap)
			set_optimizer_attribute(model,"reltol",prb.ϵ_gap)
		elseif slvr == :gurobi
			model= Model(Gurobi.Optimizer)
			set_optimizer_attribute(model,"Presolve",0)
			set_optimizer_attribute(model,"FeasibilityTol",prb.ϵ_pd)
			set_optimizer_attribute(model,"OptimalityTol",prb.ϵ_pd)
			set_optimizer_attribute(model,"BarConvTol",prb.ϵ_gap)
			set_optimizer_attribute(model,"BarQCPConvTol",prb.ϵ_gap)
		elseif slvr == :mosek
			model= Model(Mosek.Optimizer)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_DFEAS",prb.ϵ_pd)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_PFEAS",prb.ϵ_pd)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_REL_GAP",prb.ϵ_gap)
		elseif slvr == :cosmo
			model= Model(COSMO.Optimizer)
			set_optimizer_attribute(model,"eps_abs",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_rel",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_prim_inf",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_dual_inf",prb.ϵ_pd)
		elseif slvr == :scs
			model= Model(SCS.Optimizer)
			set_optimizer_attribute(model,"eps",prb.ϵ_pd)
		else
			error("Invalid solver choice for JuMP.")
		end
		if ~verbosity
			set_silent(model)
		end

		@variable(model,zz[1:nxnu*(N-1)])

		asm.set_constraint_jump_vectorized!(model,zz)

		# dynamics
		@constraint(model,asm.H[1:nx*(N-1),1:end]*zz .- asm.g[1:nx*(N-1)] .== 0)

		# affine inequality constraints
		@constraint(model,asm.H[nx*(N-1)+1:end,1:end]*zz .- asm.g[nx*(N-1)+1:end] .≥ 0)

		# cost function
		@objective(model,Min,0.5 * sum(zz .* (asm.G*zz)) + asm.h⋅zz)

		optimize!(model)

		exit_status = string(termination_status(model))
		print("$(string(slvr)) termination status: $(exit_status)\n")
		@assert (exit_status ∈ ("OPTIMAL","SLOW_PROGRESS","ALMOST_OPTIMAL")) "Problem not solved correctly."

		z .= value.(zz)[1:end]
	end

	solver!(solver_name=:ecos,verbosity=true) = solver_proxy!(z,solver_name,verbosity)

	function solver_xu_proxy!(z,slvr=:ecos,verbosity=true)
		if slvr == :ecos
			model= Model(ECOS.Optimizer)
			set_optimizer_attribute(model,"printlevel",0)
			set_optimizer_attribute(model,"feastol",prb.ϵ_pd)
			set_optimizer_attribute(model,"abstol",prb.ϵ_gap)
			set_optimizer_attribute(model,"reltol",prb.ϵ_gap)
		elseif slvr == :gurobi
			model= Model(Gurobi.Optimizer)
			set_optimizer_attribute(model,"Presolve",0)
			set_optimizer_attribute(model,"FeasibilityTol",prb.ϵ_pd)
			set_optimizer_attribute(model,"OptimalityTol",prb.ϵ_pd)
			set_optimizer_attribute(model,"BarConvTol",prb.ϵ_gap)
			set_optimizer_attribute(model,"BarQCPConvTol",prb.ϵ_gap)
		elseif slvr == :mosek
			model= Model(Mosek.Optimizer)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_DFEAS",prb.ϵ_pd)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_PFEAS",prb.ϵ_pd)
			set_optimizer_attribute(model,"INTPNT_CO_TOL_REL_GAP",prb.ϵ_gap)
		elseif slvr == :cosmo
			model= Model(COSMO.Optimizer)
			set_optimizer_attribute(model,"eps_abs",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_rel",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_prim_inf",prb.ϵ_pd)
			set_optimizer_attribute(model,"eps_dual_inf",prb.ϵ_pd)
		elseif slvr == :scs
			model= Model(SCS.Optimizer)
			set_optimizer_attribute(model,"eps",prb.ϵ_pd)
		else
			error("Invalid solver choice for JuMP.")
		end
		if ~verbosity
			set_silent(model)
		end

		@variable(model,x[1:nx,1:N])
		@variable(model,u[1:nu,1:N-1])
		@variable(model,β1[1:N-1] ≥ 0)        			# aux variable for defining objective
		@variable(model,β2[1:N-1] ≥ 0)					# aux variable for defining objective

		# dynamics constraint
		Inx = Array(Diagonal(ones(nx)))
		@constraint(model,[t=1:N-1],[Inx -prb.A[t] -prb.B[t] -Inx]*cat(x[:,t+1],x[:,t],u[:,t],prb.e[t],dims=1) .== 0)

		# initial condition
		@constraint(model,x[:,1] .== prb.x0)

		# set problem constraints
		asm.set_constraint_jump_decomposed!(model,x,u)

		# define objective terms
		@constraint(model,[t=1:N-1],[β1[t],(sqrt(prb.R[t])*(u[:,t] .- prb.uref[t]))...] in SecondOrderCone())
		@constraint(model,[t=1:N-1],[β2[t],(sqrt(prb.Q[t])*(x[:,t+1] .- prb.xref[t]))...] in SecondOrderCone())

		@objective(model,Min,0.5*sum(β1 .^ 2) + 0.5*sum(β2 .^ 2))

		optimize!(model)

		exit_status = string(termination_status(model))
		print("$(string(slvr)) termination status: $(exit_status)\n")
		@assert (exit_status ∈ ("OPTIMAL","SLOW_PROGRESS","ALMOST_OPTIMAL")) "Problem not solved correctly."

		z .= cat([cat(value.(u)[:,t],value.(x)[:,t+1],dims=1) for t=1:N-1]...,dims=1)
	end

	solver_xu!(solver_name=:ecos,verbosity=true) = solver_xu_proxy!(z_xu,solver_name,verbosity)

	compute_nrm_z() = norm(z)^2
	compute_fopt() = 0.5 * sum(z .* (asm.G*z)) + asm.h⋅z

end

# pipg solver and variable containers
module pipg

	using ..prb									# problem definition
	using ..prj
	using ..asm									# assembly of auxiliary parameters
	using ..jump
	using LinearAlgebra, StaticArrays, Plots
	pgfplotsx()

	# useful constants
	const nx = copy(prb.nx)
	const nu = copy(prb.nu)
	const nxnu = nx+nu
	const nH = copy(asm.nH)
	const N = copy(prb.N)

	const μb2σ = 0.5*asm.μ/asm.σ
	const μp2λ = asm.μ+2*asm.λ

	const nrm_z_jump = jump.compute_nrm_z()
	const fopt = jump.compute_fopt()

	# initialize main variables
	z = randn(nxnu*(N-1))
	v = randn(nx*(N-1)+nH)
	w = randn(nx*(N-1)+nH)

	zeq = randn(nxnu*(N-1))
	yeq = randn(nx*(N-1)+nH)
	veq = randn(nx*(N-1)+nH)
	weq = randn(nx*(N-1)+nH)

	# solution stats
	rd2K = randn(prb.kmax)
	rd2o = randn(prb.kmax)
	rfval = randn(prb.kmax)

	rd2K_restrt = randn(prb.kmax)
	rd2o_restrt = randn(prb.kmax)
	rfval_restrt = randn(prb.kmax)

	rd2K_eq = randn(prb.kmax)
	rd2o_eq = randn(prb.kmax)
	rfval_eq = randn(prb.kmax)

	# temporary variables
	κK = randn(nH)
	const κx = MVector{nx}(randn(nx))
	const κu = MVector{nu}(randn(nu))

	# pipg+ solver (version 1)
	function solver_v1_proxy!(z,v,w,κK,κu,κx,rd2o,rd2K,rfval,restart_idx=prb.kmax)

		for k = 1:prb.kmax
			# βᵏ = (k+1)*μb2σ
			βᵏ = (mod(k,restart_idx)+1)*μb2σ

			# feasibility measure
			rd2K[k] = prj.compute_rd2K(z,nrm_z_jump,κK)

			# relative objective value decrease
			# rfval[k] = prj.compute_relfval(z,fopt)
			rfval[k] = abs(z⋅(asm.G*z)+z⋅asm.h - jump.z⋅(asm.G*jump.z)-jump.z⋅asm.h)/abs(jump.z⋅(asm.G*jump.z)+jump.z⋅asm.h)

			v .= w .- βᵏ .* (asm.H*z .- asm.g)
			prj.project_Kdual!(v,κK)

			# αᵏ = 2/(k*asm.μ+μp2λ)
			αᵏ = 2/(mod(k,restart_idx)*asm.μ+μp2λ)

			z .= z .- αᵏ .* (asm.G*z .+ asm.h .- asm.HT*v)
			prj.project_Z!(z,κu,κx)

			w .= w .- βᵏ .* (asm.H*z .- asm.g)
			# prj.project_Kdual!(w,κK)

			rd2o[k] = prj.compute_relerr(z,jump.z,nrm_z_jump)

		end
		if restart_idx == prb.kmax
			print("PIPG+ relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
		else
			print("PIPG+ restarted w/ $(restart_idx) relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
		end

	end

	# pipg+ solver
	function solver_proxy!(z,v,w,κK,κu,κx,rd2o,rd2K,rfval,restart_idx=prb.kmax)

		zkm1 = randn(length(z))
		for k = 1:prb.kmax
			# βᵏ = (k+1)*μb2σ
			βᵏ = (mod(k,restart_idx)+1)*μb2σ

			# feasibility measure
			rd2K[k] = prj.compute_rd2K(z,nrm_z_jump,κK)

			# relative objective value decrease
			# rfval[k] = prj.compute_relfval(z,fopt)
			rfval[k] = abs(z⋅(asm.G*z)+z⋅asm.h - jump.z⋅(asm.G*jump.z)-jump.z⋅asm.h)/abs(jump.z⋅(asm.G*jump.z)+jump.z⋅asm.h)

			w .= v .+ βᵏ .* (asm.H*z .- asm.g)
			prj.project_Kpolar!(w,κK)

			# αᵏ = 2/(k*asm.μ+μp2λ)
			αᵏ = 2/(mod(k,restart_idx)*asm.μ+μp2λ)

			zkm1 .= z
			z .= z .- αᵏ .* (asm.G*z .+ asm.h .+ asm.HT*w)
			prj.project_Z!(z,κu,κx)

			v .= w .+ βᵏ .* (asm.H)*(z .- zkm1)

			rd2o[k] = prj.compute_relerr(z,jump.z,nrm_z_jump)

		end
		if restart_idx == prb.kmax
			print("PIPG+ relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
		else
			print("PIPG+ restarted w/ $(restart_idx) relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
		end

	end

	solver!(restart_idx=prb.kmax) = solver_proxy!(z,v,w,κK,κu,κx,rd2o,rd2K,rfval,restart_idx)

	# pipg solver (non-strongly convex formulation; affine equality constraint replaces the affine conic constraint)
	function solver_eq_proxy!(z,y,v,w,κK,κu,κx,rd2o,rd2K,rfval)
		β = 2.0
		α = 1/(β*asm.σ_eq+asm.λ)
		for k = 1:prb.kmax

			# feasibility measure
			rd2K[k] = prj.compute_rd2K(z,nrm_z_jump,κK)

			# relative objective value decrease
			# rfval[k] = prj.compute_relfval(z,fopt)
			rfval[k] = abs(z⋅(asm.G*z)+z⋅asm.h - jump.z⋅(asm.G*jump.z)-jump.z⋅asm.h)/abs(jump.z⋅(asm.G*jump.z)+jump.z⋅asm.h)

			v .= w .+ β .* (asm.H*z .- y .- asm.g)

			z .= z .- α .* (asm.G*z .+ asm.h .+ asm.HT*v)
			prj.project_Z!(z,κu,κx)

			y .= y .+ α .* v
			prj.project_K!(y,κK)

			w .= w .+ β .* (asm.H*z .- y .- asm.g)

			rd2o[k] = prj.compute_relerr(z,jump.z,nrm_z_jump)
		end

		print("PIPGeq relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")

	end

	solver_eq!() = solver_eq_proxy!(zeq,yeq,veq,weq,κK,κu,κx,rd2o_eq,rd2K_eq,rfval_eq)

	plot_solstat() = begin plt_d2o = plot(1:prb.kmax,rd2o,line=:solid,color=:magenta,title="Relative distance to optimum",lab="PIPG",lw=1.1,yaxis=:log,xlabel="Iterations")
			plot!(1:prb.kmax,rd2o_restrt,line=:solid,color=:red,lab="PIPG restarted",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2o_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(size=(400,500),legend=:none)
			plt_fval = plot(1:prb.kmax,rfval,line=:solid,color=:magenta,title="Relative objective value",lab="PIPG",lw=1.1,yaxis=:log,xlabel="Iterations")
			plot!(1:prb.kmax,rfval_restrt,line=:solid,color=:red,lab="PIPG restarted",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rfval_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(size=(400,500),legend=:none)
			plt_feas = plot(1:prb.kmax,rd2K,line=:solid,color=:magenta,title="Relative feasibility",lab="PIPG",lw=1.1,xlabel="Iterations",yaxis=:log)
			plot!(1:prb.kmax,rd2K_restrt,line=:solid,color=:red,lab="PIPG restarted",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2K_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(size=(400,500))
			λ = @layout [a b c]
			display(plot(plt_d2o,plt_fval,plt_feas,layout=λ,size=(4*1200/5,4*500/5)))
	end

	function reset_var!(flg=:both)
		if flg == :plus
			z .= randn(nxnu*(N-1))
			v .= randn(nx*(N-1)+nH)
			w .= randn(nx*(N-1)+nH)
		elseif flg == :eq
			zeq .= randn(nxnu*(N-1))
			yeq .= randn(nx*(N-1)+nH)
			veq .= randn(nx*(N-1)+nH)
			weq .= randn(nx*(N-1)+nH)
		elseif flg == :both
			z .= randn(nxnu*(N-1))
			v .= randn(nx*(N-1)+nH)
			w .= randn(nx*(N-1)+nH)

			zeq .= randn(nxnu*(N-1))
			yeq .= randn(nx*(N-1)+nH)
			veq .= randn(nx*(N-1)+nH)
			weq .= randn(nx*(N-1)+nH)
		else
			error("Invalid reset flag.")
		end
	end

end

# rival solvers
module rival

	using ..prb									# problem definition
	using ..prj
	using ..asm									# assembly of auxiliary parameters
	using ..jump
	using LinearAlgebra, StaticArrays, Plots, LaTeXStrings
	gr()
	default(xtickfontsize=11,ytickfontsize=11,ztickfontsize=11,legendfontsize=9,labelfontsize=11)

	# useful constants
	const nx = copy(prb.nx)
	const nu = copy(prb.nu)
	const nxnu = nx+nu
	const nH = copy(asm.nH)
	const N = copy(prb.N)

	const nrm_z_jump = jump.compute_nrm_z()
	const fopt = jump.compute_fopt()

	# initialize main variables
	z = randn(nxnu*(N-1))
	x = randn(nxnu*(N-1))
	y = randn(nx*(N-1)+nH)
	v = randn(nx*(N-1)+nH)
	w = randn(nx*(N-1)+nH)

	# solution stats
	rd2K = randn(prb.kmax)
	rd2o = randn(prb.kmax)
	rfval = randn(prb.kmax)
	nest_iter_count = zeros(prb.kmax)

	# temporary variables
	# temporary variables
	κK = randn(nH)
	const κx = MVector{nx}(randn(nx))
	const κu = MVector{nu}(randn(nu))
	κz = randn(nxnu*(N-1))

	# ADMM
	function solver_admm_proxy!(z,x,y,v,w,κz,κK,κx,κu,rd2o,rd2K,rfval,nest_iter_count)
		α = copy(prb.α_admm)
		λ̃ = asm.λ + α*asm.σ
		λ̃inv = 1/λ̃
		β = (sqrt(λ̃) - sqrt(asm.μ))/(sqrt(λ̃) + sqrt(asm.μ))

		inner_term_crit = 100.0
		k_nest = 0
		for k=1:prb.kmax

			# feasibility measure
			rd2K[k] = prj.compute_rd2K(z,nrm_z_jump,κK)

			# relative objective value decrease
			# rfval[k] = prj.compute_relfval(z,fopt)
			rfval[k] = abs(z⋅(asm.G*z)+z⋅asm.h - jump.z⋅(asm.G*jump.z)-jump.z⋅asm.h)/abs(jump.z⋅(asm.G*jump.z)+jump.z⋅asm.h)

			while inner_term_crit ≥ prb.ϵ_nest_admm
				v .= asm.H*x .- y .- asm.g .+ w

				κz .= z
				z .= x .- λ̃inv .* (asm.G*x .+ asm.h .+ α .* asm.HT*v)
				prj.project_Z!(z,κu,κx)

				x .= z .+ β .* (z .- κz)

				k_nest > prb.kmax_nest && break

				inner_term_crit = norm(z .- x)/norm(z)
				k_nest+=1
				nest_iter_count[k] += 1
			end
			k_nest = 0
			inner_term_crit = 100.0

			y .= asm.H*z .- asm.g .+ w
			prj.project_K!(y,κK)
			w .= w .+ asm.H*z .- y .- asm.g

			rd2o[k] = prj.compute_relerr(z,jump.z,nrm_z_jump)
		end
		print("ADMM relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
	end

	solver_admm!() = solver_admm_proxy!(z,x,y,v,w,κz,κK,κx,κu,rd2o,rd2K,rfval,nest_iter_count)

	z_2 = randn(nxnu*(N-1))
	y_2 = randn(nxnu*(N-1))
	w_2 = randn(nx*(N-1)+nH)

	# solution stats
	rd2K_2 = randn(prb.kmax)
	rd2o_2 = randn(prb.kmax)
	rfval_2 = randn(prb.kmax)

	# Chambolle and Pock variable step size
	function solver_cp_proxy!(z,y,w,κK,κx,κu,rd2o,rd2K,rfval)

	    βᵏ = prb.β⁰
	    γᵏ = prb.γ⁰
	    αᵏ = 1/(asm.λ+asm.σ*βᵏ)

		for k=1:prb.kmax

			# feasibility measure
			rd2K[k] = prj.compute_rd2K(z,nrm_z_jump,κK)

			# relative objective value decrease
			# rfval[k] = prj.compute_relfval(z,fopt)
			rfval[k] = abs(z⋅(asm.G*z)+z⋅asm.h - jump.z⋅(asm.G*jump.z)-jump.z⋅asm.h)/abs(jump.z⋅(asm.G*jump.z)+jump.z⋅asm.h)

	        w .= -1 .* prj.project_Kdual( -1 .* w .- βᵏ .* (asm.H*(z .+ γᵏ .* (z .- y)) .- asm.g),κK)
	        y .= z
	        α_step = αᵏ/(asm.μ*αᵏ+1)
	        z .= prj.project_Z(z .- α_step .* (asm.G*z .+ asm.h .+ asm.HT*w),κu,κx)

	        # step-size update
	        term_1 = (0.5*asm.λ*αᵏ + sqrt((1+asm.μ*αᵏ)*(1-asm.λ*αᵏ) + (0.5*asm.λ*αᵏ)^2 ))
	        γᵏ = term_1/(1+asm.μ*αᵏ)
	        αᵏ = αᵏ/term_1
	        βᵏ = (-asm.λ + 1/αᵏ)/asm.σ

			rd2o[k] = prj.compute_relerr(z,jump.z,nrm_z_jump)
		end
		print("Chambolle and Pock relative distance to optimum: $(prj.compute_relerr(z,jump.z,nrm_z_jump))\n")
	end

	solver_cp!() = solver_cp_proxy!(z_2,y_2,w_2,κK,κx,κu,rd2o_2,rd2K_2,rfval_2)

	function reset_var!(flg::Symbol)
		if flg == :admm
			z .= randn(nxnu*(N-1))
			x .= randn(nxnu*(N-1))
			y .= randn(nx*(N-1)+nH)
			v .= randn(nx*(N-1)+nH)
			w .= randn(nx*(N-1)+nH)
		elseif flg == :cp
			z_2 .= randn(nxnu*(N-1))
			y_2 .= randn(nxnu*(N-1))
			w_2 .= randn(nx*(N-1)+nH)
		else
			error("Invalid reset flag.")
		end
	end

	using ..pipg

	plot_solstat() = begin plt_d2o = plot(1:prb.kmax,pipg.rd2o,line=:solid,color=:magenta,title=L"\frac{\|z-z^{\star}\|_2^2}{\|z^{\star}\|_2^2}",lab="PIPG",lw=1.1,yaxis=:log,xlabel="Iterations")
			plot!(1:prb.kmax,pipg.rd2o_restrt,line=:solid,color=:red,lab="PIPG w/ restart",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,pipg.rd2o_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2o,line=:solid,color=:green,lab="ADMM",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2o_2,line=:solid,color=:orange,lab="PDHG",lw=1.1,yaxis=:log)
			plot!(size=(400,500),legend=:none,xlim=[1,prb.kmax],ylim=[1e-15,1e1])

			plt_rfval = plot(1:prb.kmax,pipg.rfval,line=:solid,color=:magenta,title=L"\frac{|f(z)-f(z^{\star})|}{|f(z^{\star})|}",lab="PIPG",lw=1.1,yaxis=:log,xlabel="Iterations")
			plot!(1:prb.kmax,pipg.rfval_restrt,line=:solid,color=:red,lab="PIPG w/ restart",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,pipg.rfval_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rfval,line=:solid,color=:green,lab="ADMM",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rfval_2,line=:solid,color=:orange,lab="PDHG",lw=1.1,yaxis=:log)
			plot!(size=(400,500),legend=:none,xlim=[1,prb.kmax],ylim=[1e-15,1e1])

			plt_nest_iter = plot(1:prb.kmax,nest_iter_count,line=:solid,color=:purple,title="Nesterov iteration count in ADMM",lw=1.1,xlabel="Outer iterations",lab="")
			plot!(1:prb.kmax,prb.kmax_nest .* ones(prb.kmax),line=:dash,lw=1.5,lab="Max iterations",color=:black,size=(1200,500))

			plt_feas = plot(1:prb.kmax,pipg.rd2K,line=:solid,color=:magenta,title=L"\frac{d_{K}(Hz-g)}{\|z^{\star}\|_2^2}",lab="PIPG",lw=1.1,xlabel="Iterations",yaxis=:log)
			plot!(1:prb.kmax,pipg.rd2K_restrt,line=:solid,color=:red,lab="PIPG w/ restart",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,pipg.rd2K_eq,line=:solid,color=:purple,lab="PIPGeq",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2K,line=:solid,color=:green,lab="ADMM",lw=1.1,yaxis=:log)
			plot!(1:prb.kmax,rd2K_2,line=:solid,color=:orange,lab="PDHG",lw=1.1,yaxis=:log)
			plot!(size=(400,500),legend=:topright,xlim=[1,prb.kmax],ylim=[1e-15,1e1])

			λ = @layout [a b c ; d]
			display(plot(plt_d2o,plt_rfval,plt_feas,plt_nest_iter,layout=λ,size=(4*1200/5,4*1000/5)))
	end

end

max_stars_repo_name: purnanandelango/pipg-demo

max_issues_repo_name: UW-ACL/pipg-demo

max_forks_repo_name: UW-ACL/pipg-demo


text:
using DSGE, ModelConstructors, JLD2, CSV, DataFrames, Dates, Test

regenerate_reference_output = false

if VERSION < v"1.5"
    ver = "111"
elseif VERSION < v"1.6"
    ver = "150"
else
    ver = "160"
end

# Data vintages and settings
fcast_date      = DSGE.quartertodate("2021-Q1")
date_fcast_end  = iterate_quarters(fcast_date, 40)
cond_type       = :full

# Flexible AIT Rule
tvcred_n_qtrs          = 24 # length of tv cred period
start_tvcred_level     = 0.
end_tvcred_level       = 1.
φ_π                    = 15.
φ_y                    = 3.
Thalf                  = 10.
ρ_smooth               = .75
pgap_init              = 0.125 # NOTE: positive pgap and ygap values imply
ygap_init              = 12.   # NEGATIVE initial pgap and ygaps (i.e. ygap_init = 12 implies
pgap_ygap_init_date    = Date(2020, 6, 30)  # an initial ygap of -12 in 2020:Q2)
start_zlb_date         = Date(2020, 12, 31) # start and end dates should be inclusive, i.e. first period w/zero rate
end_zlb_date           = Date(2023, 9, 30)  # is start_zlb_date, and the last period is end_zlb_date

# additional settings to implement Flexible AIT rule
flexait_custom = Dict{Symbol, Setting}(:forecast_horizons => Setting(:forecast_horizons,
                                                                     subtract_quarters(date_fcast_end, fcast_date)),
                                       :add_iid_cond_obs_gdp_meas_err =>
                                       Setting(:add_iid_cond_obs_gdp_meas_err, true),
                                       :add_iid_cond_obs_corepce_meas_err =>
                                       Setting(:add_iid_cond_obs_corepce_meas_err, true),
                                       :add_pseudo_corepce => Setting(:add_pseudo_corepce, true))

# Initialize model objects
m = Model1002("ss62", custom_settings = flexait_custom)
m <= Setting(:date_forecast_start,  fcast_date)
m <= Setting(:date_conditional_end, fcast_date)

m2p = get_setting(m, :model2para_regime)
m2p[:σ_condgdp] = Dict(1 => 1, 2 => 1, 3 => 1, 4 => 1, 5 => 1, 6 => 2)
m2p[:σ_condcorepce] = Dict(1 => 1, 2 => 1, 3 => 1, 4 => 1, 5 => 1, 6 => 2)
set_regime_val!(m[:σ_condgdp], 1, 0.)
set_regime_val!(m[:σ_condgdp], 2, 2. * m[:σ_gdp].value)
set_regime_val!(m[:σ_condcorepce], 1, 0.)
set_regime_val!(m[:σ_condcorepce], 2, 1.25 * m[:σ_corepce].value)

# Load data
df = CSV.read(joinpath(dirname(@__FILE__), "..", "reference", "ss62_modal_forecast_data.csv"), DataFrame)

# Set up Flexible AIT with a temporary ZLB. Note that df will be altered
start_reg_baseline = 6

θ = (φ_π = φ_π, φ_y = φ_y, Thalf = Thalf, pgap = pgap_init,
     ygap = ygap_init, ρ_smooth = ρ_smooth,
     historical_policy = default_policy()) # NamedTuple with parameters of the Flexible AIT policy
m <= Setting(:flexible_ait_φ_π, θ[:φ_π])
m <= Setting(:flexible_ait_φ_y, θ[:φ_y])
m <= Setting(:ait_Thalf, θ[:Thalf])
m <= Setting(:gdp_Thalf, θ[:Thalf])
m <= Setting(:pgap_value, θ[:pgap])
m <= Setting(:pgap_type, :flexible_ait)
m <= Setting(:ygap_value, θ[:ygap])
m <= Setting(:ygap_type, :flexible_ait)
m <= Setting(:flexible_ait_ρ_smooth, θ[:ρ_smooth])
m <= Setting(:alternative_policies, AltPolicy[θ[:historical_policy]])

function set_regime_vals_fnct(m, n)
    m2p = get_setting(m, :model2para_regime)

    for (k, v) in m2p
        for i in 7:n
            v[i] = 1 # For regimes 7 onward, use parameter regime 1
        end
    end

    m
end

# Set up credibility for ZLB and alternative policy
m <= Setting(:replace_eqcond, true)
m <= Setting(:gensys2, true)
m <= Setting(:uncertain_temporary_altpolicy, true)
m <= Setting(:uncertain_altpolicy, true)

# Find the first regime in which gensys2 should apply
gensys2_first_regime = findfirst([get_setting(m, :regime_dates)[i] .== start_zlb_date for i in 1:get_setting(m, :n_regimes)])
n_zlb_reg = DSGE.subtract_quarters(end_zlb_date, start_zlb_date) + 1
if end_zlb_date < get_setting(m, :regime_dates)[get_setting(m, :n_regimes)]
    set_regime_vals_fnct(m, get_setting(m, :n_regimes))
else
    # technically gensys2_first_regime + (n_zlb_reg - 1) + 1 b/c
    # gensys2_first_regime is the first regime for the temporary ZLB, so need to subtract 1
    # from n_zlb_reg, but need to also add 1 for the lift-off regime
    set_regime_vals_fnct(m, gensys2_first_regime + n_zlb_reg)
end

# Set up regime_eqcond_info
reg_dates = deepcopy(get_setting(m, :regime_dates))
regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
weights = [end_tvcred_level, 1. - end_tvcred_level]
for (regind, date) in zip(gensys2_first_regime:(n_zlb_reg - 1 + gensys2_first_regime), # See comments starting at line 75
                          DSGE.quarter_range(reg_dates[gensys2_first_regime],
                                                iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg - 1)))
    reg_dates[regind] = date
    regime_eqcond_info[regind] = DSGE.EqcondEntry(DSGE.zlb_rule(), weights) # weights are justdummy weights for now
end
reg_dates[n_zlb_reg + gensys2_first_regime] = iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg)
regime_eqcond_info[n_zlb_reg + gensys2_first_regime] = DSGE.EqcondEntry(flexible_ait(), weights)
m <= Setting(:regime_dates,               reg_dates)
m <= Setting(:regime_eqcond_info,         regime_eqcond_info)
m <= Setting(:temporary_altpolicy_length, n_zlb_reg)
m <= Setting(:zlb_rule_value,             .1)
m <= Setting(:temporary_altpolicy_names,  [:zlb_rule])

# Set up regime indices
setup_regime_switching_inds!(m; cond_type = cond_type)

# Set up time-varying credibility
tvcred_dates = (start_zlb_date, DSGE.iterate_quarters(start_zlb_date, tvcred_n_qtrs))
n_tvcred_reg = DSGE.subtract_quarters(tvcred_dates[2], tvcred_dates[1]) + 1 # number of regimes for time-varying credibility
reg_start = DSGE.subtract_quarters(tvcred_dates[1], start_zlb_date) + gensys2_first_regime # first regime for time-varying credibility
for (regind, date) in zip(reg_start:(reg_start + n_tvcred_reg),
                          DSGE.quarter_range(reg_dates[reg_start], DSGE.iterate_quarters(reg_dates[reg_start], n_tvcred_reg - 1)))
    reg_dates[regind] = date
    if date > end_zlb_date # just use dummy weights here, actual weights done later
        regime_eqcond_info[regind] = EqcondEntry(flexible_ait(), [end_tvcred_level, 1. - end_tvcred_level])
    end
end
m <= Setting(:regime_dates,             reg_dates)
m <= Setting(:regime_eqcond_info, regime_eqcond_info)
setup_regime_switching_inds!(m; cond_type = cond_type)
set_regime_vals_fnct(m, get_setting(m, :n_regimes))

# Set up TVIS information set
m <= Setting(:tvis_information_set, vcat([i:i for i in 1:(gensys2_first_regime - 1)],
                                         [i:17 for i in gensys2_first_regime:17],
                                         [i:i for i in 18:get_setting(m, :n_regimes)]))

# Actually populate weights
credvec = collect(range(start_tvcred_level, stop = end_tvcred_level, length = n_tvcred_reg))
for (i, k) in enumerate(sort!(collect(keys(regime_eqcond_info))))
    if tvcred_dates[1] <= reg_dates[k] <= tvcred_dates[2]
        get_setting(m, :regime_eqcond_info)[k].weights = [credvec[i], 1. - credvec[i]]
    end
end

m <= Setting(:zlb_rule_value, 0.1)

# Run forecasts
output_vars = [:histobs, :histpseudo, :forecastpseudo, :forecastobs]
modal_params = map(x -> x.value, m.parameters)

# Create two-rule as alternative policy
tworule_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
for (reg, eq_entry) in tworule_eqcond_info
    if reg >= start_reg_baseline + 5
        tworule_eqcond_info[reg] = DSGE.EqcondEntry(DSGE.default_policy(), [1.])
    end
    tworule_eqcond_info[reg].weights = [1.]
end
tworule = MultiPeriodAltPolicy(:two_rule, get_setting(m, :n_regimes), tworule_eqcond_info, gensys2 = true,
                               temporary_altpolicy_names = [:zlb_rule],
                               temporary_altpolicy_length = 6,
                               infoset = copy(get_setting(m, :tvis_information_set)))
delete!(DSGE.get_settings(m), :alternative_policies)
m <= Setting(:alternative_policies, DSGE.AbstractAltPolicy[tworule])

# Run forecast
fcast = DSGE.forecast_one_draw(m, :mode, cond_type, output_vars, modal_params, df,
                               regime_switching = true, n_regimes = get_setting(m, :n_regimes))

if regenerate_reference_output
    JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "ss62_modal_forecast_output_version=$(ver).jld2"), true, true, true, IOStream) do file
        write(file, "forecastobs", fcast[:forecastobs])
        write(file, "forecastpseudo", fcast[:forecastpseudo])
        write(file, "histpseudo", fcast[:histpseudo])
    end
end

@testset "Compare Modal Forecast for Model 1002 ss62 against reference" begin
    ref_out = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "ss62_modal_forecast_output_version=$(ver).jld2"), "r")
    @test maximum(abs.(fcast[:forecastobs] - ref_out["forecastobs"])) < 5e-3
    @test maximum(abs.(fcast[:forecastpseudo] - ref_out["forecastpseudo"])) < 5e-3
    @test maximum(abs.(fcast[:histpseudo] - ref_out["histpseudo"])) < 5e-3
end

nothing


text:
# E = exponent bits = 7 = 1 byte
# L = local var bits = 8
# global: x0, x1, ..., x(2^16 - 1)
# compression: global x -> local l_0, ..., l_(L-1)
# the max exp is 2^E-1
# each variable takes (E) = 8 bits
# all variabes take (L+1)*(E+1) = 256 bytes
# Monomial with capacity L+1: [tot degree bits; l_{L-1} exp bits, ..., l_1 exp bits, l_0 exp bits]
# Polynomial: [metadata, coeff n, monomial n, coeff n-1, monomial n-1, ..., coeff 0, monomial 0]
# Metadata: x_i <-> l_i bijection
# local to global: bitdict; Vector( global id ) ->
# global to local: hash map? searchsortedfirst the local to global?
# searchsorted is log2(2^L) = L; given that L is bounded, likely faster
# than hashmap (and takes up less space, as we reuse same mem)
# x_i < x_j <=> l_i < l_j
# (L + 1) * (E + 1)
# cld(ans, 64) # implement via + 63 and then >> 6
# align to 8 UInt64s
var_per_UInt64(::Val{E}) where {E} = 64 ÷ (new_E(Val(E)) + 1)
function calc_K(::Val{L},::Val{E}) where {L,E}
    @assert E <= 63
    # step 1: figure out how many E+1 we can fit inside a UInt64
    needed = cld((L + 1), var_per_UInt64(Val(E)))
    if needed <= 2
        needed
    elseif needed <= 4
        4
    else
        (needed + 7) & -8
    end
end
new_E(::Val{E}) where {E} = max(Base._nextpow2(1+E)-1, 7)

"""
    PackedMonomial{L,E}

Bit packed monomial with maximum of L variables and E bits of exponents.
"""
struct PackedMonomial{L,E,K} <: AbstractMonomial
    bits::NTuple{K,UInt64}
    function PackedMonomial{L,E}() where {L,E}
        EN = new_E(Val(E))
        K = calc_K(Val(L),Val(EN))
        PackedMonomial{L,EN,K}()
    end
    function PackedMonomial{L,E,K}() where {L,E,K}
        PackedMonomial{L,E,K}(ntuple(Ret(zero(UInt64)), Val(K)))
    end
    PackedMonomial{L,E}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
    PackedMonomial{L,E,K}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
end

PackedMonomial{L,E,K}(i::Integer) where {L,E,K} = PackedMonomial{L,E}(i)
function PackedMonomial{L,OE}(i::Integer) where {L,OE} # x_i
    @assert i < L
    E = new_E(Val(OE))
    t = PackedMonomial{L,E}().bits
    vpu = var_per_UInt64(Val(E))
    d, r = divrem(i+1, vpu)
    # o = [ 1, 0, 0, 0]
    o = one(UInt64) << ((E+1) * (vpu-1)) # `1` in `0`th typewriter (most significant) chunk
    b = o >> (r*(E+1)) # sets the bit for `i`th variable
    if d == 0 # degree and exponent on `i` are in same `UInt64`
        b |= o # set the total degree bit in same `UInt64` as `i`th monomial
    else
        t = Base.setindex(t, o, 1) # set the tot deg bit in first `UInt64`
    end
    t = Base.setindex(t, b, d+1)
    PackedMonomial{L,E}(t)
end

function degree(p::PackedMonomial{L,E}) where {L,E}
    (p.bits[1] >> ((E+1) * (var_per_UInt64(Val(E))-1))) % Int
end

function uint_type(::Val{E}) where {E}
    if E < 8
        return UInt8
    elseif E < 16
        return UInt16
    elseif E < 32
        return UInt32
    else
        return UInt64
    end
end
using SIMD

function calc_degree(p::PackedMonomial{L,E,K}) where {L,E,K}
    U = uint_type(Val(E))
    vpu = var_per_UInt64(Val(E))
    KS = vpu * K
    vp = reinterpret(Vec{KS,U},Vec{K,UInt64}(p.bits))
    sum(p) - p[vpu]
end
# we have 4 slots
# [ , , , ]
# we only need two of them: 1 for total order, 1 for a chunk
# var_per_UInt64 == 4
# L = 2
#
# [total degree, chunk, chunk, __]
# [__, total_degree, chunk, chunk] >> E+1
# Iterate right to left:
# [ a, b, c, d] mask top 3 => d; iter by shift right 1
# [ 0, a, b, c] mask top 3 => c; iter by shift right 1
# [ 0, 0, a, b] mask top 3 => b; iter by shift right 1
# [ 0, 0, 0, a] mask top 3 => a; iter by shift right 1
# number == 0 -> terminate
# Iterate left to right (typewriter)
# [ a, b, c, d] shift by 3 => a; iter by shift left 1
# [ b, c, d, 0] shift by 3 => b; iter by shift left 1
# [ c, d, 0, 0] shift by 3 => c; iter by shift left 1
# [ d, 0, 0, 0] shift by 3 => d; iter by shift left 1
# number == 0 -> terminate
function Base.gcd(p::T, q::T) where {L,E,K,T<:PackedMonomial{L,E,K}}
    U = uint_type(Val(E))
    vpu = var_per_UInt64(Val(E))
    KS = vpu * K
    vp = reinterpret(Vec{KS,U},Vec{K,UInt64}(p.bits))
    vq = reinterpret(Vec{KS,U},Vec{K,UInt64}(q.bits))
    g = min(vp, vq)
    total_degree = sum(g) - g[vpu]
    g = Base.setindex(g, total_degree, vpu)
    T(Tuple(reinterpret(Vec{K,UInt64}, g)))
end

# the first chunk is the total degree
Base.isone(x::PackedMonomial) = iszero(first(x.bits))
Base.isless(x::T, y::T) where {T<:PackedMonomial} = x.bits < y.bits
@noinline overflowed_error() = throw(Base.OverflowError("Try increasing E."))
@inline function check_zero_mask(::Val{E}) where {E}
    shifted = E
    x = one(UInt64) << E
    while shifted < 64
        x <<= 1
        x |= one(UInt64)
        x <<= E
        shifted += E
    end
    return x
end
@inline function zero_bits(x::UInt64, ::Val{E}) where {E}
    # every E+1st bit should be zero, so we create a mask, and then check for
    # zero
    (check_zero_mask(Val(E)) & x)
end
@inline function check_zero_bits(x::UInt64, ::Val{E}) where {E}
    # every E+1st bit should be zero, so we create a mask, and then check for
    # zero
    zero_bits(x, Val(E)) != zero(UInt64)
end

function Base.:*(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}}
    xys = _fmap(+, x.bits, y.bits)
    o = reduce_tup(|, _fmap(Base.Fix2(zero_bits, Val(E)), xys))

    o != zero(UInt64) && overflowed_error()
    return T(xys)
end
function Base.:/(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}}
    xys = _fmap(-, x.bits, y.bits)
    o = reduce_tup(|, _fmap(Base.Fix2(zero_bits, Val(E)), xys))

    #o != zero(UInt64) && overflowed_error()
    return T(xys), o != zero(UInt64)
end
function Base.:^(x::T, y::Integer) where {L,E,T<:PackedMonomial{L,E}}
    xys = _fmap(*, x.bits, y)
    o = reduce_tup(|, _fmap(Base.Fix2(zero_bits, Val(E)), xys))

    o != zero(UInt64) && overflowed_error()
    return T(xys)
end
zero_nondegreemask(::Val{7 }) = 0xff00000000000000
zero_nondegreemask(::Val{15}) = 0xffff000000000000
zero_nondegreemask(::Val{31}) = 0xffffffff00000000
zero_nondegreemask(::Val{63}) = 0xffffffffffffffff
function firstid(m::PackedMonomial{L,E,K}) where {L,E,K}
    b = m.bits[1] & (~zero_nondegreemask(Val(E)))
    if (b != zero(b))
        return (leading_zeros(b) ÷ (E+1)) - 1
    end
    vpu = var_per_UInt64(Val(E))
    b = vpu - 1
    for k in 2:K
        bk = m.bits[k]
        if (bk != zero(bk))
            return (leading_zeros(bk) ÷ (E+1)) + b
        end
        b += vpu
    end
    error("firstTermID should only be called if degree > 0.");
    return 0
end

function degree(m::PackedMonomial{L,E,K}, id) where {L,E,K}
    vpu = var_per_UInt64(Val(E))
    d = K == 1 ? 0 : (id + 1) ÷ vpu
    r = (id + 1) - (d * vpu)
    b = m.bits[d+1] << (r * (E + 1))
    return b >> ((E + 1) * (vpu - 1))
end

function rmid(x::T, id) where {L,E,K,T<:PackedMonomial{L,E,K}}
    bits = x.bits
    vpu = var_per_UInt64(Val(E))
    d = K == 1 ? 0 : (id + 1) ÷ vpu
    r = (id + 1) - (d * vpu)
    m = zero_nondegreemask(Val(E)) >> (r*(E+1))
    oldBits = bits[d+1]
    b = oldBits & (~m)
    remDegree = (oldBits & m) >> ((E+1)*(vpu-1-r))

    o = remDegree << ((E + 1) * (vpu - 1))
    if d != zero(d)
        bits = Base.setindex(bits, b, d+1)
        bits = Base.setindex(bits, bits[1] - o, 1)
    else
        bits = Base.setindex(bits, b - o, 1)
    end
    return T(bits)
end

function Base.show(io::IO, m::PackedMonomial{L,E}) where {L,E}
    tup = m.bits
    k = 0
    vpu = var_per_UInt64(Val(E))
    isone(m) && (print(io, '1'); return)
    for i in eachindex(tup)
        int = tup[i] << ((i == 1) * (E+1))
        for j in 1:vpu - (i == 1)
            exponent = int >> ((E+1)*(vpu-1))
            if exponent > 0
                print(io, 'x')
                print(io, int2lowerscript(k))
                exponent > 1 && print(io, int2superscript(exponent))
            end
            k += 1
            int <<= (E+1)
            L == k && break
        end
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ThreePartCost <: OperationalCost
        variable::VariableCost
        fixed::Float64
        startup::Float64
        shutdn::Float64
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure Operational Cost Data in Three parts fixed, variable cost and start - stop costs.

# Arguments
- `variable::VariableCost`: variable cost
- `fixed::Float64`: fixed cost
- `startup::Float64`: startup cost, validation range: (0, nothing), action if invalid: warn
- `shutdn::Float64`: shutdown cost, validation range: (0, nothing), action if invalid: warn
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ThreePartCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "fixed cost"
    fixed::Float64
    "startup cost"
    startup::Float64
    "shutdown cost"
    shutdn::Float64
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ThreePartCost(variable, fixed, startup, shutdn, forecasts=InfrastructureSystems.Forecasts(), )
    ThreePartCost(variable, fixed, startup, shutdn, forecasts, InfrastructureSystemsInternal(), )
end

function ThreePartCost(; variable, fixed, startup, shutdn, forecasts=InfrastructureSystems.Forecasts(), )
    ThreePartCost(variable, fixed, startup, shutdn, forecasts, )
end

# Constructor for demo purposes; non-functional.
function ThreePartCost(::Nothing)
    ThreePartCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
        startup=0.0,
        shutdn=0.0,
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get ThreePartCost variable."""
get_variable(value::ThreePartCost) = value.variable
"""Get ThreePartCost fixed."""
get_fixed(value::ThreePartCost) = value.fixed
"""Get ThreePartCost startup."""
get_startup(value::ThreePartCost) = value.startup
"""Get ThreePartCost shutdn."""
get_shutdn(value::ThreePartCost) = value.shutdn

InfrastructureSystems.get_forecasts(value::ThreePartCost) = value.forecasts
"""Get ThreePartCost internal."""
get_internal(value::ThreePartCost) = value.internal

"""Set ThreePartCost variable."""
set_variable!(value::ThreePartCost, val::VariableCost) = value.variable = val
"""Set ThreePartCost fixed."""
set_fixed!(value::ThreePartCost, val::Float64) = value.fixed = val
"""Set ThreePartCost startup."""
set_startup!(value::ThreePartCost, val::Float64) = value.startup = val
"""Set ThreePartCost shutdn."""
set_shutdn!(value::ThreePartCost, val::Float64) = value.shutdn = val

InfrastructureSystems.set_forecasts!(value::ThreePartCost, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set ThreePartCost internal."""
set_internal!(value::ThreePartCost, val::InfrastructureSystemsInternal) = value.internal = val


text:
# boot.m
# Luca Brugnolini 2014
# Based on Kilian code
# This function gives the Kilian bootstrap inverval for IRF by VAR

function boot(A, u, y, p, H)

# Input : Parameter estimates(intercept,slope), estimated residual, data
# Output: Bias-corrected bootstrap 95% confidence interval (Kilian 1999)

(T,K) = size(y);
nrep = 100;                              # number of bootstrap replications 2000
IRFrmat = zeros(nrep, K^2*(H+1));       # IRF for each bootstrap replication
u = u*sqrt((T-p)/(T-p-K*p-1));          # rescaling residual (Stine, JASA 1987)

#  Bootstrap simulation

 for j = 1:nrep

    # Draw block of initial pre-sample values
    yr = zeros(K,T);                             # bootstrap data
    pos = itrunc((rand()*(T-p+1)))+1;            # position of initial draw
    yr[:,1:p] = y[pos:pos+p-1,:]';               # drawing pre-sample obs

    # Draw innovations
    ur = zeros(K,T);                             # bootstrap innovations
    index = itrunc((rand(T-p)*(T-p))+1);         # index for innovation draws
    ur[:, p+1:T] = u[:, index];                  # drawing innovations

    # Recursively construct sample

    for i = p+1:T
            yr[:,i] = A[:,1]+ur[:,i];

        for jj = 1:p
            yr[:,i] = yr[:,i] + A[:,(jj-1)*K+2:jj*K+1]*yr[:,i-jj];

        end

    end

    yr = yr';

    for i=1:K                                   # demean yr bootstrap data
        yr[:,i] = yr[:,i] - mean(yr[:,i]);

    end

    # pbar = 4;
    # pr = varlagorder(yr,pbar);                 # lag-order selection for bootstrap data
    pr = p;
    (A1r,SIGMA1r) = olsvar(yr,pr);
    Ar = [A1r[:,2:end]; [eye(K*(pr-1)) zeros(K*(pr-1),K)]];   # Slope parameter estimates in companion form
    SIGMAr = [SIGMA1r zeros(K,K*pr-K); zeros(K*pr-K,K*pr)];   # variance of residuals in companion form

# Bias correction: if the largest root of the companion matrix
# is less than 1, do BIAS correction

if ~ any(abs(eigvals(Ar)).>=1)
     Ar = asybc(Ar,SIGMAr,T,K,pr);

end

    Ar = real(Ar)
    IRFr = irfvar1(Ar[1:K,1:K*pr], SIGMAr[1:K,1:K], pr, H);
    IRFrmat[j,:] = vec(IRFr')';

    end;                              # end bootstrap

# Calculate 95 perccent interval endpoints
CILv = zeros(1,size(IRFrmat,2));
CIHv = zeros(1,size(IRFrmat,2));

for i = 1:size(IRFrmat,2)
    CILv[:,i] = quantile(vec(IRFrmat[:,i]),0.025);
    CIHv[:,i] = quantile(vec(IRFrmat[:,i]),0.975);

end

CI = [CILv; CIHv];

                return CI

end


text:
# fixed value for Gh0 -- allow users to specify this value?
const Gh0 = 4/3  # this value results in pump gain crossing zero at the same flow rate as
                 # for efficiency, JJS 9/11/20

"""
Calculate slope and intercept for normalized power curve
"""
function power_coeffs(Gh0::Float64)
    # with fixed Gh0 value, then we have fixed values for slope and intercept for normalized
    # power curve

    # a least-squares fit -- would be more mathematically rigorouse to use an analytical
    # formula for the fit (TBD)
    N = 100
    Qhat = Array(range(0.5, 1.5, length=N))
    Phat = ((1 - Gh0)*Qhat.^2 .+ Gh0)./(2 .- Qhat)
    A = hcat(Qhat, ones(N))
    slope, intcpt = A\Phat
    return slope, intcpt
end

const pow_slp, pow_int = power_coeffs(Gh0)

"""
Estimate single-point head curve for a pump of 'POWER' type in the .inp file.
"""
function head_curve_from_power(power::Float64, efficiency::Float64, flows::Vector{Float64})
    flows_nz = flows[flows .> 0] # nonzero flow values
    Q_nominal = sum(flows_nz)/length(flows_nz) # m^3/s
    G_nominal = efficiency/100*power/(density*gravity*Q_nominal) # m
    head_curve = Vector{Tuple{Float64,Float64}}(undef,1)
    head_curve[1] = (Q_nominal,G_nominal)
    return head_curve
end

"""
Compute the pump parameters for normalized pump curves. Returns (Qbep, etabep, Gbep, slope,
intcpt).
"""
function norm_pump_params(pump::Dict{String,Any}, c_dict::Dict{String,Curve})

    # efficiency
    Qbep = NaN
    if pump["efficiency"] isa String
        eff_curve_tuples = c_dict[pump["efficiency"]].points
        eff_curve = array_from_tuples(eff_curve_tuples)

        # fit efficiency curve to determine Qbep and etabep
        A = hcat(eff_curve[:,1].^2, eff_curve[:,1])
        a, b = A\eff_curve[:,2]
        Qbep = -0.5*b/a
        etabep = -0.25*b^2/a
    else
        # presume the single-value for efficiency is at the BEP
        etabep = pump["efficiency"]
    end

    # head
    head_curve_tuples = c_dict[pump["head_curve_name"]].points
    head_curve = array_from_tuples(head_curve_tuples)
    if isnan(Qbep)
        # an efficiency curve was not provided, and so we must determine Qbep from the head
        # curve
        if size(head_curve)[1] == 1
            Qbep = head_curve[1,1]
            Gbep = head_curve[1,2]
        else
            A = hcat(head_curve[:,1].^2, ones(size(head_curve,1)))
            a, c = A\head_curve[:,2]
            Gbep = c/Gh0
            Qbep = sqrt((1 - Gh0)*Gbep/a)
        end
    else
        if size(head_curve)[1] == 1
            # probably rare that there is an efficiency curve and single-point or no head
            # curve -- nonetheless, it is a corner case that may be worth investigating
            # because Qbep would be overdetermined, and the two values may not agree
            #Qbep_alt = head_curve[1,1]
            Gbep = head_curve[1,2]
        else
            # again, we could determine Qbep from the head curve and check wether it agrees
            # with the value determined from the efficiency curve -- will presume that the
            # one from the efficiency curve is more accurate
            A = (1 - Gh0)/Qbep^2*head_curve[:,1].^2 .+ Gh0
            Gbep = A\head_curve[:,2]
        end
    end

    # Power at the BEP
    Pbep = density*gravity/(etabep/100)*Gbep*Qbep

    return Qbep, etabep, Gbep, Pbep
end

# function PumpPowerCoefs_Flow(wn_dict::Dict{Any,Any}, flow::Dict{String,Array{Float64,1}}, power::Dict{String,Array{Float64,1}})
#     aPumpPower = Dict{String,Float64}()
#     bPumpPower = Dict{String,Float64}()
#     for (name, pump) in wn_dict["pumps"]
#         coef = linear_fit(flow[name], power[name])
#         aPumpPower[name] = coef[1]
#         bPumpPower[name] = coef[2]
#         if bPumpPower[name] * minimum(flow[name]) + aPumpPower[name] < 0.0
#             aPumpPower[name], bPumpPower[name] = linear_fit([maximum(flow[name]), minimum(flow[name])], [maximum(power[name]), minimum(power[name])])
#         end
#     end
#     return(aPumpPower, bPumpPower)
# end

# function PumpPowerCoefs_Head(wn_dict::Dict{Any,Any},head::Dict{String,Array{Float64,1}}, power::Dict{String,Array{Float64,1}})
#     aPumpPower = Dict{String,Float64}()
#     bPumpPower = Dict{String,Float64}()
#     for (name, pump) in wn_dict["pumps"]
#         coef = linear_fit(head[name], power[name])
#         aPumpPower[name] = coef[1]
#         bPumpPower[name] = coef[2]
#         if bPumpPower[name] * minimum(head[name]) + aPumpPower[name] < 0.0
#             aPumpPower[name], bPumpPower[name] = linear_fit([maximum(head[name]), minimum(head[name])], [maximum(power[name]), minimum(power[name])])
#         end
#     end
#     return(aPumpPower, bPumpPower)
# end


text:
using PurlinLine
using StructuresKit

design_code = "AISI S100-16 nominal"

# Define the properties of each purlin segment along the line.
                #length, section_properties, material_properties
segments = [(25.0*12, 1, 1)]

# Define the purlin spacing.
spacing = 60;  #in.

# Define the roof slope.
roof_slope = 0.0;   #degrees

#Define the purlin cross-section type.
cross_section_dimensions = [("C", 0.059, 0.91, 2.5, 8.0, 2.5, 0.91, -90.0, 0.0, 90.0, 180.0, -90.0, 3*0.059, 3*0.059, 3*0.059, 3*0.059)]

material_properties = [(29500.0, 0.30, 55.0, 70.0)]

#type="screw-fastened", thickness, fastener spacing, fastener diameter, fastener_shear_strength
deck_details = ("screw-fastened", 0.018, 6.0, 0.212, 2.50)

deck_material_properties = (29500.0, 0.30, 55.0, 70.0)

frame_flange_width = 24.0

support_locations = [0.0, 25.0*12]

bridging_locations =[ ]

#Calculate purlin line design variables from user inputs and store them in the data structure.
purlin_line = PurlinLine.build(design_code, segments, spacing, roof_slope, cross_section_dimensions, material_properties, deck_details, deck_material_properties, frame_flange_width, support_locations, bridging_locations)
# purlin_line.applied_pressure = 10 / 1000 / 144 #psf to kips/in^2

purlin_line.cross_section_data[1].node_geometry

using Plots

plot(purlin_line.cross_section_data[1].node_geometry[:,1], purlin_line.cross_section_data[1].node_geometry[:,2], seriestype = :scatter)

plot(purlin_line.free_flange_cross_section_data[1].node_geometry[:,1], purlin_line.free_flange_cross_section_data[1].node_geometry[:,2], seriestype = :scatter)

purlin_line.loading_direction = "gravity"

#Perform a test to collapse.
purlin_line = PurlinLine.test(purlin_line)

purlin_line.applied_pressure * 1000 * 144

purlin_line.cross_section_data[1].plastic_section_properties

using Plots
plot(purlin_line.model.z, purlin_line.free_flange_model.u)

plot(purlin_line.model.z, purlin_line.free_flange_model.u)

plot(purlin_line.model.z, purlin_line.flexure_torsion_demand_to_capacity.action_Mxx)
plot!(purlin_line.model.z, purlin_line.flexure_torsion_demand_to_capacity.action_Myy)
plot!(purlin_line.model.z, purlin_line.flexure_torsion_demand_to_capacity.action_B)
plot!(purlin_line.model.z, purlin_line.flexure_torsion_demand_to_capacity.action_Myy_freeflange)

purlin_line.failure_limit_state

purlin_line.failure_location

plot(purlin_line.model.z, purlin_line.free_flange_internal_forces.P)

plot(purlin_line.model.z, purlin_line.free_flange_model.kx)

calculate_free_flange_shear_flow_properties(purlin_line)

#' Plot the purlin local buckling mode shape.
using StructuresKit
mode_index = 2
scale_x = 1.0
scale_y = 1.0
node = purlin_line.local_buckling_xx_pos[1].CUFSM_data.node
elem = purlin_line.local_buckling_xx_pos[1].CUFSM_data.elem
shapes = purlin_line.local_buckling_xx_pos[1].CUFSM_data.shapes

StructuresKit.CUFSM.view_multi_branch_section_mode_shape(node, elem, shapes, mode_index, scale_x, scale_y)

#' Plot the purlin local buckling signature curve.
num_lengths = length(purlin_line.local_buckling_xx_pos[1].CUFSM_data.lengths)
half_wavelengths = [purlin_line.local_buckling_xx_pos[1].CUFSM_data.curve[i,1][1] for i=1:num_lengths]
load_factors = [purlin_line.local_buckling_xx_pos[1].CUFSM_data.curve[i,1][2] for i=1:num_lengths]
plot(half_wavelengths, load_factors, legend = false)

#' Plot the purlin local buckling mode shape.
using StructuresKit
mode_index = 3
scale_x = 1.0
scale_y = 1.0
node = purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.node
elem = purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.elem
shapes = purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.shapes

StructuresKit.CUFSM.view_multi_branch_section_mode_shape(node, elem, shapes, mode_index, scale_x, scale_y)

#' Plot the purlin local buckling signature curve.
num_lengths = length(purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.lengths)
half_wavelengths = [purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.curve[i,1][1] for i=1:num_lengths]
load_factors = [purlin_line.distortional_buckling_xx_pos[1].CUFSM_data.curve[i,1][2] for i=1:num_lengths]
plot(half_wavelengths, load_factors, legend = false)

max_stars_repo_name: runtosolve/PurlinLine

max_issues_repo_name: runtosolve/PurlinLine

max_forks_repo_name: runtosolve/PurlinLine


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module SparsePerf
import Perftests: @perf, meta

## Sparse matrix performance
include("fem.jl")
@perf run_fem(256) meta("fem", "Finite Elements modeling")

# matsize = [log2(size(A,1)), log2(nnz(A))]
matsize = ([12,0], [12,6], [18,0], [18,9])

# indsize = log2(length(I))
indsize = [0, 8, 4]

function skinny_index(A,I)
    A[I,1]
end

for (m_log2,nz_log2) in matsize
    m = 2.^m_log2
    nz = 2.^nz_log2
    A = sprand(m,1,nz/m)
    for n_log2 in indsize
        n = 2.^n_log2
        I = rand(1:m,n)
        @perf A[I,1] meta("getindex_skinny", "m$(m)_nz$(nz)_n$(n)", "Sparse 2^$(m_log2)x1 matrix with 2^$(nz_log2) nonzero elements, indexed by 2^$(n_log2) elements")
    end
end

end # module


text:
abstract type AbstractInstruction end

struct Nop <: AbstractInstruction
    arg::Int64
end

struct Acc <: AbstractInstruction
    arg::Int64
end

struct Jmp <: AbstractInstruction
    arg::Int64
end

execute(acc::Int64, pc, ins::AbstractInstruction) = (acc, pc + 1)
execute(acc::Int64, pc, ins::Acc) = (acc += ins.arg, pc + 1)
execute(acc::Int64, pc, ins::Jmp) = (acc, pc += ins.arg)

function parse_instruction(str)
    rx = r"([a-z]+) ([\+|-][0-9]+)"
    m = match(rx, str)
    ins = m.captures[1]
    arg = parse(Int64, m.captures[2])

    ins == "nop" && return Nop(arg)
    ins == "acc" && return Acc(arg)
    ins == "jmp" && return Jmp(arg)

    error("could not parse instruction: '$(str)'")
end

parse_instructions(str) = map(parse_instruction, split(str, "\n", keepempty=false))

struct FlippedInstructionVector <: AbstractVector{AbstractInstruction}
    orig::Vector{AbstractInstruction}
    flipped_idx::Int64
end

Base.size(A::FlippedInstructionVector) = size(A.orig)

function Base.getindex(A::FlippedInstructionVector, i::Int)
    if i == A.flipped_idx
        arg = A.orig[i].arg
        return typeof(A.orig[i]) == Jmp ? Nop(arg) : Jmp(arg)
    end
    A.orig[i]
end

#=
parse_instruction("nop +0")
parse_instruction("jmp -3")
parse_instruction("acc +99")

ins_str = """
nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6
"""
ins = parse_instructions(ins_str)
=#

function has_duplicate_pc(pc, tracked_pcs)
    findfirst(isequal(pc), tracked_pcs) != nothing
end

function run_instructions(ins::AbstractVector{AbstractInstruction})
    tracked_pcs = Vector{Int64}(undef, 0)
    push!(tracked_pcs, 1)
    (acc, pc) = execute(0, 1, ins[1])
    while !has_duplicate_pc(pc, tracked_pcs)
        if pc == length(ins) + 1
            break
        end
        push!(tracked_pcs, pc)
        curr = ins[pc]
        (acc, pc) = execute(acc, pc, curr)
    end
    terminates = !has_duplicate_pc(pc, tracked_pcs)
    (acc, terminates)
end

function gen_flipped_instructions(ins::AbstractVector{AbstractInstruction})
    idxs = findall(x-> typeof(x) == Jmp || typeof(x) == Nop, ins)
    all_possible = Vector{FlippedInstructionVector}(undef, length(idxs))
    for (i, flipped_idx) in enumerate(idxs)
        all_possible[i] = FlippedInstructionVector(ins, flipped_idx)
    end
    all_possible
end

function find_fixed_instructions(ins::AbstractVector{AbstractInstruction})
    vv = gen_flipped_instructions(ins)
    for v in vv
        (acc, terminated) = run_instructions(v)
        if terminated
            return (acc, v)
        end
    end
    nothing
end

#=
run_instructions(ins)

gen_flipped_instructions(ins)
find_fixed_instructions(ins)
=#

function day08part1(ins)
    (acc, _) = run_instructions(ins)
    acc
end

function day08part2(ins)
    (acc, _) = find_fixed_instructions(ins)
    acc
end

function day08(file)
    s = read(file, String)
    ins = parse_instructions(s)
    (day08part1(ins), day08part2(ins))
end

#=
day08("input.txt")
=#


text:
# Single-line output in two-argument `show`
Base.show(io::IO, m::Monetary) = print(io, m / majorunit(m), currency(m))

# # Multi-line output with MIME types
# function Base.show(io::IO, ::MIME"text/plain", m::Monetary)
#     if get(io, :compact, false)
#         print(io, m / majorunit(m), currency(m))
#     else
#         print(io, format(m; styles=[:plain]))
#     end
# end

# function Base.show(io::IO, ::MIME"text/latex", m::Monetary)
#     print(io, string('$', format(m; styles=[:latex]), '$'))
# end


text:
module IonSim

using QuantumOptics
import OrderedCollections: OrderedDict

export OrderedDict, analytical
# Export some commonly used QuantumOptics.jl functions
export embed, ⊗, tensor, ⊕, dagger, normalize, normalize!, expect, tr, ptrace, tracenorm,
       tracedistance, entropy_vn, fidelity, diagonaljumps, dm, exp, norm

# used for copying composite types
Base.copy(x::T) where T = T([getfield(x, k) for k ∈ fieldnames(T)]...)
"""
    IonSimBasis
An abstract type for specialized bases, which are unique to IonSim.
"""
abstract type IonSimBasis <: Basis end
export IonSimBasis

include("constants.jl")
include("ions.jl")
include("vibrational_modes.jl")
include("lasers.jl")
include("ion_configurations.jl")
include("traps.jl")
include("operators.jl")
include("hamiltonians.jl")
include("time_evolution.jl")

module analytical
    include("analytic_functions.jl")
end

end  # main module


text:
macro bin_str(x)
    a = BitVector((i for i in x) .== '1')
    return a
end

tobig(x) = parse(BigInt, bytes2hex(reverse(x)), base=16)

function int2bytes(x::Integer)
    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end

    return reverse(hex2bytes(hex))
end

struct StaticBigInt{N} <: Integer
    x::NTuple{N, UInt8}
end

# function StaticBigInt(x::Integer)
#     bytes = int2bytes(x)
#     N = length(bytes)
#     sv = NTuple{N, UInt8}(bytes)
#     return StaticBigInt(sv)
# end

function StaticBigInt(x::Integer; n = bitlength(x))
    bytes = int2bytes(x)

    N = div(n, 8, RoundUp)

    append!(bytes, UInt8[0 for i in 1:N-length(bytes)])

    #N = length(bytes)

    sv = NTuple{N, UInt8}(bytes)
    return StaticBigInt(sv)
end

Base.convert(::Type{BigInt}, x::StaticBigInt) = tobig(x.x)
Base.BigInt(x::StaticBigInt) = tobig(x.x)

style(x, n) = "\33[1;$(n)m$x\33[0m"

function Base.show(io::IO, a::StaticBigInt)
    show(io, BigInt(a))
end

Base.display(x::StaticBigInt) = show(x)

nbits(x::UInt8) = 8
nbits(x::UInt64) = 64
nbits(x::UInt32) = 32
nbits(x::BigInt) = 64 * x.size

nbits(x::Int64) = 64
nbits(x::Int32) = 32
nbits(x::Int128) = 128

function tobin(a::Integer)
    s = BitVector(a << -i % Bool for i in 0:(nbits(a)-1))
    n = findlast(x->x==1, s)
    if n == nothing
        n = 0
    end
    return reverse(s[1:n])
end

function tobin(a::Integer, n::Integer)
    s = BitVector(a << -i % Bool for i in 0:(nbits(a)-1))

    if length(s) < n
        resize!(s, n)
        return reverse(s)
    else
        return reverse(s[1:n])
    end
end

function _hex2bytes(x::String)

    normalized = join(split(x, " "), "")

    N = length(normalized)

    if mod(N, 2) != 0
        normalized = "0" * normalized
    end

    return hex2bytes(normalized)
end

function bytes2bits(x::Vector{UInt8})
    bv = BitVector(u << -i % Bool for u in x for i in 7:-1:0)
    return bv
end

hex2bits(x::String) = bytes2bits(_hex2bytes(x))

struct StaticBitVector{N}
    len::Int
    chunks::NTuple{N, UInt64}
end

StaticBitVector(x::BitVector) = StaticBitVector(x.len, Tuple(x.chunks))

function Base.convert(::Type{BitVector}, x::StaticBitVector)

    a = BitVector(undef, x.len)
    a.chunks = UInt64[x.chunks...]

    return a
end

Base.show(io::IO, x::StaticBitVector) = print(io, join(i ? "1" : "0" for i in convert(BitVector, x)))

modinv(s, q) = mod(gcdx(s, q)[2], q)

function uint2string(x::Unsigned)
    bytes = reinterpret(UInt8, [x])

    N = findlast(x->x!=0, bytes)

    if N == nothing
        return String("")
    else
        return String(bytes[1:N])
    end
    #trimmed = bytes[1:findlast(x->x!=0, bytes)]
end

function string2uint(x::String)
    bytes = Vector{UInt8}(x)
    padded_bytes = append!(bytes, UInt8[0 for i in 1:16-length(bytes)])
    uint = reinterpret(UInt128, padded_bytes)
    return uint[1]
end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
# module Bukdu

abstract ApplicationView

immutable View <: ApplicationView
end


text:
# function numdecdigits(x)
#     d = 0
#     while round(x-round(x),15-d) != 0
#         x *= 10
#         d += 1
#     end
#     d
# end

# function numintdigits(x)
#     d = 0
#     while round(Int,x) != 0
#         x /= 10
#         d += 1
#     end
#     d
# end

# returns the exponent if x is expressed in scientific notation (base 10)
# e.g.
# 0.00123 = 1.23e-3
# so exponent is -3
function exponent10(x)
    if x == zero(x)
        return 0
    elseif x > zero(x)
        return floor(Int,log10(x))
    else
        return ceil(Int,log10(-x))
    end
end

# return the part of x without the exponent if expressed in scientific notation
# e.g.
# 0.00123 = 1.23e-3
# mantissa = 1.23
mantissa10(x) = x*exp10(-exponent10(x))

latex_expfmt(e) = string(" × 10^{",e,'}')
standard_expfmt(e) = string('e',e)

function formatter(x::AbstractFloat;
                   n::Int = 2,                      # number of significant figures
                   expfmt = latex_expfmt)           # exponent formatter

    (isnan(x) || isinf(x)) && return (string(x),' ',"")

    e = exponent10(x)                               # exponent
    m = round(x*exp10(-e), n-1)                     # mantissa

    # case 1: use scientific notation if number is extra long
    if e > 4 || e < -4
        s = string(format("{:.$(n-1)f}", m), expfmt(e))
        l,r = split(s, '.')
        return (l,'.',r)
    else
        xf = m*exp10(e)

        if n > e+1
            l,r = split(format("{:.$(n-e-1)f}", xf), '.')
            return (l,'.',r)
        elseif n == e+1
            l = format("{:d}", round(Int,xf))
            return (l,'.',"")
        else
            l = format("{:d}", round(Int,xf))
            return (l,' ',"")
        end
    end
end

formatter(x::Integer; kwargs...) = (string(x),' ',"")

formatter(x::AbstractString; kwargs...) = (x,' ',"")

function writetable(io::IO,
                    table::Matrix;
                    header = nothing,                    # column names
                    notes = nothing,
                    expfmt = latex_expfmt,
                    n = 2)                               # number of significant digits

    s = map(x -> formatter(x,n=n,expfmt=expfmt), table)

    if !is(notes, nothing)
        for (k,v) in notes
            i,j = k
            x = s[i,j]
            s[i,j] = (x[1],x[2],x[3] * " ($v)")
        end
    end

    ll,lr = map(x -> length(x[1]), s), map(x -> length(x[3]), s)                    # lengths of left, right parts of columns
    wl,wr = maximum(ll,1), maximum(lr,1)                                            # left, right column part widths
    w = wl+wr+1                                                                     # column widths

    addcols = cols -> string("| ", join(cols, " | "), " |")                         # combine columns in one row
    addparts = (s,wl,wr) -> string(lpad(s[1],wl), s[2], rpad(s[3],wr))              # combine left, middle, right parts of table entry

    s = [addparts(s[i,j],wl[j],wr[j]) for i in 1:size(s,1), j in 1:size(s,2)]

    if !is(header, nothing)
        wh = map(length, header)                                                    # column header widths
        s = [lpad(s[i,j], wh[j]) for i in 1:size(s,1), j in 1:size(s,2)]            # pad column to header width
        pheader = [rpad(header[i], w[i]) for i in 1:length(header)]                 # padded column header
        hs = addcols(pheader)                                                       # header string
        w = map(max,w,wh)
        hline = addcols([repeat("-", l) for l in w])
        ts = join([addcols(s[i,:]) for i in 1:size(s,1)], "\n")
        return write(io, hs, '\n', hline, '\n', ts)
    else
        ts = join([addcols(s[i,:]) for i in 1:size(s,1)], "\n")
        return write(io, ts)
    end
end


text:
comp = by(b,[:Protocol,:Wall,:MouseID]) do dd
    m = fit(FlippingModel.PoissonLapseUniform, dd)
    T, r, lapse = params(m)
    (T=T, r=r, lapse=lapse, model=m)
end

Simulated = by(comp,[:Protocol,:Wall,:MouseID]) do dd
    (simulation = FlippingModel.simulate(dd[1,:model], 500),)
end
##
comp = by(b,[:Protocol,:Wall,:MouseID]) do dd
    m = fit(FlippingModel.PoissonLapseExponential, dd)
    T, r, lapse = params(m)
    (T=T, r=r, lapse=lapse, model=m)
end
##
using Distributions, StatsBase, StructArrays

cm = countmap(StructArray((Omissions_plus_one = b.Omissions_plus_one, Leave = b.Leave)))
pE = FlippingModel.init(PoissonLapseExponential,cm)

##
e = Exponential(3)
g = Gamma(10,1/2)
plot([cdf(g,x) - cdf(g,x-1) for x in 1:20])
plot!([cdf(e,x) - cdf(e,x-1) for x in 1:20])
plot!([cdf(pE,om) - cdf(pE,om-1) for om in 1:20])
##
plot([cdf(g,x) for x in 1:20])
plot!([cdf(e,x) for x in 1:20])
plot!([cdf(pE,om) for om in 1:20])
##