text:
using FastGaussQuadrature, Printf

rtol(x, y) = maximum(abs.(x - y)) / maximum(abs.(x))

# SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, grid, type, dlr.symmetry, rtol = dlr.rtol, degree = 24, regularized = true)
SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree = 24, regularized = true)

function MultiPole(dlr, grid, type)
    Euv = dlr.Euv
    poles = [-Euv, -0.2 * Euv, 0.0, 0.8 * Euv, Euv]
    # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true)
    return Sample.MultiPole(dlr, type, poles, grid; regularized = true)
end

function compare(case, a, b, eps, requiredratio, para = "")
    err = rtol(a, b)
    ratio = isfinite(err) ? Int(round(err / eps)) : 0
    if ratio > 50
        printstyled("$case, $para err: ", color = :white)
        printstyled("$(round(err, sigdigits=3)) = $ratio x rtol\n", color = :green)
    end
    @test rtol(a, b) .< requiredratio * eps # dlr should represent the Green's function up to accuracy of the order eps
end

function compare_atol(case, a, b, atol, para = "")
    err = rtol(a, b)
    err = isfinite(err) ? round(err, sigdigits = 3) : 0
    if err > 100 * atol
        printstyled("$case, $para err: ", color = :white)
        printstyled("$err = $(Int(round(err/atol))) x atol\n", color = :blue)
    end
    @test rtol(a, b) .< 5000 * atol # dlr should represent the Green's function up to accuracy of the order eps
end

@testset "Correlator Representation" begin

    function test(case, isFermi, symmetry, Euv, β, eps)
        # println("Test $case with isFermi=$isFermi, Symmetry = $symmetry, Euv=$Euv, β=$β, rtol=$eps")

        para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps"

        dlr = DLRGrid(Euv, β, eps, isFermi, symmetry) #construct dlr basis
        dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose

        #=========================================================================================#
        #                              Imaginary-time Test                                        #
        #=========================================================================================#
        # get imaginary-time Green's function
        Gdlr = case(dlr, dlr.τ, :τ)
        # get imaginary-time Green's function for τ sample
        τSample = dlr10.τ
        Gsample = case(dlr, τSample, :τ)

        ########################## imaginary-time to dlr #######################################
        coeff = tau2dlr(dlr, Gdlr)
        Gfitted = dlr2tau(dlr, coeff, τSample)
        compare("dlr τ → dlr → generic τ $case", Gsample, Gfitted, eps, 100, para)
        # for (ti, t) in enumerate(τSample)
        #     @printf("%32.19g    %32.19g   %32.19g   %32.19g\n", t / β, Gsample[1, ti],  Gfitted[1, ti], Gsample[1, ti] - Gfitted[1, ti])
        # end

        compare("generic τ → dlr → τ $case", tau2tau(dlr, Gsample, dlr.τ, τSample), Gdlr, eps, 1000, para)
        #=========================================================================================#
        #                            Matsubara-frequency Test                                     #
        #=========================================================================================#
        # #get Matsubara-frequency Green's function
        Gndlr = case(dlr, dlr.n, :n)
        nSample = dlr10.n
        Gnsample = case(dlr, nSample, :n)

        # #Matsubara frequency to dlr
        coeffn = matfreq2dlr(dlr, Gndlr)
        Gnfitted = dlr2matfreq(dlr, coeffn, nSample)
        #     for (ni, n) in enumerate(nSample)
        #     @printf("%32.19g    %32.19g   %32.19g   %32.19g\n", n, real(Gnsample[1, ni]),  real(Gnfitted[1, ni]), abs(Gnsample[1, ni] - Gnfitted[1, ni]))
        # end

        compare("dlr iω → dlr → generic iω $case ", Gnsample, Gnfitted, eps, 100, para)
        compare("generic iω → dlr → iω $case", matfreq2matfreq(dlr, Gnsample, dlr.n, nSample), Gndlr, eps, 1000, para)

        #=========================================================================================#
        #                            Fourier Transform Test                                     #
        #=========================================================================================#
        Gnfourier = tau2matfreq(dlr, Gdlr, nSample)
        compare("τ→dlr→iω $case", Gnsample, Gnfourier, eps, 1000, para)
        # for (ti, t) in enumerate(nSample)
        #     @printf("%32.19g    %32.19g   %32.19g   %32.19g\n", t / β, imag(Gnsample[2, ti]), imag(Gnfourier[2, ti]), abs(Gnsample[2, ti] - Gnfourier[2, ti]))
        # end

        Gfourier = matfreq2tau(dlr, Gndlr, τSample)
        compare("iω→dlr→τ $case", Gsample, Gfourier, eps, 1000, para)
        # for (ti, t) in enumerate(τSample)
        #     @printf("%32.19g    %32.19g   %32.19g   %32.19g\n", t / β, Gsample[2, ti],  real(Gfourier[2, ti]), abs(Gsample[2, ti] - Gfourier[2, ti]))
        # end

        #=========================================================================================#
        #                            Noisey data Test                                             #
        #=========================================================================================#

        # err = 10 * eps
        atol = eps
        noise = atol * rand(eltype(Gsample), length(Gsample))
        GNoisy = Gsample .+ noise
        compare_atol("noisy generic τ → dlr → τ $case", tau2tau(dlr, GNoisy, dlr.τ, τSample; error = abs.(noise)), Gdlr, atol, para)

        noise = atol * rand(eltype(Gnsample), length(Gnsample))
        GnNoisy = Gnsample .+ noise
        compare_atol("noisy generic iω → dlr → iω $case", matfreq2matfreq(dlr, GnNoisy, dlr.n, nSample, error = abs.(noise)), Gndlr, atol, para)
    end

    # the accuracy greatly drops beyond Λ >= 1e8 and rtol<=1e-6
    cases = [SemiCircle, MultiPole]
    Λ = [1e3, 1e5, 1e7]
    rtol = [1e-8, 1e-10]
    for case in cases
        for l in Λ
            for r in rtol
                test(case, true, :none, 1.0, l, r)
                test(case, false, :none, 1.0, l, r)
                test(case, false, :ph, 1.0, l, r)
                test(case, true, :ph, 1.0, l, r)
                test(case, false, :pha, 1.0, l, r)
                test(case, true, :pha, 1.0, l, r)
            end
        end
    end

end

@testset "Tensor ↔ Matrix Mapping" begin
    a = rand(3)
    acopy = deepcopy(a)
    b, psize = Lehmann._tensor2matrix(a, 1)
    anew = Lehmann._matrix2tensor(b, psize, 1)
    @test acopy ≈ anew

    a = rand(3, 4)
    acopy = deepcopy(a)
    for axis = 1:2
        b, psize = Lehmann._tensor2matrix(a, axis)
        anew = Lehmann._matrix2tensor(b, psize, axis)
        @test acopy ≈ anew
    end

    a = rand(3, 4, 5)
    acopy = deepcopy(a)
    for axis = 1:3
        b, psize = Lehmann._tensor2matrix(a, axis)
        anew = Lehmann._matrix2tensor(b, psize, axis)
        @test acopy ≈ anew
    end
end

@testset "Tensor DLR" begin
    Euv, β = 1.0, 1000.0
    eps = 1e-10
    isFermi = true
    symmetry = :none
    # symmetry = :ph
    weight = π / 2 * Euv
    para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps"
    dlr = DLRGrid(Euv, β, eps, isFermi, symmetry) #construct dlr basis

    Gτ = SemiCircle(dlr, dlr.τ, :τ)
    Gn = SemiCircle(dlr, dlr.n, :n)
    Gτ_copy = deepcopy(Gτ)
    Gn_copy = deepcopy(Gn)

    n1, n2 = 1, 1
    a = rand(n1, n2)
    # a = ones(n1, n2)
    sumrule_τ = zeros(n1, n2)
    sumrule_n = zeros(n1, n2)
    tensorGτ = zeros(eltype(Gτ), (n1, n2, length(Gτ)))
    for (gi, g) in enumerate(Gτ)
        tensorGτ[:, :, gi] = a .* g
        sumrule_τ = a * weight
    end
    tensorGn = zeros(eltype(Gn), (n1, n2, length(Gn)))
    for (gi, g) in enumerate(Gn)
        tensorGn[:, :, gi] = a .* g
        sumrule_n = a * weight
    end
    tensorGτ_copy = tensorGτ
    tensorGn_copy = tensorGn

    compare("τ ↔ iω tensor", tau2matfreq(dlr, Gτ_copy), Gn, eps, 1000.0, para)
    @test Gτ ≈ Gτ_copy #make sure there is no side effect on G
    compare("iω ↔ τ tensor", matfreq2tau(dlr, Gn_copy), Gτ, eps, 1000.0, para)
    @test Gn ≈ Gn_copy #make sure there is no side effect on G

    compare("τ ↔ iω tensor", tau2matfreq(dlr, Gτ_copy; sumrule = weight), Gn, eps, 1000.0, para)
    @test Gτ ≈ Gτ_copy #make sure there is no side effect on G
    compare("iω ↔ τ tensor", matfreq2tau(dlr, Gn_copy; sumrule = weight), Gτ, eps, 1000.0, para)
    @test Gn ≈ Gn_copy #make sure there is no side effect on G

    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis = 3), tensorGn, eps, 1000.0, para)
    @test tensorGτ ≈ tensorGτ_copy #make sure there is no side effect on G
    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis = 3), tensorGτ, eps, 1000.0, para)
    @test tensorGn ≈ tensorGn_copy #make sure there is no side effect on G

    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis = 3, sumrule = sumrule_τ), tensorGn, eps, 1000.0, para)
    @test tensorGτ ≈ tensorGτ_copy #make sure there is no side effect on G
    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis = 3, sumrule = sumrule_n), tensorGτ, eps, 1000.0, para)
    @test tensorGn ≈ tensorGn_copy #make sure there is no side effect on G

end

@testset "Least square fitting" begin
    Gτ = [6.0, 5.0, 7.0, 10.0]
    kernel = zeros(4, 2)
    kernel[:, 1] = [1.0, 1.0, 1.0, 1.0]
    kernel[:, 2] = [1.0, 2.0, 3.0, 4.0]
    dlrGrid = DLRGrid(β = 10.0, isFermi = true)
    coeff = Lehmann._weightedLeastSqureFit(dlrGrid, Gτ, nothing, kernel, nothing)
    @test coeff ≈ [3.5, 1.4]
end

@testset "DLR io" begin
    function finddlr(folder, filename)
        searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
        file = searchdir(folder, filename)
        if length(file) == 1
            #dlr file found
            println("find dlr file: ", file[1])
            return joinpath(folder, file[1])
        end
        return nothing
    end

    folder = "./"
    Euv, β, rtol = 1.5, 110.0, 4.3e-8

    #save dlr to a local file
    dlr = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild = true, folder = folder, verbose = false)
    file = finddlr(folder, ".dlr")
    @test isnothing(file) == false
    #load dlr from the local file
    dlr_load = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild = false, folder = folder, verbose = false)
    @test dlr_load.τ ≈ dlr.τ
    @test dlr_load.n ≈ dlr.n
    @test dlr_load.ω ≈ dlr.ω
    @test dlr_load.ωn ≈ dlr.ωn
    @test dlr_load.Euv ≈ dlr.Euv
    @test dlr_load.β ≈ dlr.β
    @test dlr_load.rtol ≈ dlr.rtol
    @test dlr_load.Λ ≈ dlr.Λ
    rm(file)
end

@testset "JLD2 IO" begin
    dlr = DLRGrid(isFermi = true, beta = 10.0)
    save("test.jld2", Dict("dlr" => dlr), compress = true)
    dlr_new = load("test.jld2")["dlr"]
    println(dlr_new)
    @test dlr.isFermi == dlr_new.isFermi
    @test dlr.Euv ≈ dlr_new.Euv
    @test dlr.β ≈ dlr_new.β
    @test dlr.τ ≈ dlr_new.τ
    @test dlr.n ≈ dlr_new.n
    @test dlr.ωn ≈ dlr_new.ωn
    @test dlr.symmetry == dlr_new.symmetry
    rm("test.jld2")
end


text:
module Day6

export day6

function answersbygroup(filename)
    groups = Vector{Vector{String}}()
    group = Vector{String}()
    for line in readlines(filename)
        line = strip(line)
        if line == ""
            push!(groups, group)
            group = String[]
        else
            push!(group, line)
        end
    end
    push!(groups, group)
    groups
end

function day6(::Val{:a}, filename)
    groups = answersbygroup(filename)
    sum(length(unique(join(group))) for group in groups)
end

function check(c, group)
    all(in(c,str) for str in group)
end

function day6(::Val{:b}, filename)
    groups = answersbygroup(filename)
    sum(filter(c -> check(c,group), group[1]) |> unique |> length for group in groups)
end

end

ext: jl


text:
using BitsMNIST
using Test
using Random

const SEED = 10

include("./datasets.jl")
include("./statistics.jl")
include("./zerone.jl")
include("./io.jl")


text:
#==============================================================================#
# ServerlessApplicationRepository.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/serverlessrepo-2017-09-08.normal.json
#==============================================================================#

__precompile__()

module ServerlessApplicationRepository

using AWSCore

"""
    using AWSSDK.ServerlessApplicationRepository.create_application
    create_application([::AWSConfig], arguments::Dict)
    create_application([::AWSConfig]; author=, description=, name=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "POST", "/applications", arguments::Dict)
    serverlessrepo([::AWSConfig], "POST", "/applications", author=, description=, name=, <keyword arguments>)

# CreateApplication Operation

Creates an application, optionally including an AWS SAM file to create the first application version in the same call.

# Arguments

## `author = ::String` -- *Required*
The name of the author publishing the app.

Minimum length=1. Maximum length=127.

Pattern "^[a-z0-9](([a-z0-9]|-(?!-))*[a-z0-9])?\$";

## `description = ::String` -- *Required*
The description of the application.

Minimum length=1. Maximum length=256

## `homePageUrl = ::String`
A URL with more information about the application, for example the location of your GitHub repository for the application.

## `labels = [::String, ...]`
Labels to improve discovery of apps in search results.

Minimum length=1. Maximum length=127. Maximum number of labels: 10

Pattern: "^[a-zA-Z0-9+\\\\-_:\\\\/@]+\$";

## `licenseBody = ::String`
A local text file that contains the license of the app that matches the spdxLicenseID value of your application. The file is of the format file://<path>/<filename>.

Maximum size 5 MB

Note: Only one of licenseBody and licenseUrl can be specified, otherwise an error will result.

## `licenseUrl = ::String`
A link to the S3 object that contains the license of the app that matches the spdxLicenseID value of your application.

Maximum size 5 MB

Note: Only one of licenseBody and licenseUrl can be specified, otherwise an error will result.

## `name = ::String` -- *Required*
The name of the application that you want to publish.

Minimum length=1. Maximum length=140

Pattern: "[a-zA-Z0-9\\\\-]+";

## `readmeBody = ::String`
A local text readme file in Markdown language that contains a more detailed description of the application and how it works. The file is of the format file://<path>/<filename>.

Maximum size 5 MB

Note: Only one of readmeBody and readmeUrl can be specified, otherwise an error will result.

## `readmeUrl = ::String`
A link to the S3 object in Markdown language that contains a more detailed description of the application and how it works.

Maximum size 5 MB

Note: Only one of readmeBody and readmeUrl can be specified, otherwise an error will result.

## `semanticVersion = ::String`
The semantic version of the application:

[https://semver.org/](https://semver.org/)

## `sourceCodeUrl = ::String`
A link to a public repository for the source code of your application.

## `spdxLicenseId = ::String`
A valid identifier from [https://spdx.org/licenses/](https://spdx.org/licenses/).

## `templateBody = ::String`
The local raw packaged AWS SAM template file of your application. The file is of the format file://<path>/<filename>.

Note: Only one of templateBody and templateUrl can be specified, otherwise an error will result.

## `templateUrl = ::String`
A link to the S3 object cotaining the packaged AWS SAM template of your application.

Note: Only one of templateBody and templateUrl can be specified, otherwise an error will result.

# Returns

`CreateApplicationResponse`

# Exceptions

`TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException`, `ConflictException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/CreateApplication)
"""
@inline create_application(aws::AWSConfig=default_aws_config(); args...) = create_application(aws, args)

@inline create_application(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "POST", "/applications", args)

@inline create_application(args) = create_application(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.create_application_version
    create_application_version([::AWSConfig], arguments::Dict)
    create_application_version([::AWSConfig]; applicationId=, semanticVersion=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "PUT", "/applications/{applicationId}/versions/{semanticVersion}", arguments::Dict)
    serverlessrepo([::AWSConfig], "PUT", "/applications/{applicationId}/versions/{semanticVersion}", applicationId=, semanticVersion=, <keyword arguments>)

# CreateApplicationVersion Operation

Creates an application version.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `semanticVersion = ::String` -- *Required*
The semantic version of the new version.

## `sourceCodeUrl = ::String`
A link to a public repository for the source code of your application.

## `templateBody = ::String`
The raw packaged AWS SAM template of your application.

## `templateUrl = ::String`
A link to the packaged AWS SAM template of your application.

# Returns

`CreateApplicationVersionResponse`

# Exceptions

`TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException`, `ConflictException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/CreateApplicationVersion)
"""
@inline create_application_version(aws::AWSConfig=default_aws_config(); args...) = create_application_version(aws, args)

@inline create_application_version(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "PUT", "/applications/{applicationId}/versions/{semanticVersion}", args)

@inline create_application_version(args) = create_application_version(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.create_cloud_formation_change_set
    create_cloud_formation_change_set([::AWSConfig], arguments::Dict)
    create_cloud_formation_change_set([::AWSConfig]; applicationId=, stackName=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "POST", "/applications/{applicationId}/changesets", arguments::Dict)
    serverlessrepo([::AWSConfig], "POST", "/applications/{applicationId}/changesets", applicationId=, stackName=, <keyword arguments>)

# CreateCloudFormationChangeSet Operation

Creates an AWS CloudFormation change set for the given application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `parameterOverrides = [[ ... ], ...]`
A list of parameter values for the parameters of the application.
```
 parameterOverrides = [[
        "name" => <required> ::String,
        "value" => <required> ::String
    ], ...]
```

## `semanticVersion = ::String`
The semantic version of the application:

[https://semver.org/](https://semver.org/)

## `stackName = ::String` -- *Required*
The name or the unique ID of the stack for which you are creating a change set. AWS CloudFormation generates the change set by comparing this stack's information with the information that you submit, such as a modified template or different parameter input values.

Constraints: Minimum length of 1.

Pattern: ([a-zA-Z][-a-zA-Z0-9]*)|(arn:\\b(aws|aws-us-gov|aws-cn)\\b:[-a-zA-Z0-9:/._+]*)

# Returns

`CreateCloudFormationChangeSetResponse`

# Exceptions

`TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/CreateCloudFormationChangeSet)
"""
@inline create_cloud_formation_change_set(aws::AWSConfig=default_aws_config(); args...) = create_cloud_formation_change_set(aws, args)

@inline create_cloud_formation_change_set(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "POST", "/applications/{applicationId}/changesets", args)

@inline create_cloud_formation_change_set(args) = create_cloud_formation_change_set(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.delete_application
    delete_application([::AWSConfig], arguments::Dict)
    delete_application([::AWSConfig]; applicationId=)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "DELETE", "/applications/{applicationId}", arguments::Dict)
    serverlessrepo([::AWSConfig], "DELETE", "/applications/{applicationId}", applicationId=)

# DeleteApplication Operation

Deletes the specified application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

# Exceptions

`BadRequestException`, `InternalServerErrorException`, `ForbiddenException`, `NotFoundException`, `TooManyRequestsException` or `ConflictException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/DeleteApplication)
"""
@inline delete_application(aws::AWSConfig=default_aws_config(); args...) = delete_application(aws, args)

@inline delete_application(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "DELETE", "/applications/{applicationId}", args)

@inline delete_application(args) = delete_application(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.get_application
    get_application([::AWSConfig], arguments::Dict)
    get_application([::AWSConfig]; applicationId=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}", arguments::Dict)
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}", applicationId=, <keyword arguments>)

# GetApplication Operation

Gets the specified application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `semanticVersion = ::String`
The semantic version of the application to get.

# Returns

`GetApplicationResponse`

# Exceptions

`NotFoundException`, `TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/GetApplication)
"""
@inline get_application(aws::AWSConfig=default_aws_config(); args...) = get_application(aws, args)

@inline get_application(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "GET", "/applications/{applicationId}", args)

@inline get_application(args) = get_application(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.get_application_policy
    get_application_policy([::AWSConfig], arguments::Dict)
    get_application_policy([::AWSConfig]; applicationId=)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}/policy", arguments::Dict)
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}/policy", applicationId=)

# GetApplicationPolicy Operation

Retrieves the policy for the application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

# Returns

`GetApplicationPolicyResponse`

# Exceptions

`NotFoundException`, `TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/GetApplicationPolicy)
"""
@inline get_application_policy(aws::AWSConfig=default_aws_config(); args...) = get_application_policy(aws, args)

@inline get_application_policy(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "GET", "/applications/{applicationId}/policy", args)

@inline get_application_policy(args) = get_application_policy(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.list_application_versions
    list_application_versions([::AWSConfig], arguments::Dict)
    list_application_versions([::AWSConfig]; applicationId=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}/versions", arguments::Dict)
    serverlessrepo([::AWSConfig], "GET", "/applications/{applicationId}/versions", applicationId=, <keyword arguments>)

# ListApplicationVersions Operation

Lists versions for the specified application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `maxItems = ::Int`
The total number of items to return.

## `nextToken = ::String`
A token to specify where to start paginating.

# Returns

`ListApplicationVersionsResponse`

# Exceptions

`NotFoundException`, `TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/ListApplicationVersions)
"""
@inline list_application_versions(aws::AWSConfig=default_aws_config(); args...) = list_application_versions(aws, args)

@inline list_application_versions(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "GET", "/applications/{applicationId}/versions", args)

@inline list_application_versions(args) = list_application_versions(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.list_applications
    list_applications([::AWSConfig], arguments::Dict)
    list_applications([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "GET", "/applications", arguments::Dict)
    serverlessrepo([::AWSConfig], "GET", "/applications", <keyword arguments>)

# ListApplications Operation

Lists applications owned by the requester.

# Arguments

## `maxItems = ::Int`
The total number of items to return.

## `nextToken = ::String`
A token to specify where to start paginating.

# Returns

`ListApplicationsResponse`

# Exceptions

`NotFoundException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/ListApplications)
"""
@inline list_applications(aws::AWSConfig=default_aws_config(); args...) = list_applications(aws, args)

@inline list_applications(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "GET", "/applications", args)

@inline list_applications(args) = list_applications(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.put_application_policy
    put_application_policy([::AWSConfig], arguments::Dict)
    put_application_policy([::AWSConfig]; applicationId=, statements=)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "PUT", "/applications/{applicationId}/policy", arguments::Dict)
    serverlessrepo([::AWSConfig], "PUT", "/applications/{applicationId}/policy", applicationId=, statements=)

# PutApplicationPolicy Operation

Sets the permission policy for an application. See [Application Permissions](https://docs.aws.amazon.com/serverlessrepo/latest/devguide/access-control-resource-based.html#application-permissions) for the list of supported actions that can be used with this operation.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `statements = [[ ... ], ...]` -- *Required*
An array of policy statements applied to the application.
```
 statements = [[
        "actions" => <required> [::String, ...],
        "principals" => <required> [::String, ...],
        "statementId" =>  ::String
    ], ...]
```

# Returns

`PutApplicationPolicyResponse`

# Exceptions

`NotFoundException`, `TooManyRequestsException`, `BadRequestException`, `InternalServerErrorException` or `ForbiddenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/PutApplicationPolicy)
"""
@inline put_application_policy(aws::AWSConfig=default_aws_config(); args...) = put_application_policy(aws, args)

@inline put_application_policy(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "PUT", "/applications/{applicationId}/policy", args)

@inline put_application_policy(args) = put_application_policy(default_aws_config(), args)

"""
    using AWSSDK.ServerlessApplicationRepository.update_application
    update_application([::AWSConfig], arguments::Dict)
    update_application([::AWSConfig]; applicationId=, <keyword arguments>)

    using AWSCore.Services.serverlessrepo
    serverlessrepo([::AWSConfig], "PATCH", "/applications/{applicationId}", arguments::Dict)
    serverlessrepo([::AWSConfig], "PATCH", "/applications/{applicationId}", applicationId=, <keyword arguments>)

# UpdateApplication Operation

Updates the specified application.

# Arguments

## `applicationId = ::String` -- *Required*
The Amazon Resource Name (ARN) of the application.

## `author = ::String`
The name of the author publishing the app.

Minimum length=1. Maximum length=127.

Pattern "^[a-z0-9](([a-z0-9]|-(?!-))*[a-z0-9])?\$";

## `description = ::String`
The description of the application.

Minimum length=1. Maximum length=256

## `homePageUrl = ::String`
A URL with more information about the application, for example the location of your GitHub repository for the application.

## `labels = [::String, ...]`
Labels to improve discovery of apps in search results.

Minimum length=1. Maximum length=127. Maximum number of labels: 10

Pattern: "^[a-zA-Z0-9+\\\\-_:\\\\/@]+\$";

## `readmeBody = ::String`
A text readme file in Markdown language that contains a more detailed description of the application and how it works.

Maximum size 5 MB

## `readmeUrl = ::String`
A link to the readme file in Markdown language that contains a more detailed description of the application and how it works.

Maximum size 5 MB

# Returns

`UpdateApplicationResponse`

# Exceptions

`BadRequestException`, `InternalServerErrorException`, `ForbiddenException`, `NotFoundException`, `TooManyRequestsException` or `ConflictException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/serverlessrepo-2017-09-08/UpdateApplication)
"""
@inline update_application(aws::AWSConfig=default_aws_config(); args...) = update_application(aws, args)

@inline update_application(aws::AWSConfig, args) = AWSCore.Services.serverlessrepo(aws, "PATCH", "/applications/{applicationId}", args)

@inline update_application(args) = update_application(default_aws_config(), args)

end # module ServerlessApplicationRepository

#==============================================================================#
# End of file
#==============================================================================#


text:
using ASE, JuLIP, Test, Printf, DataDeps
using JuLIP.Testing

h0("   JuLIP Tests   ")

@info("preparing the tests...")

verbose=true
isCI = haskey(ENV, "CI")
notCI = !isCI

## ------ loading some prototype potentials via DataDeps

register(DataDep(
    "JuLIP_testpots",
    "A few EAM potentials for testing",
    "https://www.dropbox.com/s/leub1c9ft1mm9fg/JuLIP_data.zip?dl=1",
    post_fetch_method = file -> run(`unzip $file`)
    ))

test_pots = joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"

## ===== some prototype potentials ======
@info("Loading some interatomic potentials . .")
eam_Fe = JuLIP.Potentials.EAM(test_pots * "pfe.plt",
                              test_pots * "ffe.plt",
                              test_pots * "F_fe.plt")
print(" .")
# eam_W = JuLIP.Potentials.FinnisSinclair(
#       test_pots * "W-pair-Wang-2014.plt", test_pots * "W-e-dens-Wang-2014.plt")
# print(" .")
eam_PdAgH = JuLIP.Potentials.eam_from_ase(test_pots * "PdAgH_HybridPd3Ag.eam.alloy")
print(".")
eam_W4 = JuLIP.Potentials.EAM(test_pots * "w_eam4.fs")
println(" done.")

##

julip_tests = [
   ("testaux.jl", "Miscellaneous"),
   ("test_atoms.jl", "Atoms"),
   ("test_build.jl", "Build"),
   ("test_fio.jl", "File IO"),
   ("testanalyticpotential.jl", "Analytic Potential"),
   ("testpotentials.jl", "Potentials"),
   ("test_ad.jl", "AD Potentials"),
   ("testvarcell.jl", "Variable Cell"),
   ("testhessian.jl", "Hessian"),
   ("test_onebody.jl", "One-Body"),
   ("test_eam.jl", "EAM"),
   ]

# add solver tests if not on travis
if !isCI
   push!(julip_tests, ("testsolve.jl", "Solve"))
else
   @info("on CI : don't run solver tests")
end

# "testexpvarcell.jl";  # USE THIS TO WORK ON EXPCELL IMPLEMENTATION

@testset "JuLIP" begin
   for (testfile, testid) in julip_tests
      h1("Testset $(testid)")
      @testset "$(testid)" begin include(testfile); end
   end
end

# TODO:
# - fix performance regression in morseold / testanalyticpotential
# - some other analytic EAM potentials
# - finish proper EAM implementation
#   do this via ASE.jl and importing the data from ASE!!
# - maybe simplify the hessian assembly by skipping the step via
#   hessian_pos


text:
using LocalMemoize

abstract type Camera end

struct PerspectiveCamera <: Camera
    eye::Vector3
    front::Vector3
    ref_up::Vector3
    fov::Float64
end
get_right(pc::PerspectiveCamera) = pc.front × pc.ref_up
get_up(pc::PerspectiveCamera) = get_right(pc) × pc.front
fov_scale(pc::PerspectiveCamera) = tan(pc.fov / 2) * 2

@memoize function generate_ray(@memo(pc::PerspectiveCamera), x::Real, y::Real)
    @assert 0 <= x <= 1
    @assert 0 <= y <= 1
    @memo 𝐫 = get_right(pc) * fov_scale(pc)
    @memo 𝐮 = get_up(pc) * fov_scale(pc)
    Ray3(pc.eye, normalize(pc.front + 𝐫 * (x - 0.5) + 𝐮 * (y - 0.5)))
end


text:
##
#
# Generate all non-zero roots of polynomials with coefficients 1 or -1 of
# degree at most 25.
#
# Warning: this may take a while (a few hours in my case) and generate almost
#          50 GB file.
#

using ChristensenSet

polynomials = PolynomialIterator([-1.0, 1.0], 25)
image = RootsImage(-2.0-2.0im, 2.0+2.0im, 2000, 2000)
image_names = n -> "animation/img_$(lpad(string(n), 2, "0")).png"
DATA_FILE = "ones_Float64.dat"

if !isdir("animation")
  mkdir("animation")
end

@time animate_growth(polynomials, image, DATA_FILE, output_filename=image_names)


text:
@testset begin "Infos database"
    db = SQLite.DB("test/data/TEST_PowerFlowNetworks_SQLite.sqlite")
    SQLite.drop!(db, "instances")
    @test !has_opf_tables(db)
    db = setup_db("test/data/TEST_PowerFlowNetworks_SQLite.sqlite"; delete_if_exists=true)
    @test table_count(db, "instances") == 0
    @test table_count(db, "decompositions") == 0
    @test table_count(db, "merges") == 0
    @test table_count(db, "combinations") == 0
    @test has_opf_tables(db)
    DBInterface.execute(db, "DELETE FROM instances")
    states = count_missing_columns(db, "instances", ["name", "scenario", "nbus", "nbranch", "nbranch_unique", "ngen"])
    @test collect(values(states)) == zeros(Int, length(states))

    name = "case6ww"
    scenario = 0
    source_path = "./test/data/case6ww.m"
    source_type = "MATPOWERM"
    date = Dates.now()
    load_instance_in_db!(db, name, scenario, source_path, source_type, date)
    @test table_count(db, "instances") == 1
    @test table_count(db, "decompositions") == 0
    @test table_count(db, "merges") == 0
    @test table_count(db, "combinations") == 0

    states = count_missing_columns(db, "instances", ["name", "scenario", "nbus", "nbranch", "nbranch_unique", "ngen"])
    @test states == Dict("name" => 0,
                         "scenario" => 0,
                         "nbus" => 1,
                         "nbus" => 1,
                         "nbranch" => 1,
                         "nbranch_unique" => 1,
                         "ngen" => 1,
                         "total" => 1)

    save_basic_features_instances!(db)
    states = count_missing_columns(db, "instances", ["name", "scenario", "nbus", "nbranch", "nbranch_unique", "ngen"])
    @test collect(values(states)) == zeros(Int, length(states))

    states = count_missing_columns(db, "instances", ["global_clustering_coefficient"])
    @test states["global_clustering_coefficient"] == 1
    @test states["total"] == 1
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
    compute_loss_weights(data_presample::SubArray{Union{Float64,Missing}}, n::Int64, standardise_error::Bool, weights::Nothing)
    compute_loss_weights(data_presample::SubArray{Union{Float64,Missing}}, n::Int64, standardise_error::Bool, weights::FloatVector)

Compute weights for the forecast error.
"""
compute_loss_weights(data_presample::SubArray{Union{Float64,Missing}}, n::Int64, standardise_error::Bool, weights::Nothing) = standardise_error ? std_skipmissing(data_presample).^2 : ones(n);
compute_loss_weights(data_presample::SubArray{Union{Float64,Missing}}, n::Int64, standardise_error::Bool, weights::FloatVector) = standardise_error ? std_skipmissing(data_presample).^2 : weights;

"""
    compute_loss(weighted_resid::AbstractArray{Float64})
    compute_loss(weighted_resid::AbstractArray{Missing})
    compute_loss(weighted_resid::AbstractArray{Union{Float64, Missing}})

Compute loss function.
"""
compute_loss(weighted_resid::AbstractArray{Float64}) = [mean(mean(weighted_resid, dims=1), dims=2)[1], 0.0];
compute_loss(weighted_resid::AbstractArray{Missing}) = [0.0, 1.0];
function compute_loss(weighted_resid::AbstractArray{Union{Float64, Missing}})
    loss = 0.0;
    inactive_periods = 0.0;
    T = size(weighted_resid,2);

    for t=1:T
        weighted_resid_t = @view weighted_resid[:,t];
        if sum(.~ismissing.(weighted_resid_t)) > 0
            loss += mean_skipmissing(weighted_resid_t);
        else
            inactive_periods += 1.0;
        end
    end

    if inactive_periods == T
        return [0.0, 1.0];
    else
        return [loss/(T-inactive_periods), 0.0];
    end
end

"""
    fc_err(validation_settings::ValidationSettings, p::Int64, λ::Number, α::Number, β::Number)

Compute the in-sample / out-of-sample error associated with the candidate hyperparameters

# Arguments
- `validation_settings`: ValidationSettings struct
- `p`: (candidate) number of lags in the vector autoregression
- `λ`: (candidate) overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: (candidate) weight associated to the LASSO component of the elastic-net penalty
- `β`: (candidate) additional shrinkage for distant lags (p>1)
- `jth_jackknife_data`: j-th jackknife sample (default: nothing)

# References
Pellegrino (2020)
"""
function fc_err(validation_settings::ValidationSettings, p::Int64, λ::Number, α::Number, β::Number; jth_jackknife_data::Union{JArray{Float64}, Nothing}=nothing)

    t0 = validation_settings.T;

    # Out-of-sample error
    if validation_settings.err_type != 1
        t0 = validation_settings.t0;
    end

    # Jackknife out-of-sample
    if validation_settings.err_type > 2
        data = jth_jackknife_data;

    # Standard in-sample or out-of-sample
    else
        data = validation_settings.Y;
    end

    # Data
    data_presample = @view data[:, 1:t0];
    μ = mean_skipmissing(data_presample);
    Y = @. data_presample - μ;
    Y_output = @. data - μ;

    # Construct estim_settings
    estim_settings = EstimSettings(Y, Y_output, p, λ, α, β, ε=validation_settings.ε, tol=validation_settings.tol, max_iter=validation_settings.max_iter, prerun=validation_settings.prerun, verb=validation_settings.verb_estim);

    # Estimate penalised VAR
    kalman_settings = ecm(estim_settings);

    # Run Kalman filter
    status = KalmanStatus();
    for t=1:kalman_settings.T
        kfilter!(kalman_settings, status);
    end

    # Residuals
    forecast  = hcat(status.history_X_prior...)[1:validation_settings.n, :];

    # Compute weights
    w = compute_loss_weights(data_presample, validation_settings.n, validation_settings.standardise_error, validation_settings.weights);

    # Weighted error
    weighted_resid = @. w*(forecast - Y_output)^2;

    # In-sample error
    if validation_settings.err_type == 1
        return compute_loss(weighted_resid);

    # Out-of-sample error
    else
        weighted_resid_oos = @view weighted_resid[:, t0+1:end];
        return compute_loss(weighted_resid_oos);
    end
end

"""
    jackknife_err(validation_settings::ValidationSettings, jackknife_data::JArray{Float64, 3}, p::Int64, λ::Number, α::Number, β::Number)

Return the jackknife out-of-sample error.

# Arguments
- `validation_settings`: ValidationSettings struct
- `jackknife_data`: jackknife partitions
- `p`: (candidate) number of lags in the vector autoregression
- `λ`: (candidate) overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: (candidate) weight associated to the LASSO component of the elastic-net penalty
- `β`: (candidate) additional shrinkage for distant lags (p>1)

# References
Pellegrino (2020)
"""
function jackknife_err(validation_settings::ValidationSettings, jackknife_data::JArray{Float64, 3}, p::Int64, λ::Number, α::Number, β::Number)

    # Error management
    if validation_settings.err_type <= 2
        error("Wrong err_type for jackknife_err!");
    end

    # Number of jackknife samples
    samples = size(jackknife_data, 3);

    # Compute jackknife loss
    verb_message(validation_settings.verb_estim, "jackknife_err > running $samples iterations on $(nworkers()) workers");

    output_fc_err = @sync @distributed (+) for j=1:samples
        fc_err(validation_settings, p, λ, α, β, jth_jackknife_data=jackknife_data[:,:,j]);
    end

    # Compute average jackknife loss
    loss, inactive_samples = output_fc_err;
    if samples == inactive_samples
        error("All samples are inactive! Check the initial settings or try a different random seed.");
    end
    loss *= 1/(samples-inactive_samples);

    # Return output
    return loss;
end

"""
    select_hyperparameters(validation_settings::ValidationSettings, γ_grid::HyperGrid)

Select the tuning hyper-parameters for the elastic-net vector autoregression.

# Arguments
- `validation_settings`: ValidationSettings struct
- `γ_grid`: HyperGrid struct

# References
Pellegrino (2020)
"""
function select_hyperparameters(validation_settings::ValidationSettings, γ_grid::HyperGrid)

    # Check inputs
    check_bounds(validation_settings.max_iter, 3);
    check_bounds(validation_settings.max_iter, validation_settings.prerun);
    check_bounds(validation_settings.n, 2); # It supports only multivariate models (for now ...)

    if length(γ_grid.p) != 2 || length(γ_grid.λ) != 2 || length(γ_grid.α) != 2 || length(γ_grid.β) != 2
        error("The grids include more than two entries. They must include only the lower and upper bounds for the grids!")
    end

    check_bounds(γ_grid.p[2], γ_grid.p[1]);
    check_bounds(γ_grid.λ[2], γ_grid.λ[1]);
    check_bounds(γ_grid.α[2], γ_grid.α[1]);
    check_bounds(γ_grid.β[2], γ_grid.β[1]);
    check_bounds(γ_grid.p[1], 1);
    check_bounds(γ_grid.λ[1], 0);
    check_bounds(γ_grid.α[1], 0, 1);
    check_bounds(γ_grid.α[2], 0, 1);
    check_bounds(γ_grid.β[1], 1);

    # Construct grid of hyperparameters - random search algorithm
    errors     = zeros(γ_grid.draws);
    candidates = zeros(4, γ_grid.draws);

    for draw=1:γ_grid.draws
        candidates[:,draw] = [rand(γ_grid.p[1]:γ_grid.p[2]),
                              rand(Uniform(γ_grid.λ[1], γ_grid.λ[2])),
                              rand(Uniform(γ_grid.α[1], γ_grid.α[2])),
                              rand(Uniform(γ_grid.β[1], γ_grid.β[2]))];
    end

    if ~isnothing(validation_settings.log_folder_path)
        io = open("$(validation_settings.log_folder_path)/status.txt", "w+");
        global_logger(ConsoleLogger(io));
    end

    # Generate partitions for the block jackknife out-of-sample
    if validation_settings.err_type == 3
        jackknife_data = block_jackknife(validation_settings.Y, validation_settings.subsample);

    # Generate partitions for the artificial jackknife
    elseif validation_settings.err_type == 4
        jackknife_data = artificial_jackknife(validation_settings.Y, validation_settings.subsample, validation_settings.max_samples);
    end

    for iter=1:γ_grid.draws

        # Retrieve candidate hyperparameters
        p, λ, α, β = candidates[:,iter];
        p = Int64(p);

        # Update log
        if validation_settings.verb == true
            @info "$(round(now(), Dates.Second(1))) select_hyperparameters (error estimator $(validation_settings.err_type)) > running iteration $iter (out of $(γ_grid.draws)), γ=($(round(p,digits=3)), $(round(λ,digits=3)), $(round(α,digits=3)), $(round(β,digits=3)))";
            if ~isnothing(validation_settings.log_folder_path)
                flush(io);
            end
        end

        #=
        Some candidate hyperparameters the estimation of the VAR can be unstable. This generally happens when:
        - The candidates are extreme
        - validation_settings.subsample is high and there is not enough shrinkage

        A priori, there is not a simple way to construct a grid of candidates that does not result in errors in the ECM.
        The try-catch statement below handles this problem by skipping the candidate values that generate model instability.

        Note that this issue is more likely to happen with the block-jackknife.
        =#

        try

            # In-sample or standard out-of-sample
            if validation_settings.err_type <= 2
                errors[iter], inactive_sample = fc_err(validation_settings, p, λ, α, β);
                if inactive_sample == 1
                    error("The validation sample is a matrix of missings!");
                end

            # Jackknife out-of-sample
            else
                errors[iter] = jackknife_err(validation_settings, jackknife_data, p, λ, α, β);
            end

        catch error_iter

            # Extract info on the error
            error_ex, error_msg, error_stacktrace_array = error_info(error_iter);
            error_stacktrace = join([*(string(error_stacktrace_array[error_line]), "\n") for error_line=1:length(error_stacktrace_array)]);

            # Update log
            @error "$(round(now(), Dates.Second(1))) $error_msg \n $error_stacktrace";

            # The instability pops up in the update_loglik! function
            if isa(error_ex, DomainError) && occursin("logdet", "$error_stacktrace") && occursin("update_loglik!", "$error_stacktrace")
                errors[iter] = Inf;

            # Any other error
            else
                rethrow(error_iter);
            end
        end
    end

    # Return output
    return candidates, errors;
end


text:
metadata_packages = [
     "Combinatorics"]

Pkg.init()
Pkg.update()

for package=metadata_packages
    Pkg.add(package)
    # Precompile all packages, so scripts can operate in a read-only container
    Base.compilecache(package)
end

max_stars_repo_name: cmplopes/alpine-julia

max_issues_repo_name: cmplopes/alpine-julia

max_forks_repo_name: cmplopes/alpine-julia


text:
using PyPlot, PyCall
using Distributions
using StatsFuns
using Requests
PyDict(matplotlib["rcParams"])["mathtext.fontset"] = "cm"
PyDict(matplotlib["rcParams"])["mathtext.rm"] = "serif"
PyDict(matplotlib["rcParams"])["lines.linewidth"] = 1.5
PyDict(matplotlib["rcParams"])["font.family"] = "TakaoPGothic"

test = false

function VI(C, a1, b1, a2, b2, phi, max_iter)
    N = size(C, 1)
    phi_h = deepcopy(phi)

    a1_h = NaN
    b1_h = NaN
    a2_h = NaN
    b2_h = NaN
    for i in 1 : max_iter
        # update parameters
        phi_sum = sum(phi_h)

        #delta1 = reverse(cumsum(reverse(phi_h)))
        delta1 = reverse(cumsum(reverse(cat(1, phi_h[2:end], 0))))
        #tmp = cumsum(phi_h)
        #delta2 = cat(1, 0, tmp[1:end-1])
        delta2 = cumsum(phi_h)
        a1_h = a1 + sum(C.*delta1)
        b1_h = b1 + sum(delta1)
        a2_h = a2 + sum(C.*delta2)
        b2_h = b2 + sum(delta2)

        # update tau
        ln_p1 = Array{Float64, 1}(N)
        ln_p2 = Array{Float64, 1}(N)
        for n in 1 : N
            ln_p1[n] = C[n] * (digamma(a1_h) - log(b1_h)) - a1_h/b1_h
            ln_p2[n] = C[n] * (digamma(a2_h) - log(b2_h)) - a2_h/b2_h
        end

        ln_tmp = Array{Float64, 1}(N)
        for n in 1 : N
            #ln_tmp[n] = sum(ln_p1[1:n]) + sum(ln_p2[n+1:end]) + log(phi[n])
            ln_tmp[n] = sum(ln_p1[1:n-1]) + sum(ln_p2[n:end]) + log(phi[n])
        end
        phi_h = exp(ln_tmp - logsumexp(ln_tmp))
    end
    return phi_h, a1_h / b1_h, a2_h / b2_h
end

# model params
a1 = 2.0
b1 = 0.1
a2 = 2.0
b2 = 0.1

# generate or download data
if test
    N = 74
    phi = ones(N) / N

    # generate sample
    lambda1 = rand(Gamma(a1, 1.0/b1))
    lambda2 = rand(Gamma(a2, 1.0/b2))
    tau = rand(Categorical(phi))
    C = Array{Float64, 1}(N)
    for n in 1 : tau
        C[n] = rand(Poisson(lambda1))
    end
    for n in tau+1 : N
        C[n] = rand(Poisson(lambda2))
    end
else
    data_url = "https://git.io/vXTVC"
    dir_name = "data"
    file_name = "txtdata.csv"
    if !isdir(dir_name)
        mkdir(dir_name)
    end
    if isfile(dir_name * "/" * file_name)
        println("read csv file $(dir_name * "/" * file_name)")
        data = readcsv(dir_name * "/" * file_name)
    else
        println("download csv file from $data_url")
        res = get(data_url)
        data = readcsv(IOBuffer(res.data))
        writecsv(dir_name * "/" * file_name, data)
    end
    C = data[:,1]

    N = size(C, 1)
    phi = ones(N) / N
end

# visualize
if test
    figure("data samples")
    clf()
    bar(1:N, C, color="#348ABD")
    plot([1, tau], [lambda1, lambda1], "g--")
    plot([tau+1, N], [lambda2, lambda2], "g--")
    plot([tau+0.5, tau+0.5],[0, maximum(C)],"r--")
    xlim([0, N])
else
    figure("data samples")
    clf()
    bar(1:N, C, color="#348ABD")
    xlim([0, N])
end

# inference
max_iter = 100
phi_est, lambda1_est, lambda2_est = VI(C, a1, b1, a2, b2, phi, max_iter)

# calc expectation
tau_expt = sum(phi_est .* range(1,N))
lambda_est = Array{Float64, 1}(N)
for n in 1 : N
    tmp = sum(phi_est[1:n])
    lambda_est[n] =(1.0 - tmp) * lambda1_est + tmp * lambda2_est
end

figure("results")
clf()
subplot(2,1,1)
bar(1:N, C, color="#348ABD", label="受信数")
plot(1:N, lambda_est, color="#E24A33", label="期待値")
xlim([0, N])
ylabel("受信メッセージ数")
xlabel("経過日数")
legend()
subplot(2,1,2)
bar(1:N, phi_est, color="#348ABD", label="事後確率")
plot([tau_expt, tau_expt], [0, maximum(phi_est)], "r--", label="期待値")
xlim([0, N])
ylabel("事後確率")
xlabel("経過日数")
legend()

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
module List
import Base: length, (:), map, iterate, reverse, filter
export list, cons, nil, linkedlist, (^)

# if generic:
# map(f, l :: nil{T} ) where T  = nil()

abstract type linkedlist end

struct cons <: linkedlist
    head
    tail :: linkedlist
end

struct nil <: linkedlist end

length(l :: nil) = 0
length(l :: cons) = 1 + length(l.tail)

head(l :: nil )   = nothing
head(l :: cons)   = l.head

tail(l :: nil )   = nothing
tail(l :: cons)   = l.tail

iter(f, l :: nil )  = nothing
iter(f, l :: cons)  = begin f(l.head); iter(f, l.tail) end

map(f, l :: nil )   = nil()
map(f, l :: cons)   = cons(f(l.head), map(f, l.tail))

(^)(hd, tl :: linkedlist) = cons(hd, tl)

function filter(f :: Function, l :: linkedlist)
    l2 = nil()
    for h in l
        if f(h)
          l2 = cons(h, l2)
        end
    end
    reverse(l2)
end

iterate(l :: linkedlist, ::nil) = nothing
function iterate(l :: linkedlist, state::cons = l)
    state.head, state.tail
end

function reverse(l :: linkedlist)
    l2 = nil()
    for hd in l
        l2 = cons(hd, l2)
    end
    l2
end

function list(xs...)
    l = nil()
    for i=length(xs):-1:1
        l = cons(xs[i], l)
    end
    l
end

end


text:
# ---
# title: 978. Longest Turbulent Subarray
# id: problem978
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming, Sliding Window
# link: <https://leetcode.com/problems/longest-turbulent-subarray/description/>
# hidden: true
# ---
#
# Given an integer array `arr`, return _the length of a maximum size turbulent
# subarray of_ `arr`.
#
# A subarray is **turbulent** if the comparison sign flips between each adjacent
# pair of elements in the subarray.
#
# More formally, a subarray `[arr[i], arr[i + 1], ..., arr[j]]` of `arr` is said
# to be turbulent if and only if:
#
#   * For `i <= k < j`:
#     * `arr[k] > arr[k + 1]` when `k` is odd, and
#     * `arr[k] < arr[k + 1]` when `k` is even.
#   * Or, for `i <= k < j`:
#     * `arr[k] > arr[k + 1]` when `k` is even, and
#     * `arr[k] < arr[k + 1]` when `k` is odd.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [9,4,2,10,7,8,8,1,9]
#     Output: 5
#     Explanation: arr[1] > arr[2] < arr[3] > arr[4] < arr[5]
#
#
# **Example 2:**
#
#
#
#     Input: arr = [4,8,12,16]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: arr = [100]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 4 * 104`
#   * `0 <= arr[i] <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
using DiffEqSensitivity, OrdinaryDiffEq, DiffEqCallbacks, DiffEqFlux, Flux
using Random, Test

function test_hybridNODE(sensealg)
    Random.seed!(12345)
    datalength = 100
    tspan = (0.0,100.0)
    t = range(tspan[1],tspan[2],length=datalength)
    target = 3.0*(1:datalength)./datalength  # some dummy data to fit to
    cbinput = rand(1, datalength) #some external ODE contribution
    pmodel = Chain(
        Dense(2, 10, initW=zeros),
        Dense(10, 2, initW=zeros))
    p, re = Flux.destructure(pmodel)
    dudt(u,p,t) = re(p)(u)

    # callback changes the first component of the solution every time
    # t is an integer
    function affect!(integrator, cbinput)
        event_index = round(Int,integrator.t)
        integrator.u[1] += 0.2*cbinput[event_index]
    end
    callback = PresetTimeCallback(collect(1:datalength),(int)->affect!(int, cbinput))

    # ODE with Callback
    prob = ODEProblem(dudt,[0.0, 1.0],tspan,p)

    function predict_n_ode(p)
        arr = Array(solve(prob, Tsit5(),
            p=p, sensealg=sensealg, saveat=2.0, callback=callback))[1,2:2:end]
        return arr[1:datalength]
    end

    function loss_n_ode()
        pred = predict_n_ode(p)
        loss = sum(abs2,target .- pred)./datalength
    end

    cb = function () #callback function to observe training
        pred = predict_n_ode(p)
        display(loss_n_ode())
    end
    @show sensealg
    Flux.train!(loss_n_ode, Flux.params(p), Iterators.repeated((), 20), ADAM(0.005), cb = cb)
    @test loss_n_ode() < 0.4
    println("  ")
end

function test_hybridNODE2(sensealg)
    Random.seed!(12345)
    u0 = Float32[2.; 0.; 0.; 0.]
    tspan = (0f0,1f0)

    ## Get goal data
    function trueaffect!(integrator)
        integrator.u[3:4] = -3*integrator.u[1:2]
    end
    function trueODEfunc(dx,x,p,t)
        @views dx[1:2] .= x[1:2] + x[3:4]
        dx[1] += x[2]
        dx[2] += x[1]
        dx[3:4] .= 0f0
    end
    cb_ = PeriodicCallback(trueaffect!,0.1f0,save_positions=(true,true),initial_affect=true)
    prob = ODEProblem(trueODEfunc,u0,tspan)
    sol = solve(prob,Tsit5(),callback=cb_,save_everystep=false,save_start=true)
    ode_data = Array(sol)[1:2,1:end]'

    ## Make model
    dudt2 = Chain(Dense(4,50,tanh),
                    Dense(50,2))
    p,re = Flux.destructure(dudt2) # use this p as the initial condition!
    function affect!(integrator)
        integrator.u[3:4] = -3*integrator.u[1:2]
    end
    function ODEfunc(dx,x,p,t)
        dx[1:2] .= re(p)(x)
        dx[3:4] .= 0f0
    end
    z0 = u0
    prob = ODEProblem(ODEfunc,z0,tspan)
    cb = PeriodicCallback(affect!,0.1f0,save_positions=(true,true),initial_affect=true)

    ## Initialize learning functions
    function predict_n_ode()
        _prob = remake(prob,p=p)
        Array(solve(_prob,Tsit5(),u0=z0,p=p,callback=cb,save_everystep=false,save_start=true,sensealg=sensealg))[1:2,:]
    end
    function loss_n_ode()
        pred = predict_n_ode()[1:2,1:end]'
        loss = sum(abs2,ode_data .- pred)
        loss
    end
    loss_n_ode() # n_ode.p stores the initial parameters of the neural ODE
    cba = function ()  #callback function to observe training
        pred = predict_n_ode()[1:2,1:end]'
        display(sum(abs2,ode_data .- pred))
        return false
    end
    cba()

    ## Learn
    ps = Flux.params(p)
    data = Iterators.repeated((), 25)

    @show sensealg

    Flux.train!(loss_n_ode, ps, data, ADAM(0.0025), cb = cba)

    @test loss_n_ode() < 0.5

    println("  ")
end

@testset "PresetTimeCallback" begin
    test_hybridNODE(ForwardDiffSensitivity())
    test_hybridNODE(BacksolveAdjoint())
    test_hybridNODE(InterpolatingAdjoint())
    test_hybridNODE(QuadratureAdjoint())
end

@testset "PeriodicCallback" begin
    test_hybridNODE2(ReverseDiffAdjoint())
    test_hybridNODE2(BacksolveAdjoint())
    test_hybridNODE2(InterpolatingAdjoint())
    test_hybridNODE2(QuadratureAdjoint())
end


text:
"""
    WBS( segment_cost, n[,th_const = 1.3, sigma = 1.0, M = 5000, do_seeded = false, shrink = 1/sqrt(2)])

Runs the Wild Binary segmentation algorithm for the CUSUM cost function `segment_cost` for a time series of length `n`.
    The test threshold is determined by `th_const` and the known or estimated standard deviation `sigma`.
    `do_seeded` determined whether to use seeded intervals. If true, `shrink` is the decay factor for interval length; if not, `M` random intervals are drawn, and returns the position of found changepoints, and
the cost of this segmentation.

See also: [`@BS`](@ref), [`@segment_cost`](@ref)

# Returns
* `result::Array{NamedTuple{s::Int64, e::Int64, cpt::Int64, CUSUM::Float64, min_th::Float64, scale::Int64}}`:
    * `s::Int64`: Start of window in which maximum CUSUM was located
    * `e::Int64`: End of window in which maximum CUSUM was located
    * `cpt::Int64`: Location of maximum CUSUM
    * `CUSUM::Float64`: Value of maximum CUSUM
    * `min_th::Float64`: Threshold used
    * `scale::Int64`: Number of recursion iterations before the change point was detected
# Example

```julia-repl
# Sample Normal time series with changing mean
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
# Run wild binary segmentation
seg_cost = CUSUM(sample, σ)
WBS_return = WBS(seg_cost, n)
```

# References
Fryzlewicz, P. (2014) Wild binary segmentation for multiple change-point detection, Annals of Statistics 42(6), 2243-2281
Kovács, S., Li, H., Bühlmann, P., & Munk, A. (2020). Seeded Binary Segmentation: A general methodology for fast and optimal change point detection. arXiv preprint arXiv:2002.06633.
"""
#import StatsBase
#using StatsBase
# sigma = StatsBase.mad(x)
result_type = @NamedTuple{s::Int64, e::Int64, cpt::Int64, CUSUM::Float64, min_th::Float64, scale::Int64}

function WBS( segment_cost::Function , n::Int64, th_const::Float64 = 1.3, sigma::Float64 = 1.0,
    M::Int64 = 5000, do_seeded::Bool = false, shrink::Float64 = 1/sqrt(2))
    result = Array{result_type}(undef,0)

    # thresulthold
    th = sigma * th_const * sqrt(2 * log(n))
    scale = 0
    # call WBS
    result = WBS_RECUR(segment_cost, 1, n, th, result, scale, M, do_seeded, shrink)

    return result
end

function WBS_RECUR(segment_cost::Function, s::Int64, e::Int64, th::Float64,
    result::Array{result_type}, scale::Int64, M::Int64 = 5000, do_seeded::Bool = false, shrink::Float64 = 1/sqrt(2))
    #k = length(CP) - 2 #current number of change points
    minlength = 20 #minimum interval length
    scale += 1 #update scale

    # draw intervals
    if do_seeded
        tau = get_seeded_intervals(s,e, shrink, minlength)
        M = length(tau)
    else
        tau = Tuple{Int, Int}[] # Segmentations to test
        push!(tau, (s, e))
        for m in 2:M
            s_m = rand(s:(e-minlength))
            e_m = rand( (s_m+minlength):e)
            push!(tau, (s_m, e_m))
        end
    end

    # compute CUSUMs find candidate change point
    x = fill(0.0,M, e-s+1)#zeros(Float64, M, e-s+1)#Array{Float64}(NaN, M, e-s+1 ) #CUSUM array
    for m in 1:M
        a, b = pop!(tau)
        for j in (a + 1):(b - 2)
            x[m,j-s+1] = segment_cost(s,e,j)
        end
    end
    #x = replace(x, 0.0 => maximum(x))
    maxval , ind = findmax(x)

    # segment here
    if maxval > th #&& scale < 10
        cpt = ind[2] + s -1
        newresult = (s=s,e=e, cpt=cpt, CUSUM = maxval, min_th=th, scale=scale )
        push!(result, newresult)
        #push!(CP,cpt )
        #push!(CUSUM, minval)
        if cpt - s > 2*minlength
            result = WBS_RECUR(segment_cost, s, cpt, th, result, scale, M, do_seeded, shrink)
        end
        if e - cpt > 2*minlength
            result = WBS_RECUR(segment_cost, cpt+1, e, th, result,  scale, M, do_seeded, shrink)
        end
    end
    #CP = sort(CP)

    return result #CUSUM, CP
end

"""
    get_WBS_changepoints(segment_cost, object, Kmax = 1,  alpha = 1.01)

Obtains changepoints from an `object` (output by the WBS function) by minimising an information criterion.
`segment_cost` is a cost function, normally sSIC, and `alpha` is the penalty exponent. `Kmax` determines the maximum number of changepoints to detect.

See also: [`@WBS`](@ref), [`@segment_cost`](@ref)

# Returns
* `out_cps`: Vector of changepoint locations under optimal segmentation
* `optimal_number`: Integer number of changepoints which optimises the criterion
* `cost`: Vector of segmentations cost for 0, ..., Kmax  changes
* `object[1:optimal]`: WBS object truncated at optimal number
# Example

```julia-repl
# Sample Normal time series with changing mean
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)

# Run wild binary segmentation
seg_cost_CUSUM = CUSUM(sample, σ)
WBS_return = WBS(seg_cost_CUSUM, n)

# Obtain changepoints
seg_cost_sSIC = sSIC(sample)
get_WBS_changepoints(seg_cost_sSIC, WBS_return, 10)
```

# References
Fryzlewicz, P. (2014) Wild binary segmentation for multiple change-point detection, Annals of Statistics 42(6), 2243-2281
"""

function get_WBS_changepoints(segment_cost::Function, object::Array{NamedTuple{(:s, :e, :cpt, :CUSUM, :min_th, :scale),Tuple{Int64,Int64,Int64,Float64,Float64,Int64}},1}, Kmax::Int64 = 1,  alpha::Float64 = 1.01)
    #result = Array{result_type}(undef,0)
    Kmax = min(Kmax,length(object))
    if length(object) > 0
        cost = zeros(Kmax+1)
        sorted = sort(object, lt = (a, b) -> ( a[4] < b[4]) ) #a[6] < b[6] &&
        cps_from_sorted = Vector{Int64}()
        for i in 1:Kmax
            push!(cps_from_sorted, sorted[i][3])
        end

        cost[1] = segment_cost(1, segment_cost.n)
        for i in 1:Kmax
            current_cps = sort(union([0, segment_cost.n], cps_from_sorted[1:i])) #select cps at current threshold
            for j in 1:(i+1)
                cost[i+1] = cost[i+1] + segment_cost(current_cps[j]+1, current_cps[j+1])
            end
        end

        pen = log(segment_cost.n)^alpha
        penvec = collect(0:Kmax) * pen
        cost += penvec

        optimal_number = findmin(cost)[2] -1

        out_cps = Vector{Int64}()
        for i in 1:optimal_number
            push!(out_cps, sorted[i][3])
        end
        out_cps = sort(out_cps)

        return  out_cps, optimal_number, cost, object[1:optimal_number]
    else
         println("no changepoints detected")
    end

end

function get_seeded_intervals(s::Int64, e::Int64,  shrink::Float64 = 1/sqrt(2), minlength = 10)
    n = e-s+1;
    max_k = Int(ceil(log(1/shrink, float(n) )))
    out = Tuple{Int, Int}[]

    for k = 1:max_k
        n_k = Int(2 * ceil(shrink^(1-k) ) - 1)
        l_k = n*shrink^(k-1)
        s_k = (n-l_k)/(n_k -1)
        if !isnan(s_k) && l_k >= minlength
            for i in 1:n_k
                left =  s + Int( floor((i-1)*s_k) )
                right = s + Int( ceil((i-1)*s_k + l_k) ) -1
                push!(out, (left, right) )
            end
        end
    end
    return out
end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

@time begin
    import SparseArrays
    import Base.Threads: @threads, threadid, SpinLock

    #  ----------------------------------------------------------------------------
    # Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
    # package, then you must activate the package enviroment (see README)
    import Chemostat_Heerden2013
    const ChHd = Chemostat_Heerden2013
    const Hd  = ChHd.HeerdenData;
    const BD  = ChHd.BegData;
    const iJR = ChHd.iJR904

    #  ----------------------------------------------------------------------------
    # run add "https://github.com/josePereiro/Chemostat" in the
    # julia Pkg REPL for installing the package
    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChLP = Ch.LP

    import JuMP, GLPK
    using Serialization
    import UtilsJL
    const UJL = UtilsJL
end

## -----------------------------------------------------------------------------------------------
const FBA_Z_FIX_MAX_VG_MIN_COST = :FBA_Z_FIX_MAX_VG_MIN_COST
const FBA_Z_FIX_MIN_VG_MIN_COST = :FBA_Z_FIX_MIN_VG_MIN_COST
const FBA_Z_FIX_MAX_VG_MAX_COST = :FBA_Z_FIX_MAX_VG_MAX_COST
const FBA_Z_FIX_MIN_VG_MAX_COST = :FBA_Z_FIX_MIN_VG_MAX_COST
const FBA_Z_VG_FIX_MAX_COST = :FBA_Z_VG_FIX_MAX_COST
const FBA_Z_VG_FIX_MIN_COST = :FBA_Z_VG_FIX_MIN_COST
const FBA_Z_FIX_MAX_COST = :FBA_Z_FIX_MAX_COST
const FBA_Z_FIX_MIN_COST = :FBA_Z_FIX_MIN_COST
const FBA_MAX_Z_MIN_COST = :FBA_MAX_Z_MIN_COST
const FBA_MAX_Z_MAX_COST = :FBA_MAX_Z_MAX_COST

const EXPS = 1:13

## -----------------------------------------------------------------------------------------------
# Data container
LPDAT = UJL.DictTree()

## -----------------------------------------------------------------------------------------------
# FBA
let
    objider = iJR.BIOMASS_IDER
    costider = iJR.COST_IDER
    exglcider = iJR.EX_GLC_IDER
    max_sense = -1.0
    min_sense = 1.0

    for exp in EXPS

        @info("Doing ", exp); println()

        # FBA_Z_FIX_MIN_VG_MIN_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # min vg
            fbaout1 = ChLP.fba(model, exglcider; sense = min_sense)
            exglc = ChU.av(model, fbaout1, exglcider)
            ChU.bounds!(model, exglcider, exglc, exglc)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = min_sense)

            LPDAT[FBA_Z_FIX_MIN_VG_MIN_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MIN_VG_MIN_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_FIX_MIN_VG_MAX_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # min vg
            fbaout1 = ChLP.fba(model, exglcider; sense = min_sense)
            exglc = ChU.av(model, fbaout1, exglcider)
            ChU.bounds!(model, exglcider, exglc, exglc)

            # max cost
            fbaout = ChLP.fba(model, costider; sense = max_sense)

            LPDAT[FBA_Z_FIX_MIN_VG_MAX_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MIN_VG_MAX_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_FIX_MAX_VG_MIN_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # max vg
            fbaout1 = ChLP.fba(model, exglcider; sense = max_sense)
            exglc = ChU.av(model, fbaout1, exglcider)
            ChU.bounds!(model, exglcider, exglc, exglc)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = min_sense)

            LPDAT[FBA_Z_FIX_MAX_VG_MIN_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MAX_VG_MIN_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_FIX_MAX_VG_MAX_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # max vg
            fbaout1 = ChLP.fba(model, exglcider; sense = max_sense)
            exglc = ChU.av(model, fbaout1, exglcider)
            ChU.bounds!(model, exglcider, exglc, exglc)

            # max cost
            fbaout = ChLP.fba(model, costider; sense = max_sense)

            LPDAT[FBA_Z_FIX_MAX_VG_MAX_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MAX_VG_MAX_COST, :fbaout, exp] = fbaout
        end

        # FBA_MAX_Z_MIN_COST
        let
            model = iJR.load_model("fva_models", exp)

            # max z
            fbaout1 = ChLP.fba(model, objider; sense = max_sense)
            objval = ChU.av(model, fbaout1, objider)
            ChU.bounds!(model, objider, objval, objval)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = min_sense)

            LPDAT[FBA_MAX_Z_MIN_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_MAX_Z_MIN_COST, :fbaout, exp] = fbaout
        end

        # FBA_MAX_Z_MAX_COST
        let
            model = iJR.load_model("fva_models", exp)

            # max z
            fbaout1 = ChLP.fba(model, objider; sense = max_sense)
            objval = ChU.av(model, fbaout1, objider)
            ChU.bounds!(model, objider, objval, objval)

            # max cost
            fbaout = ChLP.fba(model, costider; sense = max_sense)

            LPDAT[:FBA_MAX_Z_MAX_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[:FBA_MAX_Z_MAX_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_FIX_MIN_COST
        let
            model = iJR.load_model("fva_models", exp)
            exp_growth = Hd.val("D", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = min_sense)

            LPDAT[FBA_Z_FIX_MIN_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MIN_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_FIX_MAX_COST
        let
            model = iJR.load_model("fva_models", exp)
            exp_growth = Hd.val("D", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # max cost
            fbaout = ChLP.fba(model, costider; sense = max_sense)

            LPDAT[FBA_Z_FIX_MAX_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_FIX_MAX_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_VG_FIX_MIN_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # fix vg
            exp_exglc = -abs(Hd.uval("GLC", exp))
            ChU.bounds!(model, exglcider, exp_exglc, exp_exglc)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = min_sense)

            LPDAT[FBA_Z_VG_FIX_MIN_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_VG_FIX_MIN_COST, :fbaout, exp] = fbaout
        end

        # FBA_Z_VG_FIX_MAX_COST
        let
            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Hd.val("D", exp)
            ChU.bounds!(model, objider, exp_growth, exp_growth)

            # fix vg
            exp_exglc = -abs(Hd.uval("GLC", exp))
            ChU.bounds!(model, exglcider, exp_exglc, exp_exglc)

            # max cost
            fbaout = ChLP.fba(model, costider; sense = max_sense)

            LPDAT[FBA_Z_VG_FIX_MAX_COST, :model, exp] = ChU.compressed_model(model)
            LPDAT[FBA_Z_VG_FIX_MAX_COST, :fbaout, exp] = fbaout
        end

    end
end

## -------------------------------------------------------------------
LP_DAT_FILE = iJR.procdir("lp_dat_file.bson")
ChU.save_data(LP_DAT_FILE, LPDAT)

# ## -------------------------------------------------------------------
# using Plots
# let
#     METHODS = [
#         FBA_Z_FIX_MIN_COST, FBA_Z_FIX_MAX_VG_MIN_COST,
#         FBA_MAX_Z_MIN_COST, FBA_Z_FIX_MIN_VG_COST,
#         FBA_Z_VG_FIX_MIN_COST
#     ]
#     FLX_IDERS = ["GLC", "SUCC", "AC", "FORM"]
#     rxn_map = iJR.load_rxns_map()
#     color_pool = Plots.distinguishable_colors(length(FLX_IDERS))
#     ider_color = Dict(ider => color for (ider, color) in zip(FLX_IDERS, color_pool))

#     ps = Plots.Plot[]
#     for method in METHODS
#         p = plot(;title = string(method), xlabel = "exp", ylabel = "model")

#         for exp in EXPS
#             model = LPDAT[method, :model, exp]
#             fbaout = LPDAT[method, :fbaout, exp]

#             for Hd_ider in FLX_IDERS
#                 iJR_ider = rxn_map[Hd_ider]

#                 exp_val = Hd.uval(Hd_ider, exp) |> abs
#                 model_val = ChU.av(model, fbaout, iJR_ider) |> abs

#                 scatter!(p, [exp_val], [model_val]; label = "",
#                     color = ider_color[Hd_ider], m = 8
#                 )
#             end
#         end
#             push!(ps, p)
#     end

#     figname = string("fba_corrs")
#     dir = iJR.plotsdir()
#     UJL.mysavefig(ps, figname, dir)
# end


text:
"""
    install(force=true)

Downloads and installs a local etcd server (mostly for testing).
"""
function install(force=true)
    if !isempty(ETCD_DOWNLOAD_URI)
        if !ispath(ETCD_SRC_PATH)
            info(getlogger(current_module()), "Downloading etcd...")
            download(ETCD_DOWNLOAD_URI, ETCD_SRC_PATH)
        end

        if !ispath(ETCD_DEST_PATH) || force
            try
                run(`tar -xzvf $ETCD_SRC_PATH -C $ETCD_DEPS_PATH`)
            catch
                rm(src_path, recursive=true)
            end
        end
    end
end

"""
    start(timeout=-1) -> Future

Starts up the local etcd server (mostly for testing).
The `timeout` specifies a number of seconds to run the server for
(using `timeout` or `gtimeout`).
"""
function start(timeout=-1)
    info(getlogger(current_module()), "Starting etcd server...")

    if !ispath(ETCD_DEST_PATH)
        install()
    end

    if timeout > 0
        if is_apple()
            return spawn(`gtimeout $timeout $ETCD_BIN`)
        else
            return spawn(`timeout $timeout $ETCD_BIN`)
        end
    else
        return spawn(`$ETCD_BIN`)
    end
end


text:
using Random, LinearAlgebra, Plots; pyplot()
Random.seed!()

N = 10^5
data     = [[rand(),rand()] for _ in 1:N]
indata   = filter((x)-> (norm(x) <= 1), data)
outdata  = filter((x)-> (norm(x) > 1), data)
piApprox = 4*length(indata)/N
println("Pi Estimate: ", piApprox)

scatter(first.(indata),last.(indata), c=:blue, ms=1, msw=0)
scatter!(first.(outdata),last.(outdata), c=:red, ms=1, msw=0,
	xlims=(0,1), ylims=(0,1), legend=:none, ratio=:equal)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv")
df = CSV.read(data_path, DataFrame; delim=';')

# ## Model

using Turing

dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept)))
df.male = [g == "male" ? 1 : 0 for g in df.gender]
df.dept_id = [dept_map[de] for de in df.dept]

@model m13_2(applications, dept_id, male, admit) = begin
    sigma_dept ~ truncated(Cauchy(0, 2), 0, Inf)
    bm ~ Normal(0, 1)
    a ~ Normal(0, 10)
    a_dept ~ filldist(Normal(a, sigma_dept), 6)

    logit_p = a_dept[dept_id] + bm*male

    admit .~ BinomialLogit.(applications, logit_p)
end

chns = sample(
    m13_2(df.applications, df.dept_id, df.male, df.admit),
    NUTS(),
    1000
)

# ## Output

# \defaultoutput{}

# ## Original output

"""
Inference for Stan model: 359c2483e3bdbf74fd0484be27c2909b.
    3 chains, each with iter=4500; warmup=500; thin=1;
    post-warmup draws per chain=4000, total post-warmup draws=12000.

                   mean se_mean   sd     2.5%      25%      50%      75%    97.5%
    a_dept[1]      0.67    0.00 0.10     0.48     0.61     0.67     0.74     0.87
    a_dept[2]      0.63    0.00 0.12     0.40     0.55     0.63     0.71     0.85
    a_dept[3]     -0.58    0.00 0.08    -0.73    -0.63    -0.58    -0.53    -0.44
    a_dept[4]     -0.62    0.00 0.09    -0.79    -0.67    -0.62    -0.56    -0.45
    a_dept[5]     -1.06    0.00 0.10    -1.26    -1.13    -1.06    -0.99    -0.86
    a_dept[6]     -2.61    0.00 0.16    -2.92    -2.71    -2.60    -2.50    -2.31
    a             -0.59    0.01 0.66    -1.94    -0.97    -0.59    -0.21     0.73
    bm            -0.09    0.00 0.08    -0.25    -0.15    -0.09    -0.04     0.06
    sigma_dept     1.49    0.01 0.62     0.78     1.09     1.35     1.71     3.01
    lp__       -2602.27    0.04 2.27 -2607.74 -2603.48 -2601.90 -2600.61 -2598.95
               n_eff Rhat
    a_dept[1]   8191    1
    a_dept[2]   8282    1
    a_dept[3]  10974    1
    a_dept[4]  10034    1
    a_dept[5]  10783    1
    a_dept[6]  10344    1
    a           5839    1
    bm          7137    1
    sigma_dept  4626    1
    lp__        4200    1
""";


text:
##### NewtonsMethodAD

export NewtonsMethodAD

"""
    NewtonsMethodAD(f!::F!, x_init::A) where {F!, A <: AbstractArray}

A non-linear system of equations type.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct NewtonsMethodAD{FT, F!, A, JA} <: AbstractNonlinearSolverMethod{FT}
    "Function to find the root of"
    f!::F!
    "Initial guess"
    x_init::A
    "Storage"
    x1::A
    "Storage"
    J::JA
    "Storage"
    J⁻¹::JA
    "Storage"
    F::A
    function NewtonsMethodAD(f!::F!, x_init::A) where {F!, A <: AbstractArray}
        x1 = similar(x_init)
        J = similar(x_init, (length(x_init), length(x_init)))
        J⁻¹ = similar(J)
        F = similar(x_init)
        JA = typeof(J)
        FT = eltype(x_init)
        return new{FT, F!, A, JA}(f!, x_init, x1, J, J⁻¹, F)
    end
end

method_args(m::NewtonsMethodAD) = (m.x_init, m.x1, m.f!, m.F, m.J, m.J⁻¹)
function solve!(
    ::NewtonsMethodAD,
    x0::AT,
    x1::AT,
    f!::F!,
    F::FA,
    J::JA,
    J⁻¹::J⁻¹A,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {FA, J⁻¹A, JA, F! <: Function, AT, FT}

    x_history = init_history(soltype, AT)
    F_history = init_history(soltype, AT)
    if soltype isa VerboseSolution
        f!(F, x0)
        ForwardDiff.jacobian!(J, f!, F, x0)
        push_history!(x_history, x0, soltype)
        push_history!(F_history, F, soltype)
    end
    for i in 1:maxiters
        f!(F, x0)
        ForwardDiff.jacobian!(J, f!, F, x0)
        x1 .= x0 .- J \ F
        push_history!(x_history, x1, soltype)
        push_history!(F_history, F, soltype)
        if tol(x0, x1, F)
            return SolutionResults(
                soltype,
                x1,
                true,
                F,
                i,
                x_history,
                F_history,
            )
        end
        x0 = x1
    end
    return SolutionResults(
        soltype,
        x0,
        false,
        F,
        maxiters,
        x_history,
        F_history,
    )
end


text:
struct DAGSummary
    ctx::Union{SyncContext,Nothing}
    dag::DAG{DurationNS}
    work::DurationNS
    span::DurationNS
end

TaskDAGAnalyzers.work(s::DAGSummary) = s.work
TaskDAGAnalyzers.span(s::DAGSummary) = s.span

TaskDAGAnalyzers.parallelism(s) =
    Real(TaskDAGAnalyzers.work(s)) / Real(TaskDAGAnalyzers.span(s))

TaskDAGAnalyzers.summary(ctx::SyncContext) = DAGSummary(ctx, TaskDAGAnalyzers.dag(ctx))
TaskDAGAnalyzers.summary(dag::DAG) = DAGSummary(nothing, dag)

function DAGSummary(ctx::Union{SyncContext,Nothing}, dag::DAG)
    work = TaskDAGAnalyzers.work(dag)
    span = TaskDAGAnalyzers.span(dag)
    return DAGSummary(ctx, dag, work, span)
end

function Base.show(io::IO, ::MIME"text/plain", s::DAGSummary)
    print(io, "TaskDAGAnalyzers.summary:")
    print(io, '\n', "work: ", labeltext(s.work))
    printstyled(io, " (single-thread run-time T₁)", color = :light_black)
    print(io, '\n', "span: ", labeltext(s.span))
    printstyled(io, " (theoretical fastest run-time Tₒₒ)", color = :light_black)
    print(io, '\n', "parallelism (work/span): ", TaskDAGAnalyzers.parallelism(s))
end


text:
export BoxLinQuadTransform

abstract type Transform end

transform(t::Transform, x) = x

inverse(t::Transform, x) = x

struct NoTransform <: Transform
end

struct BoxLinQuadTransform{T}
    lb::Vector{T}
    ub::Vector{T}
    al::Vector{T}
    au::Vector{T}
    scaling::Bool
end

function BoxLinQuadTransform(lb, ub; scaling = false)
    al = @. min((ub - lb) / 2, (1 + abs(lb)) / 20)
    au = @. min((ub - lb) / 2, (1 + abs(ub)) / 20)
    BoxLinQuadTransform(lb, ub, al, au, scaling)
end

function transform(tr::BoxLinQuadTransform, x)
    map(tr.scaling ? tr.lb .+ (tr.ub .- tr.lb) .* x : x,
        tr.lb, tr.ub, tr.al, tr.au) do y, lb, ub, al, au
        yl = lb - 2 * al - (ub - lb) / 2
        yu = ub + 2 * au + (ub - lb) / 2
        z, r = y, 2 * (ub - lb + al + au)
        y = ifelse(y < yl, y + r * ceil((yl - y) / r),
            ifelse(y > yu, y - r * ceil((y - yu) / r), y))
        y = ifelse(y < lb - al, y + 2 * (lb - al - y),
            ifelse(y > ub + au, y - 2 * (y - ub - au), y))
        y = ifelse(y < lb + al, lb + (y - (lb - al))^2 / 4 / al,
            ifelse(y > ub - au, ub -(y - (ub + au))^2 / 4 / au, y))
        oftype(z, y)
    end
end

function inverse(tr::BoxLinQuadTransform, x)
    y = map(x, tr.lb, tr.ub, tr.al, tr.au) do y, lb, ub, al, au
        z, y = y, clamp(y, lb, ub)
        y = ifelse(y < lb + al, (lb - al) + 2 * sqrt(al * (y - lb)),
            ifelse(y > ub - au, (ub + au) - 2 * sqrt(au * (ub - y)), y))
        oftype(z, y)
    end
    tr.scaling ? (y .- tr.lb) ./ (tr.ub .- tr.lb) : y
end


text:
using Surrogates

#1D
x = [1.0,2.0,3.0]
y = [1.5,3.5,5.3]
lb = 0.0
ub = 7.0
my_linear_surr_1D = LinearSurrogate(x,y,lb,ub)
val = my_linear_surr_1D(5.0)
add_point!(my_linear_surr_1D,4.0,7.2)
add_point!(my_linear_surr_1D,[5.0,6.0],[8.3,9.7])

#ND
lb = [0.0,0.0]
ub = [10.0,10.0]
x = sample(5,lb,ub,SobolSample())
y = [4.0,5.0,6.0,7.0,8.0]
my_linear_ND = LinearSurrogate(x,y,lb,ub)
add_point!(my_linear_ND,(10.0,11.0),9.0)
add_point!(my_linear_ND,[(8.0,5.0),(9.0,9.5)],[4.0,5.0])
val = my_linear_ND((10.0,11.0))


text:
using JLD2
using Plots
using Printf
using Statistics

using Oceananigans
using Oceananigans.Advection
using Oceananigans.AbstractOperations
using Oceananigans.OutputWriters
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Forcings

include("Bickley.jl")

using .Bickley

function run(; Nh = 128,
               output_time_interval = 2,
               stop_time = 200,
               arch = CPU(),
               ν = 0,
               advection = WENO5())

    name = "oceananigans_unstable_bickley_jet_Nh$(Nh)_$(typeof(advection).name.wrapper)"

    grid = RegularCartesianGrid(size=(Nh, Nh, 1),
                                x = (-2π, 2π), y=(-2π, 2π), z=(0, 1),
                                topology = (Periodic, Periodic, Bounded))

    model = IncompressibleModel(architecture = arch,
                                 timestepper = :RungeKutta3,
                                   advection = advection,
                                        grid = grid,
                                     tracers = :c,
                                     closure = IsotropicDiffusivity(ν=ν, κ=ν),
                                    buoyancy = nothing)

    # ** Initial conditions **
    #
    # u, v: Large-scale jet + vortical perturbations
    #    c: Sinusoid

    # Parameters
    ϵ = 0.1 # perturbation magnitude
    ℓ = 0.5 # Gaussian width
    k = 0.5 # Sinusoidal wavenumber

    # Total initial conditions
    uᵢ(x, y, z) = Bickley.U(y) + ϵ * Bickley.ũ(x, y, ℓ, k)
    vᵢ(x, y, z) = ϵ * Bickley.ṽ(x, y, ℓ, k)
    cᵢ(x, y, z) = Bickley.C(y, grid.Ly)

    set!(model, u=uᵢ, v=vᵢ, c=cᵢ)

    progress(sim) = @info(@sprintf("Iter: %d, time: %.1f, Δt: %.3f, max|u|: %.2f",
                                   sim.model.clock.iteration, sim.model.clock.time,
                                   sim.Δt.Δt, maximum(abs, u.data.parent)))

    wizard = TimeStepWizard(cfl=1.0, Δt=1e-4, max_change=1.1, max_Δt=10.0)

    simulation = Simulation(model, Δt=wizard, stop_time=stop_time,
                            iteration_interval=10, progress=progress)

    # Output: primitive fields + computations
    u, v, w, c = primitives = merge(model.velocities, model.tracers)

    ω   = ComputedField( ∂x(v) - ∂y(u)      )
    ω²  = ComputedField( (∂x(v) - ∂y(u))^2  )
    c²  = ComputedField( c^2                )

    ∇c² = @at (Cell, Cell, Cell) ∂x(c)^2 + ∂y(c)^2
    ∇c² = ComputedField(∇c²)

    computations = (ω=ω, ω²=ω², c²=c², ∇c²=∇c²)
    outputs = merge(primitives, computations)

    save_grid = (file, model) -> file["serialized/grid"] = model.grid

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers, (ω=ω, ∇c²=∇c²)),
                                schedule = TimeInterval(output_time_interval),
                                init = save_grid,
                                prefix = name * "_fields",
                                force = true)

    averages = Dict(name => mean(ϕ, dims=(1, 2, 3))
                    for (name, ϕ) in zip(keys(computations), values(computations)))

    simulation.output_writers[:averages] =
        JLD2OutputWriter(model, averages,
                         schedule = TimeInterval(output_time_interval),
                         init = save_grid,
                         prefix = name * "_statistics",
                         force = true)

    @info "Running a simulation of an unstable Bickley jet with $(Nh)² degrees of freedom..."

    run!(simulation)

    return name
end

function analyze(name)
    @info "Analyzing the results of an unstable Bickley jet simulation..."

    statistics_file = jldopen(name * "_statistics.jld2")

    iterations = parse.(Int, keys(statistics_file["timeseries/t"]))
    grid = statistics_file["serialized/grid"]

    ∇c² = [statistics_file["timeseries/∇c²/$iter"][1, 1, 1] for iter in iterations]
    c²  = [statistics_file["timeseries/c²/$iter"][1, 1, 1] for iter in iterations]
    t   = [statistics_file["timeseries/t/$iter"][1, 1, 1] for iter in iterations]

    close(statistics_file)

    plot(t, 2π / grid.Ly .* ∇c² ./ c²)
    savefig(name * ".png")

    return nothing
end

function visualize(name, contours=false)
    @info "Making a fun movie about an unstable Bickley jet..."

    fields_file = jldopen(name * "_fields.jld2")

    iterations = parse.(Int, keys(fields_file["timeseries/t"]))
    grid = fields_file["serialized/grid"]

    xu, yu, zu = nodes((Face, Cell, Cell), grid)
    xω, yω, zω = nodes((Face, Face, Cell), grid)
    xc, yc, zc = nodes((Cell, Cell, Cell), grid)

    anim = @animate for (i, iteration) in enumerate(iterations)

        @info "    Plotting frame $i from iteration $iteration..."

        t = fields_file["timeseries/t/$iteration"]
        ω = fields_file["timeseries/ω/$iteration"][:, :, 1]
        u = fields_file["timeseries/u/$iteration"][:, :, 1]
        c = fields_file["timeseries/c/$iteration"][:, :, 1]

        kwargs = Dict(:xlabel => "x",
                      :ylabel => "y",
                      :aspectratio => 1,
                      :linewidth => 0,
                      :colorbar => true,
                      :clims => (-1, 1),
                      :xlims => (-grid.Lx/2, grid.Lx/2),
                      :ylims => (-grid.Ly/2, grid.Ly/2))

        contours && (kwargs[:levels] = range(-1, 1, length=31))
        plotter = contours ? contourf : heatmap

        u_plot = plotter(xu, yu, clamp.(u, -1, 1)'; color = :balance, kwargs...)
        ω_plot = plotter(xω, yω, clamp.(ω, -1, 1)'; color = :balance, kwargs...)
        c_plot = plotter(xc, yc, clamp.(c, -1, 1)'; color = :thermal, kwargs...)

        u_title = @sprintf("u at t = %.1f", t)
        ω_title = @sprintf("ω at t = %.1f", t)
        c_title = @sprintf("c at t = %.1f", t)

        plot(u_plot, ω_plot, c_plot,
             title = [u_title ω_title c_title],
             layout = (1, 3),
             size = (1600, 400))
    end

    gif(anim, name * ".gif", fps = 8)

    return nothing
end

for Nh in (512, 1024)
    #for advection in (CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(), UpwindBiasedFifthOrder(), WENO5())
    for advection in (WENO5(),)
        name = run(Nh=Nh, advection=advection)
        analyze(name)
        visualize(name)
    end
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits.QML.Optimization

testFun1(x) = 1.0*x[1]^2 + 100.0*x[2]^2
dtestFun1(x) = [2.0*x[1], 200.0*x[2]]

const ϵ = 0.000001

val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=1.0, maxItr=10000000)
#val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=1.0, maxItr=10000000, debug=true)
# Standard Unable
# Momentum Unable
# RMSprop      22
# Adam        384
# Moj           4
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ
#@test itr == 3
@test itr < 20

val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=0.1, maxItr=10000000)
# Standard Unable
# Momentum    404
# RMSprop      94
# Adam        532
# Moj           4
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ
#@test itr == 3
@test itr < 20

val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=0.01, maxItr=10000000)
# Standard Unable
# Momentum    720
# RMSprop    2896
# Adam       4580
# Moj           4
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ
@test itr == 3

val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=0.001, maxItr=10000000)
# Standard  9549
# Momentum  9383
# RMSprop  15849
# Adam     17101
# Moj          3
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ
@test itr == 3

val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], α=0.0001, maxItr=10000000)
# Standard  95561
# Momentum  95398
# RMSprop  108080
# Adam     109290
# Moj           3
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ
@test itr == 3

# Test other methods
using QuantumCircuits.QML.Optimization: StandardGD, Momentum, RMSprop, Adam

# Standard gradient descent
opt = StandardGD(0.001)
val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], maxItr=10000000, opt=opt)
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ

# Momentum optimization algoritm
opt = Momentum{Float64}(0.1, 0.9, 2)
val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], maxItr=10000000, opt=opt)
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ

# RMSprop optimization algoritm
opt = RMSprop{Float64}(0.1, 0.9, 2)
val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], maxItr=10000000, opt=opt)
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ

# Adam optimization algoritm
opt = Adam{Float64}(0.1, 0.9, 0.999, 2)
val, x, itr = gradientDescent(testFun1, dtestFun1, [10.0, 10.0], maxItr=10000000, opt=opt)
@test abs(val) < ϵ
@test abs(transpose(x)*x) < ϵ


text:
include("misfit.jl")
include("modeling.jl")


text:
resODE = Dict{Symbol,Dict}()
const S = Solvers

for (n,tc) in totest
    res = Dict{Solver,Any}()
    for (solverfn,solver) in S.ODEsolvers
        if isapplicable(solver, tc) && isadaptive(solver)
            suite = TestSuite(tc, solver, abstols, reltols, [NaN])
            res[solver] = run_ode_testsuite(suite)
        end
    end
    resODE[n] = res
end


text:
module GUI

export HtmlElement, hyperhtml, Context, Content,
       render

using DataStructures, LinkedLists, StructArrays, NamedTupleTools, AbstractTrees
using CanonicalTraits
using Gumbo, Cascadia, EzXML
using HAML, Mustache, Hyperscript
using ReusePatterns

using HTTP, HTTP.URIs

include("type/abstract.jl")            #

include("support/macros.jl")           # @sym
include("type/maybe.jl")               # MaybeString, MaybeFloat, isstring, isfloat ...

include("support/misc.jl")             # first(nothing), isless(x, nothing) ...
include("support/strings.jl")          # isuppercase(String), uppercasefirst

include("support/linkedlists.jl")      # getindex(list, idx::Int)
include("support/dict.jl")             # getkeys::Vector, getvalues::Vector
include("type/prototype.jl")           # specific NamedTuple prototypes

include("html/hyperscript_html.jl")    # HTML5 tags
#include("css/hyperscript_css.jl")     # CSS3 tags
include("type/algebra.jl")             # Content, Context, HtmlElement

include("support/ezxml.jl")            # prettyxml
include("support/gumbo.jl")            # typeofnode(x::HTMLNode)
include("support/cascadia.jl")         #
include("support/hyperscript.jl")      # convert to Gumbo, push!, pushlast!

include("html/template.jl")            #
include("html/elements.jl")            #
include("html/attributes/boolean.jl")  #
include("html/googlefonts.jl")         # googlefontlink

end # module GUI


text:
using DimensionalData, Test

using DimensionalData: Name, NoName

@testset "Name" begin
    @test Name(:x) === Name{:x}()
    @test Symbol(Name(:x)) === :x
    @test string(Name(:x)) === "x"
    @test Name(Name(:x)) === Name(:x)
    @test convert(String, Name(:x)) === "x"
end

@testset "NoName" begin
    @test Symbol(NoName()) === Symbol("")
    @test string(NoName()) === ""
    @test Name(NoName()) === NoName()
    @test convert(String, NoName()) === ""
end


text:
using Test
using SurvivalSignature

include("percolation.jl")
include("signature.jl")


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x = y = -15:0.33:15
fz(x,y) = sin.(sqrt.(x.^2 + y.^2))./sqrt.(x.^2+y.^2)
fxy = [fz(x,y) for x in x, y in y]
@gsp " " :-  # reset session
@gsp(:-,x, y, fxy, "w pm3d", "set ylabel 'y' textcolor rgb 'red'",
    "set xlabel 'x' textcolor rgb 'blue'")
save(term="pngcairo size 800,800", output="plt3d_ex3_1.png")


text:
@system Count begin
    pheno: phenology ~ hold
    NU: nodal_units ~ hold

    leaves_initiated(pheno.leaves_initiated) ~ track::int
    leaves_appeared(pheno.leaves_appeared) ~ track::int
    # leaves_detached(pheno.leaves_detached) ~ track::int

    leaves_growing(x=NU["*"].leaf.growing) => sum(x) ~ track::int
    leaves_mature(x=NU["*"].leaf.mature) => sum(x) ~ track::int
    leaves_dropped(x=NU["*"].leaf.dropped) => sum(x) ~ track::int

    detach_to_maintain_max_leaf_number => 19 ~ preserve::int(parameter)
    detach_to_maintain_min_leaf_number => 15 ~ preserve::int(parameter)

    leaves_detached(leaves_mature, leaves_dropped, detach_to_maintain_max_leaf_number, detach_to_maintain_min_leaf_number) => begin
        if leaves_mature - leaves_dropped <= detach_to_maintain_max_leaf_number
            0
        else
            leaves_mature - leaves_dropped - detach_to_maintain_min_leaf_number
        end
    end ~ track::int

    leaves_real(a=leaves_mature, d=leaves_dropped, c=leaves_detached) => (a - d -c ) ~ track::int

end


text:
export LB1, LB2, LB3, LB4, LB5, LB6

"""
    compute_bound(LB, fpath)

Compute bound `LB` for the instance-file fpath
"""
function compute_bound(LB, fpath)
    i = read(fpath, Instance)
    return LB(i)
end

"""
    compute_bounds(boundset, fpath)

Computes the bounds in `boundset` for the instance-file `fpath` and returns a dictionary with `boundname => value` pairs.
"""
function compute_bounds(boundset, fpath)
    i = read(fpath, Instance)
    d = Dict{String, Int}()
    for lb in boundset
        v = lb(i)
        d[string(lb)] = v
    end
    return d
end

"""
    LB1(jobs)

Calculate a lower bound on the truck waiting time:

The sum of the processing times of truck-related jobs
"""
LB1(jobs) = sum(t_processing, Iterators.filter(istruck, jobs))

LB1(i::Instance) = LB1(jobs(i))

"""
    LB2(jobs)

Calculate a lower bound on the makespan:

The maximum of the minimum completion times of all jobs.
"""
LB2(jobs) = maximum(j -> t_arrival(j) + t_processing(j), jobs)

LB2(i::Instance) = LB2(jobs(i))

"""
    LB3(jobs, ncranes)

Calculate a lower bound on the makespan:

Balancing total processing time over cranes
"""
LB3(jobs, ncranes) = Int(ceil(sum(t_processing, jobs) / ncranes ))

LB3(i::Instance) = LB3(jobs(i), ncranes(i))

"""
    LB4(jobs, l0)

Calculate a lower bound on the makespan:

# TODO: Implement bound
"""

LB4(jobs, l0) = LB3 + minimaltravel(jobs, l0)

LB4(i::Instance) = LB4(jobs(i), crane_starting_pos(i))

"""
    LB5(jobs)

Calculate a lower bound on the truck waiting time:

Compute the total truck waiting time _after_ taking into account job precedences. Analogous to LB1, but with precedences.
"""
function LB5(jobs)
    Ecompl = compute_earliest_completion_times(jobs)
    V = Iterators.filter(istruck, jobs)
    mapreduce(j -> Ecompl[id(j)] - t_arrival(j), +,  V)
end

LB5(i::Instance) = LB5(jobs(i))

"""
    LB6(jobs)

Calculate a lower bound on the makespan:

Compute the maximum earliest completion time of the jobs.
"""
function LB6(jobs)
    Ecompl = compute_earliest_completion_times(jobs)
    return maximum(Ecompl)
end

LB6(i::Instance) = LB6(jobs(i))

"""
    interlocdist(jobs)

Returns the sorted distances between unique job-locations
"""
function interlocdist(jobs)
    locsorted = sort!(unique(loc.(jobs)))
    return sort!(map(l -> -(l...), Iterators.zip(locsorted[1:end-1], locsorted[2:end])))
end

"""
    minimaltravel(jobs, ncranes)

    Calculate the minimal travel distance for ncranes to visit all positions
"""
function minimaltravel(jobs, l0)
    Djobs = interlocdist(jobs)

end

function compute_earliest_completion_times(jobs)
    A = t_arrival.(jobs)
    P = t_processing.(jobs)
    C = A .+ P

    pg = precedence_graph(jobs)
    ids = topological_sort_by_dfs(pg)
    for id in ids
        I = map(inn -> C[inn] ,inneighbors(pg, id))
        if !isempty(I)
            C[id] = max(maximum(I), A[id]) + P[id]
        end
    end
    return C
end

makespan_bound_functions = [LB2, LB3, LB4, LB6]
t_truckwaiting_bound_functions = [LB1, LB5]