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

"""
    CreateContainer()

Creates a storage container to hold objects. A container is similar to a bucket in the Amazon S3 service.

# Required Parameters
- `ContainerName`: The name for the container. The name must be from 1 to 255 characters. Container names must be unique to your AWS account within a specific region. As an example, you could create a container named movies in every region, as long as you don’t have an existing container with that name.

# Optional Parameters
- `Tags`: An array of key:value pairs that you define. These values can be anything that you want. Typically, the tag key represents a category (such as \"environment\") and the tag value represents a specific value within that category (such as \"test,\" \"development,\" or \"production\"). You can add up to 50 tags to each container. For more information about tagging, including naming and usage conventions, see Tagging Resources in MediaStore.
"""
create_container(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("CreateContainer", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
create_container(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("CreateContainer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DeleteContainer()

Deletes the specified container. Before you make a DeleteContainer request, delete any objects in the container or in any folders in the container. You can delete only empty containers.

# Required Parameters
- `ContainerName`: The name of the container to delete.

"""
delete_container(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteContainer", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
delete_container(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteContainer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DeleteContainerPolicy()

Deletes the access policy that is associated with the specified container.

# Required Parameters
- `ContainerName`: The name of the container that holds the policy.

"""
delete_container_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteContainerPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
delete_container_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteContainerPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DeleteCorsPolicy()

Deletes the cross-origin resource sharing (CORS) configuration information that is set for the container. To use this operation, you must have permission to perform the MediaStore:DeleteCorsPolicy action. The container owner has this permission by default and can grant this permission to others.

# Required Parameters
- `ContainerName`: The name of the container to remove the policy from.

"""
delete_cors_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteCorsPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
delete_cors_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteCorsPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DeleteLifecyclePolicy()

Removes an object lifecycle policy from a container. It takes up to 20 minutes for the change to take effect.

# Required Parameters
- `ContainerName`: The name of the container that holds the object lifecycle policy.

"""
delete_lifecycle_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteLifecyclePolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
delete_lifecycle_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteLifecyclePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DeleteMetricPolicy()

Deletes the metric policy that is associated with the specified container. If there is no metric policy associated with the container, MediaStore doesn't send metrics to CloudWatch.

# Required Parameters
- `ContainerName`: The name of the container that is associated with the metric policy that you want to delete.

"""
delete_metric_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteMetricPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
delete_metric_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DeleteMetricPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    DescribeContainer()

Retrieves the properties of the requested container. This request is commonly used to retrieve the endpoint of a container. An endpoint is a value assigned by the service when a new container is created. A container's endpoint does not change after it has been assigned. The DescribeContainer request returns a single Container object based on ContainerName. To return all Container objects that are associated with a specified AWS account, use ListContainers.

# Optional Parameters
- `ContainerName`: The name of the container to query.
"""
describe_container(; aws_config::AWSConfig=global_aws_config()) = mediastore("DescribeContainer"; aws_config=aws_config)
describe_container(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("DescribeContainer", args; aws_config=aws_config)

"""
    GetContainerPolicy()

Retrieves the access policy for the specified container. For information about the data that is included in an access policy, see the AWS Identity and Access Management User Guide.

# Required Parameters
- `ContainerName`: The name of the container.

"""
get_container_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("GetContainerPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
get_container_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("GetContainerPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    GetCorsPolicy()

Returns the cross-origin resource sharing (CORS) configuration information that is set for the container. To use this operation, you must have permission to perform the MediaStore:GetCorsPolicy action. By default, the container owner has this permission and can grant it to others.

# Required Parameters
- `ContainerName`: The name of the container that the policy is assigned to.

"""
get_cors_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("GetCorsPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
get_cors_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("GetCorsPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    GetLifecyclePolicy()

Retrieves the object lifecycle policy that is assigned to a container.

# Required Parameters
- `ContainerName`: The name of the container that the object lifecycle policy is assigned to.

"""
get_lifecycle_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("GetLifecyclePolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
get_lifecycle_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("GetLifecyclePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    GetMetricPolicy()

Returns the metric policy for the specified container.

# Required Parameters
- `ContainerName`: The name of the container that is associated with the metric policy.

"""
get_metric_policy(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("GetMetricPolicy", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
get_metric_policy(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("GetMetricPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    ListContainers()

Lists the properties of all containers in AWS Elemental MediaStore.  You can query to receive all the containers in one response. Or you can include the MaxResults parameter to receive a limited number of containers in each response. In this case, the response includes a token. To get the next set of containers, send the command again, this time with the NextToken parameter (with the returned token as its value). The next set of responses appears, with a token if there are still more containers to receive.  See also DescribeContainer, which gets the properties of one container.

# Optional Parameters
- `MaxResults`: Enter the maximum number of containers in the response. Use from 1 to 255 characters.
- `NextToken`: Only if you used MaxResults in the first command, enter the token (which was included in the previous response) to obtain the next set of containers. This token is included in a response only if there actually are more containers to list.
"""
list_containers(; aws_config::AWSConfig=global_aws_config()) = mediastore("ListContainers"; aws_config=aws_config)
list_containers(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("ListContainers", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of the tags assigned to the specified container.

# Required Parameters
- `Resource`: The Amazon Resource Name (ARN) for the container.

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

"""
    PutContainerPolicy()

Creates an access policy for the specified container to restrict the users and clients that can access it. For information about the data that is included in an access policy, see the AWS Identity and Access Management User Guide. For this release of the REST API, you can create only one policy for a container. If you enter PutContainerPolicy twice, the second command modifies the existing policy.

# Required Parameters
- `ContainerName`: The name of the container.
- `Policy`: The contents of the policy, which includes the following:    One Version tag   One Statement tag that contains the standard tags for the policy.

"""
put_container_policy(ContainerName, Policy; aws_config::AWSConfig=global_aws_config()) = mediastore("PutContainerPolicy", Dict{String, Any}("ContainerName"=>ContainerName, "Policy"=>Policy); aws_config=aws_config)
put_container_policy(ContainerName, Policy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("PutContainerPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName, "Policy"=>Policy), args)); aws_config=aws_config)

"""
    PutCorsPolicy()

Sets the cross-origin resource sharing (CORS) configuration on a container so that the container can service cross-origin requests. For example, you might want to enable a request whose origin is http://www.example.com to access your AWS Elemental MediaStore container at my.example.container.com by using the browser's XMLHttpRequest capability. To enable CORS on a container, you attach a CORS policy to the container. In the CORS policy, you configure rules that identify origins and the HTTP methods that can be executed on your container. The policy can contain up to 398,000 characters. You can add up to 100 rules to a CORS policy. If more than one rule applies, the service uses the first applicable rule listed. To learn more about CORS, see Cross-Origin Resource Sharing (CORS) in AWS Elemental MediaStore.

# Required Parameters
- `ContainerName`: The name of the container that you want to assign the CORS policy to.
- `CorsPolicy`: The CORS policy to apply to the container.

"""
put_cors_policy(ContainerName, CorsPolicy; aws_config::AWSConfig=global_aws_config()) = mediastore("PutCorsPolicy", Dict{String, Any}("ContainerName"=>ContainerName, "CorsPolicy"=>CorsPolicy); aws_config=aws_config)
put_cors_policy(ContainerName, CorsPolicy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("PutCorsPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName, "CorsPolicy"=>CorsPolicy), args)); aws_config=aws_config)

"""
    PutLifecyclePolicy()

Writes an object lifecycle policy to a container. If the container already has an object lifecycle policy, the service replaces the existing policy with the new policy. It takes up to 20 minutes for the change to take effect. For information about how to construct an object lifecycle policy, see Components of an Object Lifecycle Policy.

# Required Parameters
- `ContainerName`: The name of the container that you want to assign the object lifecycle policy to.
- `LifecyclePolicy`: The object lifecycle policy to apply to the container.

"""
put_lifecycle_policy(ContainerName, LifecyclePolicy; aws_config::AWSConfig=global_aws_config()) = mediastore("PutLifecyclePolicy", Dict{String, Any}("ContainerName"=>ContainerName, "LifecyclePolicy"=>LifecyclePolicy); aws_config=aws_config)
put_lifecycle_policy(ContainerName, LifecyclePolicy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("PutLifecyclePolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName, "LifecyclePolicy"=>LifecyclePolicy), args)); aws_config=aws_config)

"""
    PutMetricPolicy()

The metric policy that you want to add to the container. A metric policy allows AWS Elemental MediaStore to send metrics to Amazon CloudWatch. It takes up to 20 minutes for the new policy to take effect.

# Required Parameters
- `ContainerName`: The name of the container that you want to add the metric policy to.
- `MetricPolicy`: The metric policy that you want to associate with the container. In the policy, you must indicate whether you want MediaStore to send container-level metrics. You can also include up to five rules to define groups of objects that you want MediaStore to send object-level metrics for. If you include rules in the policy, construct each rule with both of the following:   An object group that defines which objects to include in the group. The definition can be a path or a file name, but it can't have more than 900 characters. Valid characters are: a-z, A-Z, 0-9, _ (underscore), = (equal), : (colon), . (period), - (hyphen), ~ (tilde), / (forward slash), and * (asterisk). Wildcards (*) are acceptable.   An object group name that allows you to refer to the object group. The name can't have more than 30 characters. Valid characters are: a-z, A-Z, 0-9, and _ (underscore).

"""
put_metric_policy(ContainerName, MetricPolicy; aws_config::AWSConfig=global_aws_config()) = mediastore("PutMetricPolicy", Dict{String, Any}("ContainerName"=>ContainerName, "MetricPolicy"=>MetricPolicy); aws_config=aws_config)
put_metric_policy(ContainerName, MetricPolicy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("PutMetricPolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName, "MetricPolicy"=>MetricPolicy), args)); aws_config=aws_config)

"""
    StartAccessLogging()

Starts access logging on the specified container. When you enable access logging on a container, MediaStore delivers access logs for objects stored in that container to Amazon CloudWatch Logs.

# Required Parameters
- `ContainerName`: The name of the container that you want to start access logging on.

"""
start_access_logging(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("StartAccessLogging", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
start_access_logging(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("StartAccessLogging", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    StopAccessLogging()

Stops access logging on the specified container. When you stop access logging on a container, MediaStore stops sending access logs to Amazon CloudWatch Logs. These access logs are not saved and are not retrievable.

# Required Parameters
- `ContainerName`: The name of the container that you want to stop access logging on.

"""
stop_access_logging(ContainerName; aws_config::AWSConfig=global_aws_config()) = mediastore("StopAccessLogging", Dict{String, Any}("ContainerName"=>ContainerName); aws_config=aws_config)
stop_access_logging(ContainerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("StopAccessLogging", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ContainerName"=>ContainerName), args)); aws_config=aws_config)

"""
    TagResource()

Adds tags to the specified AWS Elemental MediaStore container. Tags are key:value pairs that you can associate with AWS resources. For example, the tag key might be \"customer\" and the tag value might be \"companyA.\" You can specify one or more tags to add to each container. You can add up to 50 tags to each container. For more information about tagging, including naming and usage conventions, see Tagging Resources in MediaStore.

# Required Parameters
- `Resource`: The Amazon Resource Name (ARN) for the container.
- `Tags`: An array of key:value pairs that you want to add to the container. You need to specify only the tags that you want to add or update. For example, suppose a container already has two tags (customer:CompanyA and priority:High). You want to change the priority tag and also add a third tag (type:Contract). For TagResource, you specify the following tags: priority:Medium, type:Contract. The result is that your container has three tags: customer:CompanyA, priority:Medium, and type:Contract.

"""
tag_resource(Resource, Tags; aws_config::AWSConfig=global_aws_config()) = mediastore("TagResource", Dict{String, Any}("Resource"=>Resource, "Tags"=>Tags); aws_config=aws_config)
tag_resource(Resource, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Resource"=>Resource, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

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

# Required Parameters
- `Resource`: The Amazon Resource Name (ARN) for the container.
- `TagKeys`: A comma-separated list of keys for tags that you want to remove from the container. For example, if your container has two tags (customer:CompanyA and priority:High) and you want to remove one of the tags (priority:High), you specify the key for the tag that you want to remove (priority).

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


text:
function Bpenning{T}( x1::T,x2::T,x3::T )
  B=zeros(T,3)
  JB=zeros(T,3,3)

  B[1]=1

  return B,JB
end

x1=1.;x2=2.;x3=3.;

function skew{T}(v::Array{T,1} )
  return [0  -v[3]  v[2]; v[3] 0 -v[1]; -v[2]  v[1]  0 ]
end

function odeHB{T}( x1::T,x2::T,x3::T,v1::T,v2::T,v3::T, q::T, m::T )

  B,JB=Bpenning(x1,x2,x3)
  F=[v1; v2; v3; cross([v1; v2; v3],B) ]
  DF=[zeros(T,3,3)  eye(T,3,3);skew([v1; v2; v3])*JB   -skew(B)]

 return F,DF
end

F,DF=odeHB(1.0,2.,3.,4.,5.,6.)

expm(DF*dt)*[x1;x2;x3;v1;v2;v3]

pinv(DF)*expm(DF)

deltat=0.1

cross(rand(3),rand(3,3))

ext: jl

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
#---------------------------------------------------------------------#
#This code computes the global matrices using DSS
#Written by F.X. Giraldo on April 24, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function global_matrices!(M,L,ψ,dψ,ξ_x,ξ_y,ξ_z,η_x,η_y,η_z,ζ_x,ζ_y,ζ_z,jac,ωq,intma,Ne,Np,Nq,Npoin,DFloat)

    #Initialize Matrices
    inode=zeros(Int64,Np,Np,Np)

    @inbounds for e=1:Ne

        #Store Coordinates
        for k=1:Np, j=1:Np, i=1:Np
            I=intma[i,j,k,e]
            inode[i,j,k]=I
        end #i,j,k

        #Do LGL Integration
        for n=1:Nq, m=1:Nq, l=1:Nq
            wq=ωq[l]*ωq[m]*ωq[n]*jac[l,m,n,e]
            ξ1, ξ2, ξ3 = ξ_x[l,m,n,e], ξ_y[l,m,n,e], ξ_z[l,m,n,e]
            η1, η2, η3 = η_x[l,m,n,e], η_y[l,m,n,e], η_z[l,m,n,e]
            ζ1, ζ2, ζ3 = ζ_x[l,m,n,e], ζ_y[l,m,n,e], ζ_z[l,m,n,e]

            #Loop through I points
            for k=1:Np, j=1:Np, i=1:Np
                I=inode[i,j,k]
                h_i=ψ[i,l]*ψ[j,m]*ψ[k,n]
                h_ξ=dψ[i,l]*ψ[j,m]*ψ[k,n]
                h_η=ψ[i,l]*dψ[j,m]*ψ[k,n]
                h_ζ=ψ[i,l]*ψ[j,m]*dψ[k,n]
                δhδx_i=h_ξ*ξ1 + h_η*η1 + h_ζ*ζ1
                δhδy_i=h_ξ*ξ2 + h_η*η2 + h_ζ*ζ2
                δhδz_i=h_ξ*ξ3 + h_η*η3 + h_ζ*ζ3

                #Loop through J points
                for kj=1:Np, jj=1:Np, ij=1:Np
                    J=inode[ij,jj,kj]
                    h_j=ψ[ij,l]*ψ[jj,m]*ψ[kj,n]
                    h_ξ=dψ[ij,l]*ψ[jj,m]*ψ[kj,n]
                    h_η=ψ[ij,l]*dψ[jj,m]*ψ[kj,n]
                    h_ζ=ψ[ij,l]*ψ[jj,m]*dψ[kj,n]
                    δhδx_j=h_ξ*ξ1 + h_η*η1 + h_ζ*ζ1
                    δhδy_j=h_ξ*ξ2 + h_η*η2 + h_ζ*ζ2
                    δhδz_j=h_ξ*ξ3 + h_η*η3 + h_ζ*ζ3
                    M[I,J] += wq*h_i*h_j
                    L[I,J] -= wq*( δhδx_i*δhδx_j + δhδy_i*δhδy_j + δhδz_i*δhδz_j )
                end #ij, jj, kj
            end #i, j, k
        end #l, m, n
    end #e

    return (M,L)
end


text:
###############################################################################
#
# Reduce step method to find peak for multiple variables
#
###############################################################################
function find_peak(
            f::F,
            ms::ReduceStepMethodND{FT},
            tol::SolutionToleranceND{FT}
) where {F<:Function, FT<:AbstractFloat}
    # define the initial step
    @unpack x_inis, x_maxs, x_mins, x_targ, x_temp, Δ_inis, Δ_oper, Δjd = ms;

    # initial the value
    tar_y = f(x_targ);

    # find the solution
    Nxs = length(x_targ);
    while true
        # total count of updates for all the x values
        count_all = 0;

        for ith in 1:Nxs
            # 1. increase the x by Δx till tar_y is bigger
            count_inc = 0
            while true
                x_temp[ith] = x_targ[ith] + Δ_oper[ith];
                x_temp[ith] > x_maxs[ith] ? break : nothing;
                new_y = f(x_temp);
                new_y >  tar_y ? (x_targ .= x_temp; tar_y=new_y) : break;
                count_inc += 1;
                count_all += 1;
            end
            x_temp .= x_targ;

            # 2. decrease the x by Δx till tar_y is bigger
            while count_inc==0
                x_temp[ith] = x_targ[ith] - Δ_oper[ith];
                x_temp[ith] < x_mins[ith] ? break : nothing;
                new_y = f(x_temp);
                new_y >= tar_y ? (x_targ .= x_temp; tar_y=new_y) : break;
                count_all += 1;
            end
            x_temp .= x_targ;
        end

        # 3. if no update, then 10% the Δx
        if count_all==0
            for ith in 1:Nxs
                Δ_oper[ith] > tol.tol[ith] ?
                              (Δ_oper[ith] /= 10) : (Δjd[ith] = true);
            end
        end

        # 4. judge whether to break
        all( Δjd ) ? break : nothing;
    end

    return x_targ
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# contributed by Jarret Revels and Alex Arslan
# based on the Java version

module NBody

using Printf
using LinearAlgebra

# Constants
const solar_mass = 4 * pi * pi
const days_per_year = 365.24

struct Vec3
    x::NTuple{3, Float64}
end
Vec3(x, y, z) = Vec3((x,y,z))
Base.:/(v::Vec3, n::Number) = Vec3(1/n .* v.x)
Base.:*(v::Vec3, n::Number) = Vec3(n .* v.x)
Base.:-(v1::Vec3, v2::Vec3) = Vec3(v1.x .- v2.x)
Base.:+(v1::Vec3, v2::Vec3) = Vec3(v1.x .+ v2.x)
# Todo, prettify
squarednorm(v1::Vec3) = v1.x[1]^2 + v1.x[2]^2 + v1.x[3]^2

# A heavenly body in the system
mutable struct Body
    pos::Vec3
    v::Vec3
    mass::Float64
end

function offset_momentum!(b::Body, p::Vec3)
    b.v -= p / solar_mass
end

function init_sun!(bodies::Vector{Body})
    p = Vec3((0.0, 0.0, 0.0))
    for b in bodies
        p += b.v * b.mass
    end
    offset_momentum!(bodies[1], p)
end

function advance(bodies::Vector{Body}, dt::Number)
    for i = 1:length(bodies)
        @inbounds for j = i+1:length(bodies)
            delta = bodies[i].pos - bodies[j].pos
            dsq = squarednorm(delta)
            distance = sqrt(dsq)
            mag = dt / (dsq * distance)

            bodies[i].v -= delta * (bodies[j].mass * mag)
            bodies[j].v += delta * (bodies[i].mass * mag)
        end
    end

    for b in bodies
        b.pos += b.v * dt
    end
end

function energy(bodies::Vector{Body})
    e = 0.0
    for i = 1:length(bodies)
        e += 0.5 * bodies[i].mass * squarednorm(bodies[i].v)
        @inbounds for j = i+1:length(bodies)
            delta = bodies[i].pos - bodies[j].pos
            distance = sqrt(squarednorm(delta))
            e -= (bodies[i].mass * bodies[j].mass) / distance
        end
    end
    return e
end

function perf_nbody(N::Int=1000)
    jupiter = Body( Vec3(4.84143144246472090e+00,                  # pos[1] = x
                         -1.16032004402742839e+00,                 # pos[2] = y
                         -1.03622044471123109e-01),                # pos[3] = z
                   Vec3(1.66007664274403694e-03 * days_per_year,   # v[1] = vx
                        7.69901118419740425e-03 * days_per_year,   # v[2] = vy
                        -6.90460016972063023e-05 * days_per_year), # v[3] = vz
                   9.54791938424326609e-04 * solar_mass)       # mass

    saturn = Body(Vec3(8.34336671824457987e+00,
                       4.12479856412430479e+00,
                      -4.03523417114321381e-01),
                  Vec3(-2.76742510726862411e-03 * days_per_year,
                       4.99852801234917238e-03 * days_per_year,
                       2.30417297573763929e-05 * days_per_year),
                  2.85885980666130812e-04 * solar_mass)

    uranus = Body(Vec3(1.28943695621391310e+01,
                    -1.51111514016986312e+01,
                    -2.23307578892655734e-01),
                  Vec3(2.96460137564761618e-03 * days_per_year,
                       2.37847173959480950e-03 * days_per_year,
                       -2.96589568540237556e-05 * days_per_year),
                  4.36624404335156298e-05 * solar_mass)

    neptune = Body(Vec3(1.53796971148509165e+01,
                        -2.59193146099879641e+01,
                        1.79258772950371181e-01),
                   Vec3(2.68067772490389322e-03 * days_per_year,
                        1.62824170038242295e-03 * days_per_year,
                       -9.51592254519715870e-05 * days_per_year),
                   5.15138902046611451e-05 * solar_mass)

    sun = Body(Vec3(0.0, 0.0, 0.0), Vec3(0.0, 0.0, 0.0), solar_mass)

    bodies = [sun, jupiter, saturn, uranus, neptune]

    init_sun!(bodies)
    @printf("%.9f\n", energy(bodies))
    for i = 1:N
        advance(bodies, 0.01)
    end
    @printf("%.9f\n", energy(bodies))
end

end # module

n = parse(Int,ARGS[1])
NBody.perf_nbody(n)


text:
__precompile__(true)

module BitStreamLib

import Base.push!
import Base.length
import Base.get
import Base.getindex
import Base.setindex!
import Base.string
import Base.show
import Base.iterate

export BitStream
export push!
export get
export getindex
export set
export setindex!
export length
export string
export show
export iterate
export data
export getvalue
export parity
export truecount
export falsecount

struct BitStream
	data::Array{UInt8, 1}
	indices::Array{Int, 1}
	BitStream() = new(zeros(UInt8, 16), [1, 1])
	BitStream(data::Array{UInt8, 1}, len::Int) = new(data, [fld(len, 8) + 1, (len) % 8 + 1])
	BitStream(data::Array{UInt8, 1}) = BitStream(data, length(data) * 8)
	BitStream(data::String) = BitStream(Array{UInt8, 1}(data))
	BitStream(data::Union{Signed, Unsigned, AbstractFloat}) = BitStream(Array{UInt8, 1}(reinterpret(UInt8, [data])))
end

function push!(strm::BitStream, data::Bool)
	strm.data[strm.indices[1]] |= UInt8(UInt8(data ? UInt8(1) : UInt8(0)) << (strm.indices[2] - 1))
	strm.indices[2] += 1

	if strm.indices[2] == 9
		strm.indices[2] = 1
		strm.indices[1] += 1
		if strm.indices[1] > length(strm.data)
			push!(strm.data, UInt8(0))
		end
	end
end

getindex(strm::BitStream, index::Int)::Bool = get(strm, index)

getindex(strm::BitStream, range::UnitRange{Int64}) = [strm[i] for i in range]

function get(strm::BitStream, index::Int)::Bool
	return (strm.data[cld(index, 8)] & (UInt8(1) << ((index - 1) % 8))) != 0
end

setindex!(A::BitStream, val::Bool, i1::Int) = set(A, i1, val)

function setindex!(strm::BitStream, vals::Array{Bool,1}, range::UnitRange{Int64})
	i = 1
	for j in range
		strm[j] = vals[i]
		i += 1
	end
end

function set(strm::BitStream, index::Int, val::Bool)
	strm.data[cld(index, 8)] &= ~UInt8(UInt8(1) << ((index - 1) % 8))
	strm.data[cld(index, 8)] |= UInt8((val ? UInt8(1) : UInt8(0)) << ((index - 1) % 8))
end

function length(strm::BitStream)::Int
	return (strm.indices[1] - 1) * 8 + strm.indices[2] - 1
end

function data(strm::BitStream)::Array{UInt8, 1}
	return strm.data[1:cld(length(strm), 8)]
end

function getvalue(T::DataType, strm::BitStream)
	if T <: Union{Signed, Unsigned, AbstractFloat}
		return reinterpret(T, data(strm))[1]
	elseif T <: AbstractString
		return String(data(strm))
	end
end

function string(strm::BitStream)
	ret = ""
	for i in 1 : length(strm)
		ret *= string(strm[i] ? '1' : 0)
	end
	return ret
end

function show(io::IO, strm::BitStream)
	print(io, string(strm))
end

function iterate(iter::BitStream, state=1)
	if state > length(iter)
		return nothing
	end

	return (iter[state], state + 1)
end

parity(strm::BitStream)::Bool = reduce(⊻, strm)

truecount(strm::BitStream)::Int = mapreduce(x -> x ? 1 : 0, +, strm)

falsecount(strm::BitStream)::Int = mapreduce(x -> x ? 0 : 1, +, strm)

end


text:
"""

#### References

Enhanced random search based incremental extreme learning machine.

Guang-Bin Huang and Lei Chen.

Neurocomputing, 2008 vol. 71 (16-18) pp. 3460-3468.

http://linkinghub.elsevier.com/retrieve/pii/S0925231207003633
"""
mutable struct EIELM{TA,TV} <: AbstractSLFN where TA<:AbstractActivation where TV<:AbstractArray{Float64}
    p::Int              # Number of training points
    q::Int              # Dimensionality of function domain
    Lmax::Int           # maximum number of neurons
    k::Int              # maximum number of trials in node addition
    ϵ::Float64          # expected learning accuracy
    activation::TA
    μx::Vector{Float64}
    σx::Vector{Float64}
    At::Matrix{Float64} # transpose of W matrix
    b::Vector{Float64}
    v::TV

    function EIELM(p::Int, q::Int, Lmax::Int, k::Int, ϵ::Float64, activation::TA,
                   μx, σx) where TA
        WARNINGS[1] && @warn("This is experimental and doesn't work properly")
        At = Array{Float64}(undef, q, Lmax)  # will chop later
        b = Array{Float64}(undef, Lmax)      # ditto
        v = Array{Float64}(undef, Lmax)      # ditto
        new{TA,typeof(v)}(p, q, Lmax, k, ϵ, activation, μx, σx, At, b, v)
    end
end

function EIELM(x::AbstractArray, y::AbstractActivation;
               activation::AbstractActivation=Tanh(),
               Lmax::Int=size(x, 1), k::Int=20,
               ϵ::Float64=1e-6)
    q = size(x, 2)  # dimensionality of function domain
    p = size(x, 1)  # number of training points
    Lmax = min(Lmax, p)
    xn, μx, σx = normalize(x[:, :])
    out = EIELM(p, q, Lmax, k, ϵ, activation, μx, σx)
    fit!(out, xn, y)
end

## API methods
function fit!(elm::EIELM, x::AbstractArray, y::AbstractArray)
    # Step 1: Initialization
    L = 0
    E = copy(y)
    min_Ei = copy(E)

    # Step 2: learning
    while L < elm.Lmax && norm(E) > elm.ϵ
        L += 1
        min_Ei_contrib = Inf

        for i in 1:elm.k
            Ai = 2*rand(elm.q, 1) - 1  # uniform [-1, 1]
            bi = 2*rand() - 1          # uniform [-1, 1]
            Hi = hidden_out(elm, x, Ai, bi)
            vi = dot(E, Hi) / dot(Hi, Hi)
            Ei = E - vi*Hi
            Ei_contrib = norm(Ei)

            if Ei_contrib < min_Ei_contrib
                elm.At[:, L] = Ai
                elm.b[L] = bi
                elm.v[L] = vi
                min_Ei_contrib = Ei_contrib
                copy!(min_Ei, Ei)
            end

        end

        copy!(E, min_Ei)
    end

    # chop before returning
    elm.At = elm.At[:, 1:L]
    elm.b = elm.b[1:L]
    elm.v = elm.v[1:L]
    elm
end

function Base.show(io::IO, elm::EIELM{TA}) where TA
    s =
    """
    EIELM with
      - $(TA) Activation function
      - $(elm.q) input dimension(s)
      - $(length(elm.b)) neuron(s)
      - $(elm.p) training point(s)
    """
    print(io, s)
end


text:
function mt_wavedirection(signal::AbstractMatrix{T}, n, n_overlap; fs=1,onesided=true) where {T}
    config = DSP.MTSpectrogramConfig{T}(size(signal,2), n, n_overlap; fs=fs,onesided=onesided)
    configcross = DSP.MTCrossSpectraConfig{T}(size(signal,1), n; fs = fs,onesided=onesided)
    X = DSP.allocate_output(config)
    tempoutput = DSP.allocate_output(configcross)
    epochs = [DSP.arraysplit(r,n,n_overlap) for r in eachrow(signal)]
    shortsignal = zeros(eltype(epochs[1][1]), size(signal,1), length(epochs[1][1]))
    @views for (time_index,(z,x,y)) in enumerate(zip(epochs...))
        shortsignal[1,:] .= z
        shortsignal[2,:] .= x
        shortsignal[3,:] .= y
        mt_cross_power_spectra!(tempoutput, shortsignal, configcross)
        X[:,time_index] .= angle.(imag.(tempoutput[2,1,:]) .+ 1im .* imag.(tempoutput[3,1,:]))
    end
    return DSP.Periodograms.Spectrogram(X, config.mt_config.freq, config.time)
end

function compute_spectrogram(waveseries, wavetime)
    continuouswave = reduce(vcat, w[:,1] for w in waveseries)
    spec = mt_spectrogram(continuouswave, 2304, 2304÷2, fs = 1) # fs has to be int for some reason
    spectrogram = spec.power./(2π * 1.28) ./2 # equivalent to * Δ/2π and div 2 because onesided
    spectime = wavetime[1]:(wavetime[2]-wavetime[1])/2:wavetime[end]
    Ω = collect((2π*1.28) .* spec.freq)
    return spectrogram, Ω, spectime
end

function compute_errorgram(waveseries, wavetime)
    continuouswave = [reduce(vcat, w[:,i] for w in waveseries) for i in 1:3]
    spec = [mt_spectrogram(c, 2304, 0) for c in continuouswave]
    spectrogram = [s.power./(2π * 1.28) ./2 for s in spec] # equivalent to * Δ/2π and div 2 because onesided (DSP only allows spectrogram with integer fs so rescale is necessary)
    Ω = collect((2π*1.28) .* spec[1].freq)
    errorgram = @. log(spectrogram[2]+spectrogram[3])+2*log(tanh(OceanWaveSpectralFitting.approx_dispersion(Ω,40)))-log(spectrogram[1])
    return errorgram, Ω, wavetime
end


text:
"""
    AutoBody(sdf,map=(x,t)->x) <: AbstractBody

    - sdf(x::AbstractVector,t::Real)::Real: signed distance function
    - map(x::AbstractVector,t::Real)::AbstractVector: coordinate mapping function

Define a geometry by its `sdf` and optional coordinate `map`. All other
properties are determined using Automatic Differentiation. Note: the `map`
is composed automatically if provided, ie `sdf(x,t) = sdf(map(x,t),t)`.
"""
struct AutoBody{F1<:Function,F2<:Function} <: AbstractBody
    sdf::F1
    map::F2
    function AutoBody(sdf,map=(x,t)->x)
        comp(x,t) = sdf(map(x,t),t)
        new{typeof(comp),typeof(map)}(comp, map)
    end
end

using ForwardDiff,DiffResults
using LinearAlgebra: norm2
"""
    measure(body::AutoBody,x::Vector,t::Real)

ForwardDiff is used to determine the geometric properties from the `sdf`.
Note: The velocity is determined _soley_ from the optional `map` function.
"""
function measure(body::AutoBody,x::AbstractVector,t::Real)
    # V = dot(map), ignoring any other time dependancy in sdf.
    V = -ForwardDiff.derivative(t->body.map(x,t), t)

    # Use DiffResults to get Hessian/gradient/value in one shot
    result = DiffResults.HessianResult(x)
    result = ForwardDiff.hessian!(result, x -> body.sdf(x,t), x)
    d = Float64(DiffResults.value(result))
    n̂ = DiffResults.gradient(result)
    H,K = curvature(DiffResults.hessian(result))

    # |∇d|=1 for a true distance function. If sdf or map violates this condition,
    # scaling by |∇d| gives an approximate correction.
    m = norm2(n̂)
    d/m,n̂./m,V,H/m,K/m^2
end

using LinearAlgebra: tr
"""
    curvature(A::AbstractMatrix)

Return `H,K` the mean and Gaussian curvature from `A=hessian(sdf)`.
K=tr(minor(A)) in 3D and K=0 in 2D.
"""
function curvature(A::AbstractMatrix)
    H,K = 0.5*tr(A),0
    if size(A)==(3,3)
        K = A[1,1]*A[2,2]+A[1,1]*A[3,3]+A[2,2]*A[3,3]-A[1,2]^2-A[1,3]^2-A[2,3]^2
    end
    H,K
end


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

## embedding stuff taken from my IntrinsicGeometry.jl package

function reduce_to_principal_components(data::AbstractMatrix{<:Real}, n_dims::Int)
  eigval, eigvec = data' |> cov |> Hermitian |> eigen
  means = mean(data, dims=2)
  n_dims = min(n_dims, size(data, 1))
  projection = eigvec[:, (end-n_dims+1):end]'
  reduced = projection * (data .- means)
  eigenvalues = eigval[(end-n_dims+1):end]
  (; reduced, projection, means, eigenvalues)
end

function compute_potential_distances(transitions::AbstractMatrix{T},
    gamma::T = T(1))::Matrix{T} where T
  n_points = size(transitions, 1)
  @assert size(transitions, 1) == size(transitions, 2)
  if gamma == 1
    f = (p::T) -> log(max(p, floatmin(T)))
  elseif gamma == -1
    f = (p::T) -> p
  else
    f = (p::T) -> 2/(1-gamma) * p^((1-gamma)/2)
  end
  distances = zeros(T, n_points, n_points)
  log_transitions = f.(transitions)
  for k in 1:n_points
    Threads.@threads for i in 1:n_points
      for j in 1:(i-1)
        difference = log_transitions[i, k] - log_transitions[j, k]
        distances[j, i] += difference ^ 2
      end
    end
  end
  for i in 1:n_points, j in 1:(i-1)
    distances[i, j] = distances[j, i]
  end
  distances .|> sqrt
end

using PyCall, Conda # maybe this shouldn't be necessary to import?

function perform_classic_mds(distances::AbstractMatrix, n_manifold_dims::Int)
  squared = distances .^ 2
  semi_centered = distances .- mean(distances, dims=1)
  centered = semi_centered .- mean(semi_centered, dims=2)
  embedding = reduce_to_principal_components(centered, n_manifold_dims)
  embedding.reduced
end

function perform_metric_mds(distances::AbstractMatrix{T},
    n_manifold_dims::Int)::AbstractMatrix{T} where T
  classic_result = perform_classic_mds(distances, n_manifold_dims)
  skl_manifold = pyimport_conda("sklearn.manifold", "scikit-learn")
  embedding, stress, n_iter = skl_manifold.smacof(distances, metric=true,
    n_components=n_manifold_dims, return_n_iter=true, n_init=1, init=classic_result')
  convert.(T, embedding)
end

function fit_line(x::AbstractVector{T}, y::AbstractVector{T}) where T <: Real
  x_bar, y_bar = mean(x), mean(y)
  numerator = sum((x .- x_bar) .* (y .- y_bar))
  denominator = sum((x .- x_bar).^2)
  slope = numerator / denominator # how to represent an infinnite slope? don't: filter for NaNs later
  intercept = y_bar - slope*x_bar
  fit = intercept .+ slope.*x
  error = sum((fit .- y).^2)
  (; slope, intercept, error)
end

# very simple method, which I believe is the same one harnessed by and for PHATE
function find_knee(x::AbstractVector{T}, y::AbstractVector{T})::T where T <: Real
  @assert length(x) == length(y)
  n_points = length(x)
  errors = zeros(T, n_points)
  for (i, pivot) in enumerate(x)
    # do not exclude the pivot for symmetry. we need to sum over the same set of items each time for commensurable error tallies
    left_mask = (x .<= pivot)
    right_mask = .!left_mask
    left_fit = @views fit_line(x[left_mask], y[left_mask])
    right_fit = @views fit_line(x[right_mask], y[right_mask])
    total_error = left_fit.error + right_fit.error
    @debug pivot left_fit right_fit
    errors[i] = total_error
  end
  masked_errors = ifelse.(isfinite.(errors), errors, Inf)
  # possibly disallow lines of length 1 and 2 as there are 2 degrees of freedom
  best_index = argmin(masked_errors[2:end-1]) + 1
  x[best_index]
end

function estimate_von_neumann_entropy(eigval::AbstractVector{T})::T where T <: Real
  positive_eigval = eigval[eigval .> 0]
  normalizer = sum(positive_eigval)
  probabilities = positive_eigval / normalizer
  pointwise_entropies = -probabilities .* log.(probabilities)
  sum(pointwise_entropies)
end

function optimize_von_neumann_entropy(eigval::Vector{T},
    powers::AbstractRange{T})::T where T # <:Real is implied
  n_powers = length(powers)
  x, y = zeros(T, n_powers), zeros(T, n_powers)
  for (i, t) in enumerate(powers)
    step_eigval = eigval.^t
    entropy = estimate_von_neumann_entropy(step_eigval)
    x[i] = t
    y[i] = entropy
  end
  knee = find_knee(x, y)
  #convert(Int, knee) # no InexactErrors here
end

function decompose_diffusion(diffusion::Matrix{T})::Tuple{Vector{T}, Matrix{T}, Matrix{T}} where T
  scales = sum(diffusion, dims=2).^0.5 # symmetric version
  transitions = diffusion ./ (scales .* scales')
  eigval, eigvec = transitions |> Hermitian |> eigen
  eigval, eigvec, scales
end

function embed_diffusion(diffusion::Matrix{T}, n_manifold_dims::Int;
    n_steps::T = T(-1), step_powers::AbstractRange{T} = 1:T(50)
    )::Matrix{T} where T
  eigval, eigvec, scales = decompose_diffusion(diffusion)
  transition = diffusion ./ (scales.^2)
  if n_steps <= 0
    n_steps = optimize_von_neumann_entropy(eigval, step_powers)
    @info "Taking $n_steps diffusion steps."
  end
  distances = compute_potential_distances(
    transition ^ n_steps |> real)
  perform_metric_mds(distances, n_manifold_dims)
end


text:
using Mustache

# TODO: currently auto and log aren't implemented properly
_axis_type(s::Symbol) = s == :auto ? (:BasicTicker, :BasicTickFormatter) :
                        s == :linear ? (:BasicTicker, :BasicTickFormatter) :
                        s == :log ? (:BasicTicker, :BasicTickFormatter) :
                        s == :datetime ? (:DatetimeTicker, :DatetimeTickFormatter) :
                        error("Unknown axis type")

function _genmodels(p::Plot)
    bkplot = Bokehjs.Plot()
    doc = Bokehjs.uuid4()
    obs = Dict{AbstractString, BkAny}[]

    cdss = Bokehjs.ColumnDataSource[]
    renderers = Bokehjs.PlotObject[]
    legends = Tuple[]
    for datacolumn in p.datacolumns
        cds = Bokehjs.ColumnDataSource(datacolumn.data)
        push!(cdss, cds)
        pushdict!(obs, cds, doc)

        glyph = datacolumn.glyph
        pushdict!(obs, glyph, doc)

        glyphrenderer = Bokehjs.GlyphRenderer(cds, glyph, nothing, glyph)
        push!(renderers, glyphrenderer)
        pushdict!(obs, glyphrenderer, doc)

        if datacolumn.legend != nothing
            push!(legends, (datacolumn.legend, glyphrenderer))
        end
    end

    cdss, renderers, legends

    if length(legends) > 0
        legend = Bokehjs.Legend(bkplot, legends, p.legendsgo)
        push!(renderers, legend)
        pushdict!(obs, legend, doc)
    end

    dr1x = Bokehjs.DataRange1d(cdss, AbstractString["x"])
    dr1y = Bokehjs.DataRange1d(cdss, AbstractString["y"])
    pushdict!(obs, dr1x, doc)
    pushdict!(obs, dr1y, doc)

    axes = Dict{Symbol, Vector}(:above => [],
                                :below => Bokehjs.LinearAxis[],
                                :left => Bokehjs.LinearAxis[],
                                :right => []
                               )

    for (i, loc, sym) in [(0, :below, :x), (1, :left, :y)]
        @eval begin
            tick_type, tickform_type = _axis_type(getfield($p, Symbol($"$(sym)_axis_type")))

            tick = Bokehjs.Ticker(tick_type)
            tickform = Bokehjs.TickFormatter(tickform_type)
            axis = Bokehjs.LinearAxis($i, tickform, tick, $bkplot)
            grid = Bokehjs.Grid($i, $bkplot, tick)

            for obj in (tick, axis, tickform, grid)
                pushdict!($obs, obj, $doc)
            end

            push!($renderers, axis)
            push!($renderers, grid)
            $axes[Symbol($"$loc")] = [axis]
        end
    end

    tools = Bokehjs.PlotObject[]
    if in(:pan, p.tools)
        pantool = Bokehjs.Metatool("PanTool", bkplot,
                                   AbstractString["width", "height"])
        pushdict!(obs, pantool, doc)
        push!(tools, pantool)
    end
    if in(:wheelzoom, p.tools)
        wheelzoomtool = Bokehjs.Metatool("WheelZoomTool", bkplot,
                                         AbstractString["width", "height"])
        pushdict!(obs, wheelzoomtool, doc)
        push!(tools, wheelzoomtool)
    end
    if in(:boxzoom, p.tools)
        boxzoomtool = Bokehjs.Metatool("BoxZoomTool", bkplot)
        pushdict!(obs, boxzoomtool, doc)
        push!(tools, boxzoomtool)
    end
    if in(:resize, p.tools)
        resizetool = Bokehjs.Metatool("ResizeTool", bkplot)
        pushdict!(obs, resizetool, doc)
        push!(tools, resizetool)
    end
    if in(:reset, p.tools)
        resettool = Bokehjs.Metatool("ResetTool", bkplot)
        pushdict!(obs, resettool, doc)
        push!(tools, resettool)
    end

    bkplot = Bokehjs.Plot(bkplot, dr1x, dr1y, renderers, axes, tools, p.title,
                          p.height, p.width)
    pushdict!(obs, bkplot, doc)

    plotcontext = Bokehjs.PlotContext(bkplot)
    pushdict!(obs, plotcontext, doc)

    indent = DEBUG ? 2 : 0
    json(obs, indent), plotcontext
end

function _obdict(ob::Bokehjs.PlotObject, doc::Bokehjs.UUID)
    d = Dict{AbstractString, BkAny}()
    d["id"] = ob.uuid
    extra_attrs = fieldnames(ob)
    d["type"] = in(:_type_name, extra_attrs) ? ob._type_name : typeof(ob)

    attrs = Dict{AbstractString, Any}()
    attrs["id"] = d["id"]
    attrs["doc"] = doc

    special = [:_type_name]

    for name in extra_attrs[2:end]
        in(name, special) && continue
        getfield(ob, name) == Bokehjs.omit && continue
        key = string(name)
        # key = begingswith(key, "_") ? key[2:end] : key
        attrs[key] = getfield(ob, name)
    end

    d["attributes"] = attrs

    return d
end

pushdict!(obs::Any, ob::Bokehjs.PlotObject, doc::Bokehjs.UUID) =
    push!(obs, _obdict(ob, doc))

_get_resources_dir() = Pkg.dir("Bokeh", "templates")

function _gettemplate(template::AbstractString,
                      path::Union{AbstractString, Void}=nothing)
    path = path == nothing ? _get_resources_dir() : path
    fname = joinpath(path, template)
    open(readstring, fname, "r")
end

function _bokehjs_paths(minified::Bool=true)
    dir = joinpath(dirname(@__FILE__), "..", "deps", "bokehjs")
    jspath = joinpath(dir, minified ? "bokeh.min.js" : "bokeh.js")
    csspath = joinpath(dir, minified ? "bokeh.min.css" : "bokeh.css")
    (jspath, csspath)
end

function _render_jscss(jspath::AbstractString, csspath::AbstractString,
                       buildin::Bool)
    if buildin
        css = open(readall, csspath, "r")
        js = open(readall, jspath, "r")
        return "<style>\n$css\n</style>\n<script type=\"text/javascript\">\n$js\n</script>\n"
    else
        return "<link rel=\"stylesheet\" href=\"$csspath\" type=\"text/css\" />\n"*
               "<script type=\"text/javascript\" src=\"$jspath\"></script>\n"
    end
end

function _rendertemplate(models::AbstractString, plotcon::Bokehjs.PlotContext,
                         isijulia::Bool)
    base = isijulia ? _gettemplate("ijulia.html") :
                      _gettemplate("standalone.html")

    main = _gettemplate("main.html")
    body = _gettemplate("body.html")

    jspath, csspath = _bokehjs_paths(!DEBUG)
    builtin = isijulia || INCLUDE_JS
    jscss = _render_jscss(jspath, csspath, builtin)

    if DEBUG
        open(replace(PLOTFILE, ".html", "") * ".json", "w") do f
            print(f, models)
        end
    end

    context = Dict{AbstractString, AbstractString}(
        "model_id" => string(plotcon.uuid),
        "all_models" => models,
        "div_id" => string(Bokehjs.uuid4())
    )

    main = render(main, context)
    body = render(body, context)

    maincontext = Dict{AbstractString, AbstractString}(
        "jscss" => jscss,
        "main" => main,
        "body" => body,
    )

    result = render(base, maincontext)
end

function renderplot(plot::Plot, isijulia::Bool=false)
    modelsjson, plotcontext = _genmodels(plot)
    _rendertemplate(modelsjson, plotcontext, isijulia)
end

renderplot() = renderplot(CURPLOT)

function genplot(p::Plot, filename::NullString=nothing)
    filename = filename == nothing ? p.filename : filename
    html = renderplot(p)
    if ispath(filename) && FILE_WARNINGS
        if WARN_FILE != filename
            println("$(filename) already exists, overwriting")
            global WARN_FILE = filename
        end
    end
    open(filename, "w") do f
        print(f, html)
    end
end

genplot() = genplot(CURPLOT)
genplot(filename::AbstractString) = genplot(CURPLOT, filename)


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

# Test approximate equality of vectors or columns of matrices modulo floating
# point roundoff and phase (sign) differences.
#
# This function is designed to test for equality between vectors of floating point
# numbers when the vectors are defined only up to a global phase or sign, such as
# normalized eigenvectors or singular vectors. The global phase is usually
# defined consistently, but may occasionally change due to small differences in
# floating point rounding noise or rounding modes, or through the use of
# different conventions in different algorithms. As a result, most tests checking
# such vectors have to detect and discard such overall phase differences.
#
# Inputs:
#     a, b:: StridedVecOrMat to be compared
#     err :: Default: m^3*(eps(S)+eps(T)), where m is the number of rows
#
# Raises an error if any columnwise vector norm exceeds err. Otherwise, returns
# nothing.
function test_approx_eq_modphase(a::StridedVecOrMat{S}, b::StridedVecOrMat{T},
                                 err = length(axes(a,1))^3*(eps(S)+eps(T))) where {S<:Real,T<:Real}
    @test axes(a,1) == axes(b,1) && axes(a,2) == axes(b,2)
    for i in axes(a,2)
        v1, v2 = a[:, i], b[:, i]
        @test min(abs(norm(v1-v2)),abs(norm(v1+v2))) ≈ 0.0 atol=err
    end
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using CGDycore

# Physical parameters
Param=CGDycore.PhysParameters();

# Grid
nx=60;
ny=2;
Param.Lx=300*1.e3;
Param.Ly=2*Param.Lx/nx;
x0=0;
y0=0;
Param.H=10*1.e3;
Param.OrdPoly=4;
#Horizontal grid size
dx=Param.Lx/nx/(Param.OrdPoly+1);
dz=min(dx/2,Param.H/10);
nz=ceil(Param.H/dz);
nz=40;

Boundary = (;WE="Period", BT="Period")
# Boundary.WE="Period";
# Boundary.BT="Period";
Param.hS="";
Param.Grid=CGDycore.CartGrid(nx,ny,Param.Lx,Param.Ly,x0,y0,CGDycore.OrientFaceCart,Boundary,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP[1]=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP[i]=Param.Grid.zP[i-1]+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z[i]=Param.Grid.z[i-1]+Param.Grid.dz;
end

# Model
Param.ProfRho="BaldaufCart";
Param.ProfRhoBGrd="BaldaufCart";
Param.ProfpBGrd="BaldaufCart";
Param.ProfTheta="BaldaufCart";
Param.ProfVel="Const";
Param.Damping=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.Source=false;
Param.xc=Param.Lx/3;
Param.d=5000;
Param.a=5*1.e3;
Param.T0=250;
Param.DeltaT=1.e-2;
Param.uMax=0;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation="Compressible";
Param.TopoS="";
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo="";

# Discretization
OrdPolyZ=1;
(CG,Param)=CGDycore.Discretization(Param.OrdPoly,OrdPolyZ,CGDycore.JacobiDG3,Param);
Param.HyperVisc=false;
Param.HyperDCurl=2.e5; #1.e14*(dx/LRef)^3.2;
Param.HyperDGrad=2.e5;
Param.HyperDDiv=2.2e5;
Param.Upwind=false;

# Output
Param.Flat=false;
Param.vtkFileName="BaldaufCart";
Param.vtk=0;
vtkGrid=CGDycore.vtkCGGrid(CG,CGDycore.TransCart,CGDycore.Topo,Param);
Param.cNames = Vector{String}(undef, 6)
Param.cNames[1]="Rho";
Param.cNames[2]="u";
Param.cNames[3]="w";
Param.cNames[4]="Th";
Param.cNames[5]="TPrime";
Param.cNames[6]="Pres";

# Initial conditions

# fRho
# fVel
# fTheta
# fTBGrd
# fT
# fpBGrd
# fRhoBGrd

U=zeros(CG.NumG,nz,Param.NumV);
U[:,:,Param.RhoPos]=CGDycore.Project(CGDycore.fRho,CG,Param);

(U[:,:,Param.uPos],U[:,:,Param.vPos])=CGDycore.ProjectVec(CGDycore.fVel,CG,Param);
U[:,:,Param.ThPos]=CGDycore.Project(CGDycore.fTheta,CG,Param).*U[:,:,Param.RhoPos];
Param.TBGrd=CGDycore.Project(CGDycore.fTBGrd,CG,Param);
Param.T=CGDycore.Project(CGDycore.fT,CG,Param);
pBGrd=CGDycore.Project(CGDycore.fpBGrd,CG,Param);
RhoBGrd=CGDycore.Project(CGDycore.fRhoBGrd,CG,Param);
OP=Param.OrdPoly+1;
NF=Param.Grid.NumFaces;
RhoBGrd=reshape(RhoBGrd[reshape(CG.Glob,OP*OP*NF,1),:],OP,OP,NF,nz);
Param.RhoBGrdF=0.5*(RhoBGrd[:,:,:,1:nz-1]+RhoBGrd[:,:,:,2:nz]);
Param.pBGrd=reshape(pBGrd[reshape(CG.Glob,OP*OP*NF,1),:],OP,OP,NF,nz);

# Integration
dtau=.40;
time=0;
#IntMethod="Rosenbrock"; #"RungeKutta";
IntMethod="RungeKutta";
Param.EndTime=1800;
nIter=Param.EndTime/dtau;
PrintTime=100;
PrintInt=PrintTime/dtau;
Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);

error("Success")

Param.RK=RungeKuttaMethod("RK4");
Param.ROS=RosenbrockMethod("ROSRK3");

switch IntMethod
  case "Rosenbrock"
    for i=1:nIter
      U=Rosenbrock(U,dtau,@FcnNHCurlVec,@Jac,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
      end
    end
  case "RungeKutta"
    for i=1:nIter
      i
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
      end
    end
end

end


text:
import LogDensityProblems
import TransformVariables

"""
    loglikelihood(transformed_gradient_problem, theta)

Return the value of the log likelihood function evaluated at `theta`.

# Arguments
- `transformed_gradient_problem`
- `theta`
"""
function loglikelihood(transformed_gradient_problem,
                       theta)
    transformed_log_density = parent(transformed_gradient_problem)
    transformation = transformed_log_density.transformation
    theta_inversetransformed = TransformVariables.inverse(transformation,
                                                          theta)
    result = _loglikelihood(transformed_gradient_problem,
                            theta_inversetransformed)
    return result
end

function _loglikelihood(transformed_gradient_problem,
                        theta)
    log_likelihood_value, gradient = LogDensityProblems.logdensity_and_gradient(transformed_gradient_problem,
                                                                                theta)
    return log_likelihood_value
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Cubature_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Cubature")
JLLWrappers.@generate_main_file("Cubature", UUID("7bc98958-0e37-5d67-a6ac-a3a19030071a"))
end  # module Cubature_jll


text:
using RDatasets

include("censored.jl")
include("cross_validate.jl")
include("initialize.jl")
include("precision_at_k.jl")
include("simple_glrms.jl")
include("fit_rdataset.jl")


text:
using FileIO
using IterTools
using Makie
using MakieRipserer
using MeshIO
using Ripserer
using SparseArrays
# import Pkg; Pkg.add(url="https://github.com/mtsch/TopologicalDatasets.jl")
using TopologicalDatasets

torus = sample(Torus(; r1=3), 1000)
res = ripserer(Alpha(torus); reps=true, dim_max=2)

# barcode
bcplt = plot_barcode(res; infinity=Observable(5))
save(joinpath(@__DIR__, "assets", "barcode.png"), bcplt)

# diagram
dgplt = plot(res; infinity=Observable(5))
save(joinpath(@__DIR__, "assets", "diagram.png"), dgplt)

# cocycle
repplt = plot(torus)
plot!(repplt, res[end][end - 1].representative, torus)
save(joinpath(@__DIR__, "assets", "cocycle.png"), repplt)

# app demo

time = Observable(0.0)
torus_app = MakieRipserer.app(torus; time=time)
rng = range(0, 4.5; length=200)
rng = vcat(rng, reverse(rng[2:(end - 1)]))
record(torus_app, joinpath(@__DIR__, "assets", "torus.gif"), rng) do r
    println("time: ", r)
    time[] = r
end

# sublevel cat
cat = load(joinpath(@__DIR__, "assets", "cat.obj"))
cat_pts = cat.position
cat_index = Dict([cat_pts[i] => i for i in 1:length(cat_pts)])
dist = zeros(length(cat_pts), length(cat_pts))
for t in cat
    for (p1, p2) in IterTools.subsets(t, Val(2))
        i = cat_index[p1]
        j = cat_index[p2]
        dist[i, j] = dist[j, i] = max(p1[3], p2[3])
    end
end
for point in cat_pts
    i = cat_index[point]
    dist[i, i] = point[3]
end
cat_flt = Rips(sparse(dist))
time = Observable(-1.0)
cat_app = MakieRipserer.app(cat_pts, cat_flt; time=time, dim_max=2)
rng = range(
    minimum(Ripserer.births(cat_flt)), maximum(Ripserer.births(cat_flt)); length=200
)
rng = vcat(rng, reverse(rng[2:(end - 1)]))
record(cat_app, joinpath(@__DIR__, "assets", "cat.gif"), rng) do r
    println("time: ", r)
    time[] = r
end


text:
## Factorization code

## Factor
##
##  1 0 -p1                   1 0 0 0 0       -p1
##  0 1 -p2  --> Yn + X -->   0 1 0 0 0       -p2
##  0 0 -p3                   0 0 1 0 -1  +   -p3
##                            0 0 0 1 0       - 1
## we pass in ps = [-p1, -p2, ..., -pn, -1] the whol column including 1
## We can leave as Ct * B * D or D * Ct * B (the default)
##
function R_factorization(ps::Vector{Complex{T}}, Ct, B, side=Val{:left}) where {T}

    N = length(ps) - 1 # ps has 1 in it?
    c,s,tmp = givensrot(conj(ps[N]), -one(Complex{T}))

    nrm = norm(c)
    alpha = c/nrm

    alpha = alpha

    Ct[N] = Rotator(conj(c), -s, N)
    B[N] = Rotator(-s*alpha, norm(c), N)

    ## do we leave on left? If so, it must pass through Cts and B
    gamma = side == Val{:left} ? alpha : one(Complex{T})

    @inbounds for i in (N-1):-1:1

        c,s,tmp = givensrot(ps[i], tmp)

        Ct[i] = Rotator(conj(c*gamma), -s, i)
        B[i] = Rotator(c * gamma, s, i)

    end

    side == Val{:left} ? alpha : conj(alpha)
end

## Factor R =
# 1 0 0 p1     1 0 0 0       p1
# 0 1 0 p2 --> 0 1 0 0  -->  p2
# 0 0 1 p3     0 0 0 -1      p3
#              0 0 1 0       -1
function R_factorization(ps::Vector{T}, Ct, B, side=Val{:not_applicable}) where {T <: Real}
    N = length(ps) - 1
    c, s, tmp = givensrot(ps[N], -one(T))

    Ct[N] = Rotator(c, -s, N)
    B[N] = Rotator(-s, c, N)

    @inbounds for i in (N-1):-1:1
        c,s,tmp = givensrot(ps[i], tmp)
        Ct[i] = Rotator(c, -s, i)
        B[i] = Rotator(c,s, i)
    end

    nothing
end

function Q_factorization(state::FactorizationType{T, St, P, Tw}) where {T, St, P, Tw}
    N = state.N
    Q = state.Q
    S = St == Val{:SingleShift} ? Complex{T} : T # type of cosine term in rotator
    zS, oS,zT,oT = zero(S), one(S), zero(T), one(T)
    @inbounds for ii = 1:(N-1)
        Q[ii] = Rotator(zS, oT, ii)

    end
    Q[N] = Rotator(oS, zT, N)
end

## Init State uses factorization of R and Q above
function init_state(state::FactorizationType{T, St, P, Val{:NotTwisted}}, decompose) where {T, St, P}
    # (Q,D, sigma)
    # Ct, B [Ct1, B1]
    init_triu(state, decompose)
    Q_factorization(state)
end

function init_state(state::FactorizationType{T, St, P, Val{:IsTwisted}}, decompose) where {T, St, P}
    # (Q,D, sigma)
    # Ct, B [Ct1, B1]
    alpha = init_triu(state, decompose)
    Q_factorization(state)
    ## XXX ... Need to twist the Q here .... XXXX

end

## Upper triangular piecies R and (V,W)
## populate (D,Ct,B)
function init_triu(state::FactorizationType{T, Val{:SingleShift}, Val{:NoPencil}, Tw}, decompose) where {T, Tw}
    ps = decompose(state.POLY)
    alpha =  R_factorization(ps, state.Ct, state.B)
    state.D[state.N] = alpha
    state.D[state.N+1] = conj(alpha)
end

function init_triu(state::FactorizationType{T, Val{:DoubleShift}, Val{:NoPencil}, Tw}, decompose) where {T, Tw}
    ps = decompose(state.POLY)
    alpha =  R_factorization(ps, state.Ct, state.B)
end

# Here we have V = D Ct B; W = D1 Ct1 B1
function init_triu(state::FactorizationType{T, Val{:SingleShift}, Val{:HasPencil}, Tw}, decompose) where {T, Tw}

    ps, qs = decompose(state.POLY)
    beta =  R_factorization(qs, state.Ct1, state.B1)
    state.D1[state.N] = beta
    state.D1[state.N+1] = conj(beta)

    alpha =  R_factorization(ps, state.Ct, state.B)
    state.D[state.N] = alpha
    state.D[state.N+1] = conj(alpha)
end

# Here we have V = D Ct B; W = D1 Ct1 B1
function init_triu(state::FactorizationType{T, Val{:DoubleShift}, Val{:HasPencil}, Tw}, decompose) where {T, Tw}
    ps, qs = decompose(state.POLY)
    beta =  R_factorization(qs, state.Ct1, state.B1)
    alpha =  R_factorization(ps, state.Ct, state.B)
end

# # If there is an issue, this function can be used to resetart the algorithm
# # could be merged with init_state?
# function restart(state::ShiftType{T}) where {T}
#     # try again
#     init_state(state)

#     for i in 1:state.N
#         state.REIGS[i] = state.IEIGS[i] = zero(T)
#     end
#     state.ctrs.zero_index = 0
#     state.ctrs.start_index = 1
#     state.ctrs.stop_index = state.N - 1
#     state.ctrs.it_count = 0
#     state.ctrs.tr = state.N - 2
# end


text:
using Test

using TestLib.Lens
using TestLib.Omega.Space

using Omega
using Omega.Space

@testset "Simple.memrand" begin
  ω = SimpleΩ()
  startIdx = Omega.Space.base(Vector{Int},0)
  nextIdx = increment

  @test test_resolve_repeatable(ω, startIdx, nextIdx)

  ω = SimpleΩ{Vector{Int}, Bool}()
  @test test_resolve_uniform(ω, startIdx, nextIdx)

  #TODO: Add tests for the other variants of memrand
end

idx(i::Int) = base(Vector{Int}, i)

@testset "(Un)linearize" begin
  ωs1 = SimpleΩ(); ωs1[idx(1)] = 0.535; ωs1[idx(2)] = 0.5210; ωs1[idx(3)] = 0.645;
  ωs2 = SimpleΩ(); ωs2[idx(1)] = 0.452; ωs2[idx(2)] = 0.6477; ωs2[idx(3)] = 0.732;
  ωs3 = SimpleΩ(); ωs3[idx(1)] = 0.535; ωs3[idx(2)] = 0.7607; ωs3[idx(3)] = 0.065;

  # unlinearize requires that sω and ωvec are the same length
  test_ωs_len_3   = [ωs1, ωs2, ωs3]
  test_vals_len_3 = [[0.123, 0.5, 0.7], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]

  @test test_lens_put_get(linearize, unlinearize, test_ωs_len_3, test_vals_len_3)

  @test test_lens_put_put(linearize, unlinearize, test_ωs_len_3, test_vals_len_3)
  @test test_lens_get_put(linearize, unlinearize, test_ωs_len_3, test_vals_len_3)
end


text:
module FLPDevelopment

using Reexport
@reexport using FLPprocess, DataFrames, CSV, Dates, CategoricalArrays
@reexport using StatsBase, HypothesisTests, Bootstrap, MixedModels, Random, Survival
@reexport using StatsPlots, StatsPlots.PlotMeasures, LaTeXStrings
import Statistics: median, std

include("groups.jl")
include("utilities.jl")
include("pokes_plotting.jl")
include("Analysis_fun.jl")
include("Plot_fun.jl")

export youngs, dario_youngs, females, sixty_days_old, age_exp, cnos_animals, wt_animals, below_thrs, VirusDict
export first_females_juveniles, second_females_juveniles
export first_females_group, second_females_group
export median, std
export frequency, fraction_true, bin_axis
export session_plot, maintitle!
export individual_summary, group_summary
export function_analysis, survivalrate_algorythm, cumulative_algorythm, hazardrate_algorythm, mediansurvival_analysis
export incorrect_fraction, median_ci_scatter
export add_bar!, add_info!, Difference

end # module


text:
# Affine decision function #
# ========================== #
mutable struct AffineDecisionConstraintBridge{T, S} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.ScalarAffineFunction{T}, S}
    decision_function::AffineDecisionFunction{T}
    set::S
end

function MOIB.Constraint.bridge_constraint(::Type{AffineDecisionConstraintBridge{T, S}},
                                           model,
                                           f::AffineDecisionFunction{T},
                                           set::S) where {T, S}
    # All decisions have been mapped to either the decision part constant or the variable part terms
    # at this point.
    F = MOI.ScalarAffineFunction{T}
    # Calculate total constant
    constant = f.variable_part.constant +
        f.known_part.constant
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.ScalarAffineFunction(f.variable_part.terms, zero(T)),
                                    MOIU.shift_constant(set, -constant))
    # Save the constraint index, the decision function, and the set, to allow modifications
    return AffineDecisionConstraintBridge{T, S}(constraint, f, set)
end

function MOI.supports_constraint(::Type{<:AffineDecisionConstraintBridge{T}},
                                 ::Type{<:AffineDecisionFunction{T}},
                                 ::Type{<:MOI.AbstractScalarSet}) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:AffineDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:AffineDecisionConstraintBridge{T, S}}) where {T, S}
    return [(MOI.ScalarAffineFunction{T}, S)]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:AffineDecisionConstraintBridge{T}},
                                              ::Type{<:AffineDecisionFunction{T}},
                                              S::Type{<:MOI.AbstractScalarSet}) where T
    return AffineDecisionConstraintBridge{T, S}
end

MOI.get(b::AffineDecisionConstraintBridge{T, S}, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, S}) where {T, S} = 1
MOI.get(b::AffineDecisionConstraintBridge{T, S}, ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, S}) where {T, S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::AffineDecisionConstraintBridge{T}) where T
    f = bridge.decision_function
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v, f.decision_part.terms)
    end
    g = AffineDecisionFunction(
        MOIU.filter_variables(v -> !from_decision(v), f.variable_part),
        copy(f.decision_part),
        copy(f.known_part))
    return g
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::AffineDecisionConstraintBridge{T}) where T
    return bridge.set
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::AffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::AffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.delete(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::AffineDecisionConstraintBridge{T,S},
                 f::AffineDecisionFunction{T}) where {T,S}
    # Update bridge functions and function constant
    bridge.decision_function = f
    # Change the function of the bridged constraints
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.ScalarAffineFunction(f.variable_part.terms, zero(T)))
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::AffineDecisionConstraintBridge{T,S}, change::S) where {T,S}
    f = bridge.decision_function
    bridge.set = change
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(change, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::MOI.ScalarConstantChange) where {T,S}
    f = bridge.decision_function
    # Modify variable part of decision function
    f.variable_part.constant = change.new_constant
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::MOI.ScalarCoefficientChange) where {T,S}
    f = bridge.decision_function
    # Update coefficient in variable part
    modify_coefficient!(f.variable_part.terms, change.variable, change.new_coefficient)
    # Modify the variable part of the constraint as usual
    MOI.modify(model, bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::DecisionCoefficientChange) where {T,S}
    f = bridge.decision_function
    # Update coefficient in decision part
    modify_coefficient!(f.decision_part.terms, change.decision, change.new_coefficient)
    # Update mapped variable through ScalarCoefficientChange
    MOI.modify(model, bridge, MOI.ScalarCoefficientChange(change.decision, change.new_coefficient))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::KnownCoefficientChange) where {T,S}
    f = bridge.decision_function
    i = something(findfirst(t -> t.variable_index == change.known,
                            f.known_part.terms), 0)
    # Update known part of constraint constant
    coefficient = iszero(i) ? zero(T) : f.known_part.terms[i].coefficient
    known_value = MOI.get(model, MOI.VariablePrimal(), change.known)
    f.known_part.constant +=
        (change.new_coefficient - coefficient) * known_value
    # Update coefficient in known part
    modify_coefficient!(f.known_part.terms, change.known, change.new_coefficient)
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::KnownValueChange) where {T,S}
    f = bridge.decision_function
    i = something(findfirst(t -> t.variable_index == change.known,
                            f.known_part.terms), 0)
    if iszero(i)
        # Known value not in objective, nothing to do
        return nothing
    end
    # Update known part of constraint constant
    coefficient = f.known_part.terms[i].coefficient
    f.known_part.constant +=
        coefficient * change.value_difference
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::AffineDecisionConstraintBridge{T,S}, change::KnownValuesChange) where {T,S}
    f = bridge.decision_function
    known_val = zero(T)
    for term in f.known_part.terms
        known_val += term.coefficient * MOI.get(model, MOI.VariablePrimal(), term.variable_index)
    end
    # Update known part of constraint constant
    f.known_part.constant = known_val
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end


text:
using celerite
using Base.Test

include("test_solver.jl")
include("test_cholesky.jl")


text:
blocksizes(A) = BlockSizes(vcat.(size(A))...)
getindex(a::Number, ::Block{0}) = a

blocksizes(A::AbstractTriangular) = blocksizes(parent(A))

blocksizes(A::HermOrSym) = blocksizes(parent(A))
blocksizes(A::AdjOrTrans{<:Any,<:AbstractMatrix}) = BlockSizes(reverse(cumulsizes(blocksizes(parent(A)))))
blocksizes(A::AdjOrTrans{<:Any,<:AbstractVector}) = BlockSizes(([1,2],cumulsizes(blocksizes(parent(A)))[1]))

Base.print_matrix_row(io::IO,
        X::Union{AbstractTriangular{<:Any,<:AbstractBlockMatrix},
                 Symmetric{<:Any,<:AbstractBlockMatrix},
                 Hermitian{<:Any,<:AbstractBlockMatrix}}, A::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString) =
        _blockarray_print_matrix_row(io, X, A, i, cols, sep)


text:
#= test/spfft1.jl
=#

println("spfft1.jl")
@time begin

@testset "spfft1" begin

n = 1500
k =  151

# complex transforms

r = randperm(n)[1:k]

for fcn in (:fft, :bfft)
  sf = Symbol("sp", fcn)
  for T in (ComplexF32, ComplexF64)
    println("  $sf/$T")
    Tr = real(T)

    psf = Symbol("plan_", sf)
    @eval P = $psf($T, $n, $r)

    x  = Array{T }(undef, n)
    y  = Array{T }(undef, k)
    xr = Array{Tr}(undef, size(x))
    yr = Array{Tr}(undef, size(y))

    # f2s
    f2s = Symbol(sf, "_f2s!")
    rand!(x)

    ## (c2c, f2s)
    @eval f = $fcn($x)[$r]
    @eval $f2s($y, P, $x)
    @test f ≈ y

    ## (c2r, f2s)
    @eval $f2s($yr, P, $x)
    @test real(f) ≈ yr

    ## (r2c, f2s)
    xr[:] = real(x)
    @eval f = $fcn($xr)[$r]
    @eval $f2s($y, P, $xr)
    @test f ≈ y

    # s2f
    s2f = Symbol(sf, "_s2f!")
    rand!(y)

    ## (c2c, s2f)
    fill!(x, 0)
    x[r] = y
    @eval f = $fcn($x)
    @eval $s2f($x, P, $y)
    @test f ≈ x

    ## (c2r, s2f)
    @eval $s2f($xr, P, $y)
    @test real(f) ≈ xr

    ## (r2c, s2f)
    yr[:] = real(y)
    fill!(x, 0)
    x[r] = yr
    @eval f = $fcn($x)
    @eval $s2f($x, P, $yr)
    @test f ≈ x
  end
end

# real transforms

m = div(n,2) + 1
r = vcat([1,m], 1 .+ randperm(m-2)[1:k-2])  # include edge cases

for T in (Float32, Float64)
  Tc = Complex{T}

  x  = Array{T }(undef, n)
  y  = Array{Tc}(undef, k)
  xc = Array{Tc}(undef, m)

  # (r2c, f2s)
  println("  sprfft/$T")
  P = plan_sprfft(T, n, r)
  rand!(x)
  f = rfft(x)[r]
  sprfft_f2s!(y, P, x)
  @test f ≈ y

  # (c2r, s2f)
  println("  spbrfft/$T")
  P = plan_spbrfft(T, n, r)
  rand!(y)
  y[r.==1] = real(y[r.==1])
  even = n % 2 == 0
  even && (y[r.==m] = real(y[r.==m]))
  fill!(xc, 0)
  xc[r] = y
  f = brfft(xc, n)
  spbrfft_s2f!(x, P, y)
  @test f ≈ x
end

end
end


text:
function auto_queue(mission::Mission)
    master_df = master_query_public(mission)

    queue = filter(x->x[:publicity]    == true,  master_df)
    queue = filter(x->x[:report_exists]== false, queue)
    queue = filter(x->x[:obs_type]     != "CAL", queue)
    queue = filter(x->x[:error]        == false, queue)

    return queue
end

"""
    auto_report(::Mission; limit::Union{Bool,Int}, update::Bool, nuke::Bool)

Calls `auto_queue` function to generate a queue of reports to make, the queue filters:

    * Public-only

    * Reportless

    * Not 'CAL' type observations

    * Error free

Leaving only suitable observations to be analysed

Calls 'report_all' using the queued observations

Will continue to generate reports until the `limit` is reached (if there is one)
"""
function auto_report(mission::Mission; limit=false, update=true, nuke=false)
    if update
        master(mission; update=true)
    end

    obs_queue = auto_queue(mission)

    i = 0
    if typeof(obs_queue) == DataFrames.DataFrameRow
        report_all(mission, obs_queue, update_masterpage=false)
    else
        for obs_row in DataFrames.eachrow(obs_queue)
            println(repeat("-", 8), "  ", obs_row[:obsid], "  ", repeat("-", 8), "\n")

            errors_dl = JAXTAM._log_query(mission, obs_row, "errors", :download; surpress_warn=true)

            if !ismissing(errors_dl)
                @warn "Download errors found, skipping"
                continue
            end

            report_all(mission, obs_row; nuke=nuke, update_masterpage=false)

            i+=1
            if i == limit
                break
            end
        end
    end

    master_append(mission; update=true)
    webgen_mission(mission)
end


text:
module OneBody

# include("../../parameter.jl")
include("../../utility/constant.jl")
using QuantumStatistics: Green, FastMath
using StaticArrays
const Mom = Main.Mom
include("./update.jl")

struct Group
    diag::Int
    order::Int
    neighbor::Array{Group,1} # define which groups are connected to the correct one
    reweight::Float
    innerK::Tuple{Int,Int}
    innerT::Tuple{Int,Int}
    statistics::Array{Float,2} # two dimension array of (k, tau)
    ########### extra ################
    tauGrid::Any
    KGrid::Any
    function Group(diag, order, grids)
        if diag == Diag.∅
            @assert order == 0 "∅ diagram must be order 0!"
            innerK, innerT = (2, 1), (2, 1)
            tauGrid = nothing
            KGrid = nothing
        else
            @assert order > 0 "$diag diagram must be order ≥ 0!"
            tauGrid = grids.tau
            innerK, innerT = (2, order + 1), (2, order)
            if (diag == Diag.Σ) || (diag == Diag.Δ)
                kGrid = grids.fermiK
            elseif diag == Diag.Π
                kGrid = grids.boseK
            else
                error("$diag diagram is not an implemented one-body observable!")
            end
        end

        data = zeros(Flat, tauGrid.size, kGrid.size)
        return new(diag, 0, [], 1.0, innerK, innerT, data, tauGrid, kGrid)
    end
end

"""
two groups are equal if and only if the diagram and the order are the same
"""
Base.:=(x::Group, y::Group) = (x.diag == y.diag && x.order == y.order)

mutable struct State
    group::Group
    absweight::Float
    T::Array{Float,1}
    K::Array{Mom,1}
    propose::Dict{Symbol,Float}
    accept::Dict{Symbol,Float}
    visited::Dict{Group,Float}
    #### extra #############
    extTidx::Int
    extKidx::Int

    """
    MC Variable initialize here
    """
    function State(para, groups)
        LastT = para.order + 2
        LastK = para.order + 1
        varT = rand(LastT) * para.β
        varK = [rand(Mom) * para.Kf for i = 1:LastK]

        curr = groups[1]
        varT[1] = 0.0
        varT[LastT] = curr.tauGrid[curr.extTidx]
        varK[1] = zero(Mom)
        varK[1][1] = curr.KGrid[curr.extKidx]

        state = new(0, curr, 0.0, LastT, varT, varK, {}, {}, {}, 1, 1)
        state.absweight = abs(evaluate(para, state))
        return state
    end
end

"""
    init(parameter, grids)

build Diagram groups and the State from the para.diagrams and other parameters
"""
function init(para, grids)
    # maxOrder=maximum([d[2] for d in diagList])
    lastT = para.order + 1
    groups = [Group(Diag.∅, 0, grids), ]

    for o in 1:para.order
        for diag in para.diagrams
            push!(groups, Group(diag, o, grids))
        end
    end
    state = State(para, groups)
    return state, groups, (propGroup, rejectGroup)
end

function evaluate(para, state)
    group = state.group
    if group.diag == Diag.∅
        return 1.0
    elseif group.diag == Diag.Π
        if group.order == 1
            tau = state.T[state.lastT] - state.T[1]
            q, k = state.K[1], state.K[2]
            gweight = bareFermi(para.β, tau, k, para.μ) * bareFermi(para.β, tau, k + q, para.μ)
            return -para.spin * gweight / (2π)^para.dim
        else
            error("Not implemented")
        end
    else
        error("Not implemented")
    end
end

function measure(para, state)
    group = state.group
    if group.diag == Diag.∅
    else
    end
end

end

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
using deSolveDiffEq
using Test

function lorenz(u,p,t)
 du1 = 10.0(u[2]-u[1])
 du2 = u[1]*(28.0-u[3]) - u[2]
 du3 = u[1]*u[2] - (8/3)*u[3]
 [du1, du2, du3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,deSolveDiffEq.lsoda())
sol = solve(prob,deSolveDiffEq.lsode())
sol = solve(prob,deSolveDiffEq.lsodes())
sol = solve(prob,deSolveDiffEq.lsodar())
sol = solve(prob,deSolveDiffEq.vode())
sol = solve(prob,deSolveDiffEq.daspk())
sol = solve(prob,deSolveDiffEq.euler())
sol = solve(prob,deSolveDiffEq.rk4())
sol = solve(prob,deSolveDiffEq.ode23())
sol = solve(prob,deSolveDiffEq.ode45())
sol = solve(prob,deSolveDiffEq.radau())
sol = solve(prob,deSolveDiffEq.bdf())
sol = solve(prob,deSolveDiffEq.bdf_d())
sol = solve(prob,deSolveDiffEq.adams())
sol = solve(prob,deSolveDiffEq.impAdams())
sol = solve(prob,deSolveDiffEq.impAdams_d())
#sol = solve(prob,deSolveDiffEq.iteration())

sol = solve(prob,deSolveDiffEq.lsoda(),saveat=0.1)
#using Plots; plot(sol,vars=(1,2,3))


text:
using Convex, GLPKMathProgInterface
include("antidiag.jl")

n = 8
x = Variable((n, n), :Bin)
# We can have at most one queen on any anti-diagonal
constr = Constraint[sum(antidiag(x, k)) <= 1 for k = -n+2:n-2]
# We can have at most one queen on any diagonal
constr += Constraint[sum(diag(x, k)) <= 1 for k = -n+2:n-2]
# We must have one queen per row and one queen per column
constr += Constraint[sum(x, 1) == 1, sum(x, 2) == 1]
p = satisfy(constr)
solve!(p, GLPKSolverMIP())

# Testing
for k = -n+2:n-2
	@assert evaluate(sum(antidiag(x, k))) <= 1
	@assert evaluate(sum(diag(x, k))) <= 1
end
@assert all(evaluate(sum(x, 1)) .== 1)
@assert all(evaluate(sum(x, 2)) .== 1)


text:
using ADCME
using AdFem

include("chip_unstructured_solver.jl")
include("chip_unstructured_geometry.jl")
include("plot_inverse_one_iter.jl")

trialnum = 3

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

# θ0 = [2.60475; 0.49; 1.0]
# θ = Variable([2.60475; 0.49; 1.0])
# paramA = θ[1]
# paramB = θ[2]
# paramC = θ[3]

function k_exact(x, y)
    k_mold + 1000 * k_chip_ref * (x-0.49)^2 / (1 + x^2)
end

function k_nn(x, y, θ)
    constant(@. k_exact(x, y)) +  fc(x, [20,20,20,1], θ)
end

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 14    # number of iterations for Newton's method, 8 is good for m=400

n = 901  #ae_num([20,20,20,1]) => 861?
θ = Variable(randn(n))
θ = θ ./ norm(θ)

xy = mesh.nodes
xy2 = zeros(mesh.nedge, 2)
for i = 1:mesh.nedge
    xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
end
xy = [xy;xy2]

x, y = xy[chip_fem_idx, 1], xy[chip_fem_idx, 2]
k_chip =  k_nn(x, y, θ)
k_chip_exact = constant( @. k_exact(x, y) )

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

B = constant(compute_interaction_matrix(mesh))
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
      bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

S0 = constant(zeros(nelem+3*ndof))
S = solve_navier_stokes(S0, NT, k_chip)
S_computed = S[end, :]
S_data = matread("fn1/xz_chip_unstructured_data.mat")["V"]

loss =  mean((S_computed .- S_data)^2)
loss = loss * 1e10
sess = Session(); init(sess)

# ---------------------------------------------------
n = 901
θ0 = zeros(n)
@info run(sess, loss, θ=>θ0)
lineview(sess, θ, loss, θ0, ones(n))
savefig("fn$trialnum/lineview.png")
gradview(sess, θ, loss, zeros(n))
savefig("fn$trialnum/gradview.png")


text:
# Example for using the 1x1 convolution operator to permute an image along the channel dimension.
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test
using Flux

# Dimensions
nx = 64 # no. of pixels in x dimension
ny = 64 # no. of pixels in y dimension
k = 10  # no. of channels
batchsize = 4

# Input image: nx x ny x k x batchsize
X = glorot_uniform(nx, ny, k, batchsize)

# 1x1 convolution operators
C = Conv1x1(k)
C0 = Conv1x1(k)

# Generate "true/observed" data with the same dimension as X
Y = C.forward(X)

# Predicted data
Y0 = C0.forward(X)
@test isnothing(C0.v1.grad) # after forward pass, gradients are not set

# Data residual
ΔY = Y0 - Y

# Backward pass: Pass ΔY to compute ΔX, the gradient w.r.t the input X.
# Also pass Y0 to recompute the forward state X using the inverse mapping
# and use it to compute the derivative w.r.t. the coefficients of the
# Householder matrix.
ΔX, X_ = C0.inverse((ΔY, Y0))
@test ~isnothing(C0.v1.grad)    # after inverse pass, gradients are set
@test isapprox(norm(X - X_)/norm(X), 0f0, atol=1f-6)    # X and X_ should be the same


text:
# Filesystem

@! function Base.open(filename::AbstractString, mode::AbstractString; kws...)
    io = open(filename, mode; kws...)
    @defer close(io)
    io
end

@! function Base.open(filename::AbstractString; kws...)
    io = open(filename; kws...)
    @defer close(io)
    io
end

@! function Base.mktemp(parent=tempdir())
    (path, io) = mktemp(parent, cleanup=false)
    @defer begin
        close(io)
        rm(path, force=true)
    end
    (path, io)
end

@! function Base.mktempdir(parent=tempdir(); prefix="jl_")
    path = mktempdir(parent, cleanup=false)
    @defer rm(path, recursive=true, force=true)
    path
end

@! function Base.cd(path::AbstractString=homedir())
    oldpath = pwd()
    cd(path)
    @defer cd(oldpath)
    # TODO: `Base.cd(::Function, path)` has some magic here for unix which we
    # could replicate.
end

# Locks

@! function Base.lock(lk::Base.AbstractLock)
    lock(lk)
    @defer unlock(lk)
end

# Standard streams

# Incompatibility due to
# https://github.com/JuliaLang/julia/pull/39132
@static if VERSION < v"1.7"

@! function Base.redirect_stdout(stream)
    prev_stream = stdout
    x = redirect_stdout(stream)
    @defer redirect_stdout(prev_stream)
    x
end

@! function Base.redirect_stderr(stream)
    prev_stream = stderr
    x = redirect_stderr(stream)
    @defer redirect_stderr(prev_stream)
    x
end

@! function Base.redirect_stdin(stream)
    prev_stream = stdin
    x = redirect_stdin(stream)
    @defer redirect_stdin(prev_stream)
    x
end

else

@! function (f::Base.RedirectStdStream)(stream)
    # See https://github.com/JuliaLang/julia/blob/294b0dfcd308b3c3f829b2040ca1e3275595e058/base/stream.jl#L1417
    stdold = f.unix_fd == 0 ? stdin :
             f.unix_fd == 1 ? stdout :
             f.unix_fd == 2 ? stderr :
             throw(ArgumentError("Not implemented to get old handle of fd except for stdio"))
    x = f(stream)
    @defer f(stdold)
    x
end

@! function Base.redirect_stdio(; kws...)
    @! enter_do(redirect_stdio; kws...)
end

end


text:
""" Test JuMP model with a variety of variable bounds and sizes. """
function test_model()
    model = Model()
    @variables(model, begin
        -5 <= x[1:5] <= 5
        -4 <= y[1:3] <= 1
        -30 <= z
        -2 <= a[1:2, 1:3]
    end)
    JuMP.set_upper_bound(a[2,2], 3)
    JuMP.set_upper_bound(a[1,1], 3)
    return model, x, y, z, a
end

""" Basic qp for testing. """
function test_qp(solver = OCTHaGOn.SOLVER_SILENT)
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, 2 >= x[1:2] >= 0)
    @objective(m, Min, 3*x[1]^2 + x[2]^2 + 2*x[1]*x[2] + x[1] + 6*x[2] + 2)
    @constraint(m, 2*x[1] + 3*x[2] >= 4)
    return m
end

""" Basic gqp, same as above. """
function test_gqp(solver = OCTHaGOn.SOLVER_SILENT)
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, 2 >= x[1:2] >= 0)
    @constraint(m, 2*x[1] + 3*x[2] >= 4)
    @variable(m, obj)
    @objective(m, Min, obj)
    gm = GlobalModel(model = m)
    add_nonlinear_constraint(gm, :(x -> 3*x[1]^2 + x[2]^2 + 2*x[1]*x[2] + x[1] + 6*x[2] + 2), dependent_var = obj)
    return gm
end

"""
    random_qp(dims::Int64, nconstrs::Int64, sparsity=dims, solver = OCTHaGOn.SOLVER_SILENT)

Generates a random quadratic program with the specified parameters.
"""
function random_qp(dims::Int64, nconstrs::Int64, sparsity=dims, solver = OCTHaGOn.SOLVER_SILENT)
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, -1 <= x[1:dims] <= 1)
    m[:sparse_vars] = randperm(dims)[1:sparsity]
    m[:A] = rand(dims, sparsity).*2 .- 1
    m[:b] = rand(dims)
    m[:C] = rand(nconstrs, dims).*2 .- 1
    m[:d] = rand(nconstrs)
    @constraint(m, m[:C]*x .>= m[:d])
    @objective(m, Min, sum((m[:A]*x[m[:sparse_vars]] - m[:b]).^2))
    return m
end

""" Turns random_qp into a GlobalModel. """
function gmify_random_qp(m::JuMP.Model)
    @variable(m, obj)
    @objective(m, Min, obj)
    gm = GlobalModel(model = m)
    add_nonlinear_constraint(gm, :(x -> sum(($(m[:A])*x[$(m[:sparse_vars])] - $(m[:b])).^2)), vars = m[:x][m[:sparse_vars]],
                                    dependent_var = obj, name = "qp")
    return gm
end


text:
export miccalibrate, miccalibrate!, load_angular_measurements, load_vituix_spinorama_curves, load_spinorama_from_FRDs, load_vituix_FIRs

function miccalibrate!(a::AbstractSpectrumArray, calfile; skipto=0)
    # df = CSV.read(calfile, DataFrame; header=false, skipto=skipto)
    df = CSV.read(calfile, DataFrame; header=false)
    freqs, magdb = df.Column1, df.Column2
    mag = LinearInterpolation(freqs, db2amp.(magdb); extrapolation_bc = Flat())(domain(a))
    a ./ Complex.(mag)
end

miccalibrate(a::AbstractSpectrumArray, calfile; skipto=0) = miccalibrate!(deepcopy(a), calfile; skipto=skipto)

"""
    load_directivity_spectra(rootdir, dirtemplate, loader, select=:all)

Load all directivity spectra in `irdir` with each spectrum in a separate IR file/directory matching regex `template`.

Use `loader(f)` function to get the a subtype of the `AbstractSpectrumArray`. Imported spectra can be limited to only `:positive` or `:negative` angles using `select`.
Degrees are infered as to be the first group of `ir_template`. Example: `r"ir_(\\d\\d\\d).wav"`.
"""
function load_angular_measurements(irdir, template, loader, select=:all)
    @assert select ∈ [:all, :positive, :negative]
    spectra = []
    angles = Float64[]
    for f in readdir(irdir)
        fpath = joinpath(irdir, f)
        m = match(template, f)
        if !isnothing(m)
            angle = normdeg(parse(Float64, m.captures[1]))
            if select == :positive && angle < 0
                continue
            end
            if select == :negative && angle > 0
                continue
            end
            spectrum = loader(fpath)
            names!(spectrum, [Symbol(round(angle) ≈ angle ? round(Int, angle) : angle)])
            push!(spectra, spectrum)
            push!(angles, angle)
        end
    end
    length(spectra) > 0 || error("no file imported!")
    catspectra = hcat(spectra...)
    AngularMeasurements(catspectra, angles)
end

function load_vituix_spinorama_curves(rootdir; samplerate=44100*Unitful.Hz)
    function load(file_name, channel_name)
        readfrd(joinpath(rootdir, file_name), samplerate, magf=db2amp, name=channel_name) .|> abs
    end
    pairs = [
        "VituixCAD_CTA-2034-A Directivity Index.frd" => :SPDI,
        "VituixCAD_CTA-2034-A ER Ceiling.frd" => :early_reflections_ceiling,
        "VituixCAD_CTA-2034-A ER Floor.frd" => :early_reflections_floor,
        "VituixCAD_CTA-2034-A ER Front.frd" => :early_reflections_front,
        "VituixCAD_CTA-2034-A ER Horizontal.frd" => :early_reflections_horizontal,
        "VituixCAD_CTA-2034-A ER Rear.frd" => :early_reflections_rear,
        "VituixCAD_CTA-2034-A ER Side.frd" => :early_reflections_side,
        "VituixCAD_CTA-2034-A ER Total.frd" => :early_reflections,
        "VituixCAD_CTA-2034-A ER Vertical.frd" => :early_reflections_vertical,
        "VituixCAD_CTA-2034-A ERDI Horizontal.frd" => :ERDI_horizontal,
        "VituixCAD_CTA-2034-A ERDI Total.frd" => :ERDI,
        "VituixCAD_CTA-2034-A ERDI Vertical.frd" => :ERDI_vertical,
        "VituixCAD_CTA-2034-A In-room response.frd" => :in_room,
        "VituixCAD_CTA-2034-A Listening window.frd" => :listening_window,
        "VituixCAD_CTA-2034-A Power response.frd" => :power_response,
        "VituixCAD_CTA-2034-A Reference angle.frd" => :reference,
    ]
    hcat(map(p -> load(p...), pairs)...)
end

function load_spinorama_from_FRDs(rootdir;
    template_hor=Regex(raw"VituixCAD_PolarFR hor (.*)\.frd"), template_ver=Regex(raw"VituixCAD_PolarFR ver (.*)\.frd"), samplerate=44100*Unitful.Hz)
    loader = f->readfrd(f, samplerate, magf=db2amp)
    hor = load_angular_measurements(rootdir, template_hor, loader, :all)
    ver = load_angular_measurements(rootdir, template_ver, loader, :all)
    return Spinorama(hor, ver)
end

function load_vituix_FIRs(wavs...; Δt=0.0, names=nothing)
    if !isnothing(names) && length(wavs) != length(names)
        throw(ArgumentError("the number of WAV files does not match the number of names!"))
    end
    names_ = isnothing(names) ? [Symbol("FIR_$i") for i in 1:length(wavs)] : names
    Hs = hcat((readwav_rfft(wav)[2] for wav in wavs)...)
    names!(Hs, names_)
    Hmax = maximum(abs.(Hs))
    Δt ≈ 0.0 || delay!(Hs, Δt)
    return Hs ./ Hmax
end