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

"""
Unit test runner for the TestTools.jl package.
"""

# --- Imports

# Standard library
using Test

# External packages
using Suppressor

# Local package
using TestTools.jltest

# --- Preparations

# Change to test directory
#
# Note: this is needed for consistency of results when tests are run via
# `jltest runtests.jl` and `import Pkg; Pkg.test()`
cd(@__DIR__)

# --- Normal unit tests

tests = [
    joinpath(@__DIR__, "pkg_tests.jl"),
    joinpath(@__DIR__, "jltest", "EnhancedTestSet_passing_tests.jl"),
    joinpath(@__DIR__, "jltest", "EnhancedTestSet_fail_fast_tests.jl"),
    joinpath(@__DIR__, "jlcodestyle", "cli_tests.jl"),
    joinpath(@__DIR__, "jlcoverage", "cli_tests.jl"),
    joinpath(@__DIR__, "jlcoverage", "utils_tests.jl"),
]
jltest.run_tests(tests; desc="jltest")

# --- jltest unit tests that have expected failures and errors

local log_message
local error_type, error_message

# EnhancedTestSet with failing tests
println()
test_file = joinpath(@__DIR__, "jltest", "EnhancedTestSet_failing_tests.jl")
output = strip(
    @capture_out begin
        try
            @testset EnhancedTestSet "EnhancedTestSet" begin
                global log_message = strip(
                    @capture_err begin
                        jltest.run_tests(test_file; desc="failing tests")
                    end
                )
            end
        catch error
            bt = catch_backtrace()
            global error_type = typeof(error)
            global error_message = sprint(showerror, error, bt)
        end
    end
)

print("jltest/EnhancedTestSet_failing_tests: ")
@testset EnhancedTestSet "EnhancedTestSet: check for expected test failures" begin
    @test log_message ==
        "[ Info: For EnhancedTestSet_failing_tests.jl, 6 failures and 1 error are expected."

    @test error_type == TestSetException
    @test error_message ==
        "Some tests did not pass: 7 passed, 6 failed, 1 errored, 0 broken."

    # Check output from EnhancedTestSet
    expected_output = strip(
        """
        $(joinpath("jltest", "EnhancedTestSet_failing_tests")): .......

        Test Summary:                                | Pass  Fail  Error  Total
        EnhancedTestSet                              |    7     6      1     14
          failing tests                              |    7     6      1     14
            EnhancedTestSet: Array equality test     |          1             1
            EnhancedTestSet: Dict equality test      |          1             1
            EnhancedTestSet: String equality test    |          1             1
            EnhancedTestSet: Boolean expression test |          1             1
            EnhancedTestSet: Exception test          |                 1      1
            EnhancedTestSet: inequality test         |          1             1
            EnhancedTestSet: Matrix equality test    |          1             1
        """
    )
    @test output == expected_output
end

# EnhancedTestSet with nested test sets
println()
test_file = joinpath(@__DIR__, "jltest", "EnhancedTestSet_nested_test_set_tests.jl")
output = strip(
    @capture_out begin
        try
            @testset EnhancedTestSet "EnhancedTestSet" begin
                global log_message = strip(
                    @capture_err begin
                        jltest.run_tests(test_file; desc="nested test set tests")
                    end
                )
            end
        catch error
            bt = catch_backtrace()
            global error_type = typeof(error)
            global error_message = sprint(showerror, error, bt)
        end
    end
)

print("jltest/EnhancedTestSet_nested_test_set_tests: ")
@testset EnhancedTestSet "EnhancedTestSet: check for expected test failures" begin
    @test log_message ==
        "[ Info: For EnhancedTestSet_nested_test_set_tests.jl, 2 failures and 0 error are expected."

    @test error_type == TestSetException
    @test error_message ==
        "Some tests did not pass: 2 passed, 2 failed, 0 errored, 0 broken."

    # Check output from EnhancedTestSet
    expected_output = strip(
        """
        $(joinpath("jltest", "EnhancedTestSet_nested_test_set_tests")): ..

        Test Summary:                                         | Pass  Fail  Total
        EnhancedTestSet                                       |    2     2      4
          nested test set tests                               |    2     2      4
            EnhancedTestSet: nested inherited EnhancedTestSet |          1      1
              Nested Inherited Test Set                       |          1      1
            EnhancedTestSet: nested DefaultTestSet            |          1      1
              DefaultTestSet Nested in EnhancedTestSet        |          1      1
        """
    )
    @test output == expected_output
end

# utils.jl
println()
test_file = joinpath(@__DIR__, "jltest", "utils_tests.jl")
output = strip(
    @capture_out begin
        try
            @testset EnhancedTestSet "jltest" begin
                global log_message = strip(
                    @capture_err begin
                        jltest.run_tests(test_file; desc="utils tests")
                    end
                )
            end
        catch error
            bt = catch_backtrace()
            global error_type = typeof(error)
            global error_message = sprint(showerror, error, bt)
        end
    end
)

print("jltest/utils_tests: ")
@testset EnhancedTestSet "jltest.utils: check for expected test failures" begin
    @test log_message ==
        "[ Info: For utils_tests.jl, 13 failures and 0 errors are expected."

    @test error_type == TestSetException
    @test error_message ==
        "Some tests did not pass: 105 passed, 13 failed, 0 errored, 0 broken."

    # Check output from EnhancedTestSet
    expected_output = strip(
        """
        $(joinpath("jltest", "utils_tests")): .......................................................

        Test Summary:                                    | Pass  Fail  Total
        jltest                                           |  105    13    118
          utils tests                                    |  105    13    118
            jltest.run_tests(): basic tests              |   70    11     81
              test set                                   |    2            2
              test set                                   |    2            2
              test set                                   |    8     2     10
                failing tests                            |    1     1      2
                some tests                               |    2            2
                more tests                               |    2            2
              test set                                   |   10     2     12
                failing tests                            |    1     1      2
                some tests                               |    2            2
                more tests                               |    2            2
                some tests                               |    2            2
              test set                                   |    2            2
              test description                           |    1     1      2
              test set                                   |    1     1      2
                failing tests                            |    1     1      2
              failing tests                              |    1     1      2
              test set                                   |    6     2      8
                failing tests                            |    1     1      2
                some tests                               |    2            2
              test set                                   |    8     2     10
                failing tests                            |    1     1      2
                some tests                               |    2            2
                more tests                               |    2            2
            jltest.run_tests(): log message tests        |    9            9
            jltest.run_tests(): current directory checks |    4            4
            jltest.run_tests(): JLTEST_FAIL_FAST tests   |   16     2     18
              test set                                   |    3     1      4
                failing tests                            |    1     1      2
                some tests                               |    2            2
              test set                                   |    3     1      4
                failing tests                            |    1     1      2
                some tests                               |    2            2
            jltest.find_tests()                          |    4            4
            jltest: Pkg.test() tests                     |    2            2
            """,
    )

    @test output == expected_output
end

# cli.jl
println()
test_file = joinpath(@__DIR__, "jltest", "cli_tests.jl")
output = strip(
    @capture_out begin
        try
            @testset EnhancedTestSet "jltest" begin
                global log_message = strip(
                    @capture_err begin
                        jltest.run_tests(test_file; desc="cli tests")
                    end
                )
            end
        catch error
            bt = catch_backtrace()
            global error_type = typeof(error)
            global error_message = sprint(showerror, error, bt)
        end
    end
)

print("jltest/cli_tests: ")
@testset EnhancedTestSet "jltest.cli: check for expected test failures" begin
    @test log_message == "[ Info: For cli_tests.jl, 5 failures and 0 errors are expected."

    @test error_type == TestSetException
    @test error_message ==
        "Some tests did not pass: 59 passed, 5 failed, 0 errored, 0 broken."

    # Check output from EnhancedTestSet
    expected_output = strip(
        """
        $(joinpath("jltest", "cli_tests")): ......................................

        Test Summary:                     | Pass  Fail  Total
        jltest                            |   59     5     64
          cli tests                       |   59     5     64
            jltest.cli.parse_args()       |   12           12
            jltest.cli.run(): basic tests |   46     5     51
              All tests                   |    4            4
              failing tests               |    1     1      2
              All tests                   |    6     2      8
                failing tests             |    1     1      2
                some tests                |    2            2
              All tests                   |    8     2     10
                failing tests             |    1     1      2
                some tests                |    2            2
                more tests                |    2            2
            jltest.cli.run(): error cases |    1            1
        """
    )
    @test output == expected_output
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors
using DelimitedFiles

path = joinpath("test_chaboche", "unitelement_results.rpt")
data = readdlm(path, Float64; skipstart=4)
ts = data[:,1]
s11_ = data[:,2]
s12_ = data[:,3]
s13_ = data[:,4]
s22_ = data[:,5]
s23_ = data[:,6]
s33_ = data[:,7]
e11_ = data[:,8]
e12_ = data[:,9]
e13_ = data[:,10]
e22_ = data[:,11]
e23_ = data[:,12]
e33_ = data[:,13]

strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)]

parameters = ChabocheParameterState(E = 200.0e3,
                                    nu = 0.3,
                                    R0 = 100.0,
                                    Kn = 100.0,
                                    nn = 10.0,
                                    C1 = 10000.0,
                                    D1 = 100.0,
                                    C2 = 50000.0,
                                    D2 = 1000.0,
                                    Q = 50.0,
                                    b = 0.1)
chabmat = Chaboche(parameters = parameters)
s33s = [chabmat.variables.stress[3,3]]
for i=2:length(ts)
    dtime = ts[i]-ts[i-1]
    dstrain = fromvoigt(SymmetricTensor{2,3,Float64}, strains[i]-strains[i-1]; offdiagscale=2.0)
    chabmat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain)
    integrate_material!(chabmat)
    update_material!(chabmat)
    push!(s33s, chabmat.variables.stress[3,3])
end
@test isapprox(s33s, s33_; rtol=0.05)

# function test_chab()
#     mat = Material(Chaboche, ())
#     mat.properties.youngs_modulus = 200.0e3
#     mat.properties.poissons_ratio = 0.3
#     mat.properties.yield_stress = 100.0
#     mat.properties.K_n = 100.0
#     mat.properties.n_n = 10.0
#     mat.properties.C_1 = 10000.0
#     mat.properties.D_1 = 100.0
#     mat.properties.C_2 = 50000.0
#     mat.properties.D_2 = 1000.0
#     mat.properties.Q = 50.0
#     mat.properties.b = 0.1
#
#     mat.stress = zeros(6)
#     s33s = [0.0]
#     for i=2:length(ts)
#         dtime = ts[i]-ts[i-1]
#         dstrain = strains[i]-strains[i-1]
#         mat.dtime = dtime
#         mat.dstrain = dstrain
#         integrate_material!(mat)
#         mat.time += mat.dtime
#         mat.strain .+= mat.dstrain
#         mat.stress .+= mat.dstress
#         mat.properties.plastic_strain .+= mat.properties.dplastic_strain
#         mat.properties.backstress1 .+= mat.properties.dbackstress1
#         mat.properties.backstress2 .+= mat.properties.dbackstress2
#         mat.properties.R += mat.properties.dR
#         mat.properties.cumulative_equivalent_plastic_strain += mat.properties.dcumulative_equivalent_plastic_strain
#         push!(s33s, copy(mat.stress[3]))
#     end
#     # @test isapprox(s33s, s33_; rtol=0.001)
# end

# using BenchmarkTools
# @btime test_chab()

ext: jl


text:
using NEOS.MPSWriter, Base.Test

@testset "getrowsense" begin
    # LE, GE, Eq, Ranged
    row_sense, hasranged = MPSWriter.getrowsense([-Inf, 0.], [0., Inf])
    @test row_sense == [:(<=), :(>=)]
    @test hasranged == false

    row_sense, hasranged = MPSWriter.getrowsense([1., -1.], [1., 1.])
    @test row_sense == [:(==), :ranged]
    @test hasranged == true

    @test_throws Exception MPSWriter.getrowsense([1.], [1., 1.])
    @test_throws Exception MPSWriter.getrowsense([-Inf], [Inf])
end

@testset "writecolumns!" begin
    io = IOBuffer()

    colnames  = ["V$i" for i in 1:2]
    rownames  = ["C$i" for i in 1:2]

    MPSWriter.writecolumns!(io, [0. 0.]', [:Cont], [1.], :Min, colnames, rownames)
    @test String(take!(io)) == "COLUMNS\n    V1        OBJ       1\n"

    MPSWriter.writecolumns!(io, [0. 0.]', [:Cont], [1.], :Max, colnames, rownames)
    @test String(take!(io)) == "COLUMNS\n    V1        OBJ       -1\n"

    MPSWriter.writecolumns!(io, [0. 0.], [:Bin, :Cont], [3., 4.], :Max, colnames, rownames)
    @test String(take!(io)) == "COLUMNS\n    MARKER    'MARKER'                 'INTORG'\n    V1        OBJ       -3\n    MARKER    'MARKER'                 'INTEND'\n    V2        OBJ       -4\n"

    MPSWriter.writecolumns!(io, [0. 0.], [:Fixed, :Int], [3., 4.], :Max, colnames, rownames)
    @test String(take!(io)) == "COLUMNS\n    V1        OBJ       -3\n    MARKER    'MARKER'                 'INTORG'\n    V2        OBJ       -4\n    MARKER    'MARKER'                 'INTEND'\n"

    @test_throws Exception MPSWriter.writecolumns!(io, [0. 0.]', [:Cont], [1.], :badsense, colnames, rownames)
    @test_throws Exception MPSWriter.writecolumns!(io, [0. 0.]', [:badtype], [1.], :Min, colnames, rownames)
    @test_throws Exception MPSWriter.writecolumns!(io, [0. 0.]', [:Cont, :Cont], [1.], :Min, colnames, rownames)

    close(io)
end

@testset "writecolumn!" begin

    colnames  = ["V$i" for i in 1:2]
    rownames  = ["C$i" for i in 1:2]

    for Aty in [SparseMatrixCSC{Float64, Int}, Array{Float64, 2}]
        io = IOBuffer()
        A = convert(Aty, [1. 0.; 1.5 1.4])
        MPSWriter.writecolumn!(io, A, 1, colnames, rownames)
        @test String(take!(io)) == "    V1        C1        1\n    V1        C2        1.5\n"
        MPSWriter.writecolumn!(io, A, 2, colnames, rownames)
        @test String(take!(io)) == "    V2        C2        1.4\n"
        close(io)
    end
end

@testset "writerows!" begin
    rownames  = ["C$i" for i in 1:4]
    io = IOBuffer()
    MPSWriter.writerows!(io, [:(<=), :(>=), :(==), :ranged], rownames)
    s = String(take!(io))
    @test s == "ROWS\n N  OBJ\n L  C1\n G  C2\n E  C3\n E  C4\n"
    @test_throws Exception MPSWriter.writerows!(io, [:badsym], rownames)
    close(io)
end

@testset "writerhs!" begin
    io = IOBuffer()
    rownames  = ["C$i" for i in 1:4]
    MPSWriter.writerhs!(io, [-Inf, -1., 0., 0.], [0., 1., 0., Inf], [:(<=), :ranged, :(==), :(>=)], rownames)
    s = String(take!(io))
    @test s == "RHS\n    rhs       C1        0\n    rhs       C2        -1\n    rhs       C3        0\n    rhs       C4        0\n"
    @test_throws Exception MPSWriter.writeMPS!(io, [-Inf], [0., 1], [:(<=), :(==)], rownames)
    close(io)
end

@testset "writeranges!" begin
    io = IOBuffer()
    rownames  = ["C$i" for i in 1:4]
    MPSWriter.writeranges!(io, [-Inf, -2.5, 0., 1.], [1., 1., 1., 2.], [:(<=), :ranged, :ranged, :ranged], rownames)
    s = String(take!(io))
    @test s == "RANGES\n    rhs       C2        3.5\n    rhs       C3        1\n    rhs       C4        1\n"
    close(io)
end

@testset "writebounds!" begin
    io = IOBuffer()
    colnames  = ["V$i" for i in 1:4]
    # (-Inf, Inf) = FR
    # (  * , Inf) = PL
    # (  * ,  x ) = UP
    # (-Inf,  * ) = MI
    # (  x ,  * ) = LO
    MPSWriter.writebounds!(io, [-Inf, 1., 0., -Inf], [Inf, Inf, Inf, 0.], colnames)

    s = String(take!(io))

    @test s == "BOUNDS\n FR BOUNDS    V1\n PL BOUNDS    V2\n LO BOUNDS    V2        1\n PL BOUNDS    V3\n UP BOUNDS    V4        0\n MI BOUNDS    V4\n"
    @test_throws Exception MPSWriter.writebounds!(io, [-Inf], [Inf, 0.], colnames)
    close(io)
end

@testset "writebound!" begin
    io = IOBuffer()
    colnames  = ["V$i" for i in 1:4]
    MPSWriter. writebound!(io, "FR", 1, colnames[1])
    @test String(take!(io)) == " FR BOUNDS    V1\n"
    MPSWriter. writebound!(io, "MI", 1, colnames[1])
    @test String(take!(io)) == " MI BOUNDS    V1\n"
    MPSWriter. writebound!(io, "PL", 1, colnames[1])
    @test String(take!(io)) == " PL BOUNDS    V1\n"

    @test_throws Exception MPSWriter. writebound!(io, "LO", 1, "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "UP", 1, "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "badstring", 1, "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "FR", 1, 1., "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "MI", 1, 1., "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "PL", 1, 1., "BadVar")
    @test_throws Exception MPSWriter. writebound!(io, "badstring", 1, "BadVar")

    MPSWriter. writebound!(io, "LO", 1, 1.0, colnames[1])
    @test String(take!(io)) == " LO BOUNDS    V1        1\n"
    MPSWriter. writebound!(io, "UP", 1, 1.0, colnames[1])
    @test String(take!(io)) == " UP BOUNDS    V1        1\n"
    MPSWriter. writebound!(io, "LO", 100, 1.0, "V100")
    @test String(take!(io)) == " LO BOUNDS    V100      1\n"
    MPSWriter. writebound!(io, "UP", 100, 1.0, "V100")
    @test String(take!(io)) == " UP BOUNDS    V100      1\n"
    close(io)
end

@testset "SOS" begin
    colnames  = ["V$i" for i in 1:4]
    io = IOBuffer()
    MPSWriter.writesos!(io, [(1, [1,2,3], [1.,3.,2.]), (2, [1,2,3], [2.,1.,3.])], 3, colnames)
    s = String(take!(io))
    @test s == "SOS\n S1 SOS1\n    V1        1\n    V2        3\n    V3        2\n S2 SOS2\n    V1        2\n    V2        1\n    V3        3\n"

    @test_throws Exception MPSWriter.writesos!(io, [(1, [1,2,3], [1.,2.,3.])], 2, colnames)
    @test_throws Exception MPSWriter.writesos!(io, [(1, [0,2,3], [1.,2.,3.])], 3, colnames)

    close(io)
end

@testset "Quad" begin
    io = IOBuffer()
    colnames  = ["V$i" for i in 1:4]
    MPSWriter.writequad!(io, sparse([1, 1, 2, 1], [1, 2, 2, 1], [0.5, 1.5, 1, 0.5]), :Min, colnames)
    @test String(take!(io)) == "QMATRIX\n    V1       V1        1\n    V1       V2        1.5\n    V2       V2        1\n"
    MPSWriter.writequad!(io, [1.75 0;1.4 1], :Max, colnames)
    @test String(take!(io)) == "QMATRIX\n    V1       V1        -1.75\n    V2       V1        -1.4\n    V2       V2        -1\n"
    close(io)
end

@testset "write" begin
const MPSFILE = """NAME          TestModel
ROWS
 N  OBJ
 G  CON1
 G  CON2
 E  CON3
 E  CON4
COLUMNS
    V1        CON1      1
    V2        CON2      1
    V3        CON3      1
    MARKER    'MARKER'                 'INTORG'
    V4        OBJ       1
    MARKER    'MARKER'                 'INTEND'
    V5        CON4      1
    V5        OBJ       -1
    V6        CON4      1
    V7        CON4      1
    MARKER    'MARKER'                 'INTORG'
    V8        OBJ       0
    MARKER    'MARKER'                 'INTEND'
RHS
    rhs       CON1      0
    rhs       CON2      2
    rhs       CON3      1
    rhs       CON4      1
RANGES
    rhs       CON3      1.5
BOUNDS
 UP BOUNDS    V1        3
 MI BOUNDS    V1
 UP BOUNDS    V2        3
 MI BOUNDS    V2
 UP BOUNDS    V3        3
 MI BOUNDS    V3
 PL BOUNDS    V4
 LO BOUNDS    V4        5.5
 UP BOUNDS    V5        1
 UP BOUNDS    V6        1
 UP BOUNDS    V7        1
 UP BOUNDS    V8        1
SOS
 S2 SOS1
    V5        1
    V6        2
    V7        3
QMATRIX
    V1       V1        2
    V2       V1        -1.1
    V1       V2        -1.1
    V2       V2        2
ENDATA
"""
    io = IOBuffer()
    MPSWriter.write(io,
    [
    1 0 0 0 0 0 0 0;
    0 1 0 0 0 0 0 0;
    0 0 1 0 0 0 0 0;
    0 0 0 0 1 1 1 0
    ],
    [-Inf, -Inf, -Inf, 5.5, 0, 0, 0, 0],
    [3, 3, 3, Inf, 1, 1, 1, 1],
    [0,0,0,-1,1,0,0,0],
    [0, 2, 1, 1],
    [Inf, Inf, 2.5, 1],
    :Max,
    [:Cont, :Cont, :Cont, :Int, :Cont, :Cont, :Cont, :Bin],
    MPSWriter.SOS[(2, [5,6,7], [1,2,3])],
    sparse([1,1,2,2], [1,2,1,2], [-2,1.1,1.1,-2], 8,8),
    "TestModel",
    ["V$i" for i in 1:8],
    ["CON$i" for i in 1:4]
    )
    @test String(take!(io)) == MPSFILE
    close(io)
end


text:
export Agent

"""
    Agent(;kwargs...)

One of the most commonly used [`AbstractAgent`](@ref).

Generally speaking, it does nothing but

1. Pass observation to the policy to generate an action
1. Update the buffer using the `observation => action` pair
1. Update the policy with the newly updated buffer

# Keywords & Fields

- `π`::[`AbstractPolicy`](@ref): the policy to use
- `buffer`::[`AbstractTurnBuffer`](@ref): used to store transitions between agent and environment
- `role=:DEFAULT`: used to distinguish different agents
"""
Base.@kwdef mutable struct Agent{P<:AbstractPolicy, B<:AbstractTurnBuffer, R} <: AbstractAgent
    π::P
    buffer::B
    role::R = :DEFAULT
end

(agent::Agent)(obs::Observation) = agent.π(obs)

function update!(agent::Agent, experience::Pair)
    update!(agent.buffer, experience, agent)
    update!(agent.π, agent.buffer)
end

function update!(buffer::AbstractTurnBuffer, experience::Pair, agent::Agent)
    push!(agent.buffer, experience)
end

function update!(
    buffer::AbstractTurnBuffer,
    experience::Pair,
    agent::Agent{<:QBasedPolicy{<:Union{PrioritizedDQNLearner,RainbowLearner}}},
)
    push!(priority(agent.buffer), agent.π.learner.default_priority)
    push!(agent.buffer, experience)
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

"""
Optimization interface consists of two functions `optimizationvariables` and `updateoptimizationvariables`.
`optimizationvariables` packs variables to be optimized into a vector.
`updateoptimizationvariables` receives a vector of variables and creates a new optical system with the variable values.
"""
module Optimization

using ..OpticSim
using ..OpticSim: detectorsize, temperature, pressure

"""
    optimizationvariables(a::AxisymmetricOpticalSystem{T}) -> Vector{T}

Pack variables that have been marked to be optimized into a vector in a form suitable for the optimizer.
Variables are marked for optimization by having a true value in the `:OptimizeName` column, where `Name` can be `Radius`, `Thickness` or `Conic`.
    """
function optimizationvariables(a::AxisymmetricOpticalSystem{T}) where {T<:Real}
    prescription = a.prescription

    radiusvariables = in(:OptimizeRadius, propertynames(prescription)) ? prescription[prescription[!, :OptimizeRadius], :Radius] : []
    radiuslower = [typemin(T) for _ in 1:length(radiusvariables)]
    radiusupper = [typemax(T) for _ in 1:length(radiusvariables)]
    thicknessvariables = in(:OptimizeThickness, propertynames(prescription)) ? prescription[prescription[!, :OptimizeThickness], :Thickness] : []
    thicknesslower = [zero(T) for _ in 1:length(thicknessvariables)] # TODO
    thicknessupper = [T(100) for _ in 1:length(thicknessvariables)] # TODO
    conicvariables = in(:OptimizeConic, propertynames(prescription)) ? prescription[prescription[!, :OptimizeConic], :Conic] : []
    coniclower = [typemin(T) for _ in 1:length(radiusvariables)]
    conicupper = [typemax(T) for _ in 1:length(radiusvariables)]

    #convert to Vector{T} because may have missing values and would otherwise get Array{Union{Missing,Float64}} which optimization routines won't like.
    return Vector{T}(vcat(radiusvariables, thicknessvariables, conicvariables)), Vector{T}(vcat(radiuslower, thicknesslower, coniclower)), Vector{T}(vcat(radiusupper, thicknessupper, conicupper))
end

"""
    updateoptimizationvariables(a::AxisymmetricOpticalSystem{T}, optimizationvariables::Vector{S}) -> AxisymmetricOpticalSystem{S}

Creates a new optical system with updated variables corresponding to the optimization variables.
"""
function updateoptimizationvariables(a::AxisymmetricOpticalSystem{T}, optimizationvariables::Vector{S}) where {T<:Real,S<:Real}
    prescription = a.prescription
    #assume these two types of prescription variables will always be present
    prescradii = prescription[!, :Radius]
    prescthicknesses = prescription[!, :Thickness]

    @assert reduce(&, map((x) -> !isnan(x), optimizationvariables)) #make sure no optimization variables are NaN

    conic = in(:Conic, propertynames(prescription)) ? prescription[!, :Conic] : nothing

    radii = Vector{Union{Missing,S}}(undef, 0)
    thicknesses = Vector{Union{Missing,S}}(undef, 0)
    newconic = Vector{Union{Missing,S}}(undef, 0)

    optindex = 1
    if in(:OptimizeRadius, propertynames(prescription))
        radiusflags = prescription[!, :OptimizeRadius]
        for index in eachindex(radiusflags)
            if radiusflags[index]
                push!(radii, optimizationvariables[optindex])
                optindex += 1
            else
                push!(radii, prescradii[index])
            end
        end
    end

    if in(:OptimizeThickness, propertynames(prescription))
        thicknessflags = prescription[!, :OptimizeThickness]
        for index in eachindex(thicknessflags)
            if thicknessflags[index]
                push!(thicknesses, optimizationvariables[optindex])
                optindex += 1
            else
                push!(thicknesses, prescthicknesses[index])
            end
        end
    end

    if conic !== nothing
        if in(:OptimizeConic, propertynames(prescription))
            conicflags = prescription[!, :OptimizeConic]
            for index in eachindex(conicflags)
                if conicflags[index]
                    push!(newconic, optimizationvariables[optindex])
                    optindex += 1
                else
                    push!(newconic, conic[index])
                end
            end
        end
    end

    newprescription = copy(prescription)
    newprescription.Radius = radii
    newprescription.Thickness = thicknesses
    if conic !== nothing
        newprescription.Conic = newconic
    end

    return AxisymmetricOpticalSystem{S}(newprescription, detectorsize(a)..., S, temperature = S(temperature(a)), pressure = S(pressure(a)))
end

end #module
export Optimization


text:
include("../../conf.jl")

mutable struct StatsDocument
    parts::Dict{String,String}
    order::Array{String,1}
    css::String
end
StatsDocument() = StatsDocument(Dict(), [], "")

function appendGraph(doc::StatsDocument, path::String, titles::Array{Tuple{String,String,String},1}, description::String, title::String)::StatsDocument
    if MAKE_OVERVIEW
        if "overview" ∉ keys(doc.parts)
            doc.parts["overview"] = ""
            push!(doc.order, "overview")
        end
        doc.parts["overview"] *= ""
    end

    if title ∉ keys(doc.parts)
        doc.parts[title] = ""
        push!(doc.order, title)
    end
    doc.parts[title] *= "
" * replace("$description", "\n"=>"\n
") *
            "

\n" *
            reduce(*, map(x->"
$(x[1]) ($(x[2])) $(x[3])
\n", titles)) * "

"

    return doc
end

function generateHTML(doc::StatsDocument)::String
    return """

<meta http-equiv="content-type" content="text/html; charset=UTF-8">

$(reduce(*, map(x->"
$(x)

$(doc.parts[x])
", doc.order)))

"""
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
export stack, stack!, remove_nan, remove_nan!, shorten, shorten!
export pws, pws!, robuststack, robuststack!, adaptive_filter!, adaptive_filter
export robustpws, robustpws!
"""

  stack!(C)

Stack correlation by time interval.
The default is to stack by day. Using `allstack == true` will stack all available
correlations. To use phase-weighted stack, specify the amount of
`phase_smoothing` in seconds.

# Arguments
- `C::CorrData`: Correlation data.
- `interval::Period`: Interval over which to stack `C`.
- `allstack::Bool`: If `true`, stack all data.
- `stack_type::Function`: Type of stacking. Options are mean, pws, robuststack, etc..
"""
function stack!(C::CorrData; interval::Period=Day(1),
                allstack::Bool=false,stacktype::Function=mean)

    if allstack == true
          if stacktype == mean
                C.corr = mean(C.corr,dims=2)
          else # use phase-weighted stack
                C.corr = stacktype(C.corr)
          end
          C.t = [C.t[1]]
    else # stack by interval
          C.t = d2u.(round.(u2d.(C.t),interval,RoundDown))
          stackT = unique(C.t)
          ind = indexin(stackT,C.t)
          push!(ind,length(C.t)+1)
          stack_out = Array{eltype(C.corr)}(undef,size(C.corr,1),length(stackT))

          for ii = 1:length(stackT)
                 if stacktype == mean
                       stack_out[:,ii] = mean(C.corr[:,ind[ii]:ind[ii+1]-1],dims=2)
                 else
                       stack_out[:,ii] = stacktype(C.corr[:,ind[ii]:ind[ii+1]-1])
                 end
          end

          C.t = stackT
          C.corr = stack_out
    end
    return C
end
stack(C::CorrData; interval::Period=Day(1),
      allstack::Bool=false,stacktype::Function=mean) = (U = deepcopy(C);
      stack!(U,interval=interval,allstack=allstack,
             stacktype=stacktype);return U)

"""

  robuststack(A)

Performs robust stack on array `A`.

Follows methods of Pavlis and Vernon, 2010.

# Arguments
- `A::AbstractArray`: Time series stored in columns.
- `ϵ::AbstractFloat`: Threshold for convergence of robust stack.
- `maxiter::Int`: Maximum number of iterations to converge to robust stack.
"""
function robuststack(A::AbstractArray{T};ϵ::AbstractFloat=Float32(1e-4),
                     maxiter::Int=10) where T <: AbstractFloat
    N = size(A,2)
    Bold = median(A,dims=2)
    Bold ./= norm(Bold,2)
    w = Array{T}(undef,N)
    r = Array{T}(undef,N)
    d2 = Array{T}(undef,N)

    # do 2-norm for all columns in A
    for ii = 1:N
        d2[ii] = norm(A[:,ii],2)
    end

    BdotD = sum(A .* Bold,dims=1)

    for ii = 1:N
        r[ii] = norm(A[:,ii] .- (BdotD[ii] .* Bold),2)
        w[ii] = abs(BdotD[ii]) ./ d2[ii] ./ r[ii]
    end

    Bnew = mean(A,weights(w),dims=2)
    Bnew ./= norm(Bnew,2)

    # check convergence
    ϵN = norm(Bnew .- Bold,1) / (norm(Bnew,2) * N)
    Bold = Bnew
    iter = 0
    while (ϵN > ϵ) && (iter <= maxiter)
        BdotD = sum(A .* Bold,dims=1)

        for ii = 1:N
            r[ii] = norm(A[:,ii] .- (BdotD[ii] .* Bold),2)
            w[ii] = abs(BdotD[ii]) ./ d2[ii] ./ r[ii]
        end

        Bnew = mean(A,weights(w),dims=2)
        Bnew ./= norm(Bnew,2)

        # check convergence
        ϵN = norm(Bnew .- Bold,1) / (norm(Bnew,2) * N)
        Bold = Bnew
        iter += 1
    end
    return Bnew
end
robuststack!(C::CorrData;ϵ::AbstractFloat=eps(Float32),maxiter::Int=10) =
       (C.corr = robuststack(C.corr,ϵ=ϵ,maxiter=maxiter); C.t = C.t[1:1]; return C)
robuststack(C::CorrData;ϵ::AbstractFloat=eps(Float32),maxiter::Int=10)  =
       (U = deepcopy(C); U.corr = robuststack(U.corr,ϵ=ϵ,maxiter=maxiter); U.t = U.t[1:1];
       return U)

"""

  pws(A)

Performs phase-weighted stack on array `A` of time series.

Follows methods of Schimmel and Paulssen, 1997.
If s(t) is time series data,
S(t) = s(t) + i*H(s(t)), where H(s(t)) is Hilbert transform of s(t)
S(t) = s(t) + i*H(s(t)) = A(t)*exp(i*phi(t)), where
A(t) is envelope of s(t) and phi(t) is phase of s(t)
Phase-weighted stack, g(t), is then:
g(t) = 1/N sum j = 1:N s_j(t) * | 1/N sum k = 1:N exp[i * phi_k(t)]|^v
where N is number of traces used, v is sharpness of phase-weighted stack

# Arguments
- `A::AbstractArray`: Time series stored in columns.
- `pow::Int`: Sharpness of transition from phase similarity to dissimilarity.
"""
function pws(A::AbstractArray{T}; pow::Real=2) where T <: AbstractFloat
    # preserve type-stability
    if !isa(pow,Int)
        pow = T(pow)
    end
    Nrows,Ncols = size(A)
    phase_stack = abs.(sum(exp.(im .* angle.(hilbert(A))),dims=2) ./ Ncols) .^ pow
    return mean(A .* phase_stack,dims=2)
end
pws!(C::CorrData; pow::Real=2) = (C.corr=pws(C.corr,pow=pow); C.t = C.t[1:1]; return nothing)
pws(C::CorrData; pow::Real=2) = (U = deepcopy(C);
    U.corr=pws(U.corr,pow=pow); U.t = U.t[1:1];return U)

"""

  remove_nan!(C)

Remove correlations with Nan from CorrData `C`.
"""
function remove_nan!(C::CorrData)
      ind = []
      for ii = 1:length(C.t)
            if !any(isnan,C.corr[:,ii])
                  append!(ind,ii)
            end
      end

      # throw error if all values are NaN
      if length(ind) == 0
          throw(ArgumentError("All correlation windows contain NaNs."))
      end

      # return non-NaN columns
      C.corr = C.corr[:,ind]
      C.t = C.t[ind]
      return nothing
end
remove_nan(C::CorrData) = (U = deepcopy(C);remove_nan!(U);return U)

"""

  shorten!(C,newlag)

Clip CorrData `C` from lags [in s] τ = -`newlag` to `newlag`.
"""
function shorten!(C::CorrData, newlag::Real)
    @assert newlag < C.maxlag && newlag > 0 "newlag must satisfy 0 < newlag < C.maxlag"

    # get timearray
    lags = -C.maxlag:1/C.fs:C.maxlag
    ind = findall(x -> abs(x) <= newlag, lags)
    C.maxlag = Float64(newlag)
    C.corr = C.corr[ind,:]
    return nothing
end
shorten(C::CorrData,newlag::Real) = (U = deepcopy(C); shorten!(U,newlag);return U)

"""

  adaptive_filter!(A,window,fs)

Adaptive covariance filter to enhance coherent signals.

Fellows the method of Nakata et al., 2015 (Appendix B). The filtered signal
`X` is given by `X = irfft(P*X(w))`` where `X` is the FFT'ed spectra of `A`
and `P` is the filter. `P` is constructed by using the temporal covariance matrix.

# Arguments
- `A::AbstractArray`: Array of daily/hourly cross-correlation functions
- `window::AbstractFloat`: Window length to apply adaptive filter [s].
- `fs::Real`: Sampling rate of time series `A`.
- `g::Int`: Positive number to adjust the filter harshness
- `overlap::AbstractFloat`: Percent overlap between windows
"""
function adaptive_filter!(A::AbstractArray{T}, window::AbstractFloat,
                          fs::Real; g::AbstractFloat=2., overlap::AbstractFloat=0.9) where T <: AbstractFloat
    if ndims(A) == 1
        return nothing
    end

    if overlap <= 0 || overlap >= 1
        throw(ArgumentError("overlap must be > 0 and < 1."))
    end

    Nrows, Ncols = size(A)
    window_samples = convert(Int,round(window * fs,digits=0))
    window_step = convert(Int,round(window_samples * (1. - overlap),digits=0))
    minind = 1:window_step:Nrows - window_samples
    overlap_factor = T(round(1. - overlap,digits=3))

    # allocate out array
    Aout = zeros(T,size(A))

    # loop through each window
    for ii in eachindex(minind)
        Ain = taper(A[minind[ii]:minind[ii]+window_samples-1,:],fs,
                    max_percentage=T(overlap_factor/2))
        Aout[minind[ii]:minind[ii]+window_samples-1,:] .+= ACF_kernel(@view(Ain[:,:]),g=g) .* overlap_factor
    end

    # windows at right edge
    reverseind = Nrows:-window_step:Nrows-window_samples
    for ii in eachindex(reverseind)
        Ain = taper(A[reverseind[ii]-window_samples+1:reverseind[ii],:],fs,
                    max_percentage=T(overlap_factor/2))
        Aout[reverseind[ii]-window_samples+1:reverseind[ii],:] .+= ACF_kernel(@view(Ain[:,:]),g=g) .* overlap_factor
    end

    copyto!(A,Aout)
    return nothing
end

adaptive_filter(A::AbstractArray{T}, window::AbstractFloat, fs::Real; g::AbstractFloat=2.,
                overlap::AbstractFloat=0.9) where T <: AbstractFloat = (U = deepcopy(A);
                adaptive_filter!(U,window,fs,g=g,overlap=overlap);
                return U)

adaptive_filter!(C::CorrData, window::AbstractFloat; g::AbstractFloat=2.,
                 overlap::AbstractFloat=0.9) = (adaptive_filter!(C.corr,
                 window,C.fs,g=g,overlap=overlap); return nothing)
adaptive_filter(C::CorrData, window::AbstractFloat; g::AbstractFloat=2.,
                overlap::AbstractFloat=0.9) = (U = deepcopy(C);
                adaptive_filter!(U.corr,window,U.fs,g=g,overlap=overlap);
                return U)

function ACF_kernel(A::AbstractArray{T}; g::AbstractFloat=2.) where T <: AbstractFloat
    Nrows, Ncols = size(A)
    g = T(g)
    # fft the 2D array
    spec = rfft(A,1)

    # create auto-covariance function
    Nspec = size(spec,1)
    S1 = zeros(complex(T),Nspec)
    S2 = zeros(complex(T),Nspec)
    for ii = 1:Ncols
        for jj = 1:Nspec
            S2[jj] += spec[jj,ii] .* conj(spec[jj,ii])
        end
    end

    for ii = 1:Ncols
        for jj = 1:Ncols
            if jj != ii
                for kk = 1:Nspec
                    S1[kk] += spec[kk,ii] .* conj(spec[kk,jj])
                end
            end
        end
    end

    # construct filter p
    p = ((S1 .- S2) ./ (S2 .* (Ncols -1))).^g

    # make ifft
    spec .*= p
    Aout = irfft(spec,Nrows,1)
    return Aout
end

"""

  robustpws(A)

Performs combined robust and phase-weighted stack on time series `A`.

Follows methods of Pavlis and Vernon, 2010 and Schimmel and Paulssen, 1997.
This method improves stacking by using the global convergence of the robust
stack and the phase convergence of the phase-weighted stack.
The robust stack downweights outlier correlations, while the phase-weighted
stack downweights outlier phases.

# Arguments
- `A::AbstractArray`: Time series stored in columns.
- `ϵ::AbstractFloat`: Threshold for convergence of robust stack.
- `maxiter::Int`: Maximum number of iterations to converge to robust stack.
- `pow::Int`: Sharpness of transition from phase similarity to dissimilarity.
"""
function robustpws(A::AbstractArray{T};ϵ::AbstractFloat=Float32(1e-6),
                     maxiter::Int=10,pow::Real=2) where T <: AbstractFloat
    N = size(A,2)
    Bold = median(A,dims=2)
    w = Array{T}(undef,N)
    r = Array{T}(undef,N)
    d2 = Array{T}(undef,N)

    # do 2-norm for all columns in A
    for ii = 1:N
        d2[ii] = norm(A[:,ii],2)
    end

    BdotD = sum(A .* Bold,dims=1)

    for ii = 1:N
        r[ii] = norm(A[:,ii] .- (BdotD[ii] .* Bold),2)
        w[ii] = abs(BdotD[ii]) ./ d2[ii] ./ r[ii]
    end
    w ./= sum(w)

    Bnew = mean(A,weights(w),dims=2)

    # check convergence
    ϵN = norm(Bnew .- Bold,2) / (norm(Bnew,2) * N)
    Bold = Bnew
    iter = 0
    while (ϵN > ϵ) && (iter <= maxiter)
        BdotD = sum(A .* Bold,dims=1)

        for ii = 1:N
            r[ii] = norm(A[:,ii] .- (BdotD[ii] .* Bold),2)
            w[ii] = abs(BdotD[ii]) ./ d2[ii] ./ r[ii]
        end
        w ./= sum(w)

        Bnew = mean(A,weights(w),dims=2)

        # check convergence
        ϵN = norm(Bnew .- Bold,2) / (norm(Bnew,2) * N)
        Bold = Bnew
        iter += 1
    end

    # multiply by weights
    W = A .* w'

    # return the phase weighted stack
    return pws(W,pow=pow)
end

robustpws!(C::CorrData; ϵ::AbstractFloat=Float32(1e-6),
           maxiter::Int=10,pow::Int=2) = (C.corr=robustpws(C.corr,ϵ=ϵ,
           maxiter=maxiter,pow=pow); C.t = C.t[1:1]; return nothing)
robustpws(C::CorrData; ϵ::AbstractFloat=Float32(1e-6),
           maxiter::Int=10,pow::Int=2) = (U = deepcopy(C);
           U.corr=robustpws(U.corr,ϵ=ϵ,
           maxiter=maxiter,pow=pow); U.t = U.t[1:1];return U)


text:
using NoRainbows
using Documenter

DocMeta.setdocmeta!(NoRainbows, :DocTestSetup, :(using NoRainbows); recursive=true)

makedocs(;
    modules=[NoRainbows],
    authors="Andrey Oskin",
    repo="https://github.com/Arkoniak/NoRainbows.jl/blob/{commit}{path}#{line}",
    sitename="NoRainbows.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Arkoniak.github.io/NoRainbows.jl",
        siteurl="https://github.com/Arkoniak/NoRainbows.jl",
    ),
    pages=[
        "Home" => "index.md",
        "Colors" => "colors.md",
        "Colorschemes" => "colorschemes.md",
        "Reverse stacktrace" => "reverse.md",
        "Frameline format" => "frameline.md",
    ],
)

deploydocs(;
    repo="github.com/Arkoniak/NoRainbows.jl",
)


text:
export aperture_photometry

@doc raw"""
    aperture_photometry(channel::Array,
                                 obj::Tuple,
                                 ref1::Tuple,
                                 ref2::Tuple,
                                 radius::Tuple
                                )

Aperture photometry of a given object
"""
function aperture_photometry(channel::Array,
                             obj::Tuple,
                             ref1::Tuple,
                             ref2::Tuple,
                             radius::Tuple
                            )

    f_obj  = Array{Float64, 1}(undef, 0)
    f_ref1 = Array{Float64, 1}(undef, 0)
    f_ref2 = Array{Float64, 1}(undef, 0)

    for i in 1:length(channel)

        flux_o_in = Float64(sum(channel[i][obj[1]-radius[1]:obj[1]+radius[1], obj[2]-radius[1]:obj[2]+radius[1]]))
        flux_o_out = Float64(sum(channel[i][obj[1]-radius[2]:obj[1]+radius[2], obj[2]-radius[2]:obj[2]+radius[2]]))

        bg_o = (flux_o_out - flux_o_in) / (radius[2]^2) # bg / pix
        flux_o = flux_o_in - (bg_o * (radius[1]^2))

        flux_r1_in = Float64(sum(channel[i][ref1[1]-radius[1]:ref1[1]+radius[1], ref1[2]-radius[1]:ref1[2]+radius[1]]))
        flux_r1_out = Float64(sum(channel[i][ref1[1]-radius[2]:ref1[1]+radius[2], ref1[2]-radius[2]:ref1[2]+radius[2]]))

        bg_r1 = (flux_r1_out - flux_r1_in) / (radius[2]^2)
        flux_r1 = flux_r1_in - (bg_r1 * (radius[1]^2))

        flux_r2_in = Float64(sum(channel[i][ref2[1]-radius[1]:ref2[1]+radius[1], ref2[2]-radius[1]:ref2[2]+radius[1]]))
        flux_r2_out = Float64(sum(channel[i][ref2[1]-radius[2]:ref2[1]+radius[2], ref2[2]-radius[2]:ref2[2]+radius[2]]))

        bg_r2 = (flux_r2_out - flux_r2_in) / (radius[2]^2)
        flux_r2 = flux_r2_in - (bg_r2 * (radius[1]^2))

        f_obj  = append!(f_obj, flux_o)
        f_ref1 = append!(f_ref1, flux_r1)
        f_ref2 = append!(f_ref2, flux_r2)
    end

    return f_obj, f_ref1, f_ref2

end


text:
export Stiefel

struct Stiefel{T <: Real} <: ConvexSet
    r::T
end
Stiefel() = Stiefel(0.0)

function project!(s::Stiefel{<:Real}, v::Matrix{T}, y::Matrix{T}) where {T <: Real}
    U, S, V = LinearAlgebra.svd(y)
    v = U*V'
end

function project(s::Stiefel{<:Real}, y::Matrix{<:Real})
    project!(s, similar(y), y)
end


text:
function get_current_volume(config)
    current_volume = typeof(config["total_value"]) == Float64 ? config["total_value"] : parse(Float64, config["total_value"])
    if length(config["extra_volume"]) > 0
        for (k,v) in config["extra_volume"]
            current_volume += v
        end
    end
    current_volume += config["thes_volume"]
    return current_volume
end
export get_current_volume


text:
Normal(d::Dict) = Normal(d[:μ], d[:σ])
MixtureModel(d::Dict) = MixtureModel(map(Normal, d[:components]), Vector{Float64}(d[:prior][:p]))
# HMMBase.HMM(d::Dict, T::Type) = HMM(Vector{Float64}(d[:a]), hcat(d[:A]...), map(T, d[:B]))

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk


text:
@test 1   |> Posit{8,0}  |> SigmoidNumbers.find_lsb ==  0
@test 1   |> Posit{32,2} |> SigmoidNumbers.find_lsb ==  0
@test 5   |> Posit{8,0}  |> SigmoidNumbers.find_lsb ==  0
@test 9   |> Posit{32,2} |> SigmoidNumbers.find_lsb ==  0
@test 0.5 |> Posit{8,0}  |> SigmoidNumbers.find_lsb == -1
@test 0.5 |> Posit{32,2} |> SigmoidNumbers.find_lsb == -1
@test 1.5 |> Posit{8,0}  |> SigmoidNumbers.find_lsb == -1
@test 1.5 |> Posit{32,2} |> SigmoidNumbers.find_lsb == -1

@test -1   |> Posit{8,0}  |> SigmoidNumbers.find_lsb ==  0
@test -1   |> Posit{32,2} |> SigmoidNumbers.find_lsb ==  0
@test -5   |> Posit{8,0}  |> SigmoidNumbers.find_lsb ==  0
@test -9   |> Posit{32,2} |> SigmoidNumbers.find_lsb ==  0
@test -0.5 |> Posit{8,0}  |> SigmoidNumbers.find_lsb == -1
@test -0.5 |> Posit{32,2} |> SigmoidNumbers.find_lsb == -1
@test -1.5 |> Posit{8,0}  |> SigmoidNumbers.find_lsb == -1
@test -1.5 |> Posit{32,2} |> SigmoidNumbers.find_lsb == -1

function find_lsb_q{N,ES}(x::Posit{N,ES})
  q = Quire(Posit{N,ES})
  set!(q, x)
  SigmoidNumbers.find_lsb(q)
end

@test 1   |> Posit{8,0}  |> find_lsb_q ==  0
@test 1   |> Posit{32,2} |> find_lsb_q ==  0
@test 5   |> Posit{8,0}  |> find_lsb_q ==  0
@test 9   |> Posit{32,2} |> find_lsb_q ==  0
@test 0.5 |> Posit{8,0}  |> find_lsb_q == -1
@test 0.5 |> Posit{32,2} |> find_lsb_q == -1
@test 1.5 |> Posit{8,0}  |> find_lsb_q == -1
@test 1.5 |> Posit{32,2} |> find_lsb_q == -1

@test -1   |> Posit{8,0}  |> find_lsb_q ==  0
@test -1   |> Posit{32,2} |> find_lsb_q ==  0
@test -5   |> Posit{8,0}  |> find_lsb_q ==  0
@test -9   |> Posit{32,2} |> find_lsb_q ==  0
@test -0.5 |> Posit{8,0}  |> find_lsb_q == -1
@test -0.5 |> Posit{32,2} |> find_lsb_q == -1
@test -1.5 |> Posit{8,0}  |> find_lsb_q == -1
@test -1.5 |> Posit{32,2} |> find_lsb_q == -1

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
import IntervalArithmetic: Interval, mid
import SpecialFunctions: erfc

# Monkey-patch a few functions for Intervals
# ... this is far from proper and a bit specific for our use case here
# (that's why it's not contributed upstream).
# should be done e.g. by changing  the rounding mode ...
erfc(i::Interval) = Interval(prevfloat(erfc(i.lo)), nextfloat(erfc(i.hi)))

function compute_fft_size(lattice::AbstractMatrix{T}, Ecut; kwargs...) where T <: Interval
    # This is done to avoid a call like ceil(Int, ::Interval)
    # in the above implementation of compute_fft_size,
    # where it is in general cases not clear, what to do.
    # In this case we just want a reasonable number for Gmax,
    # so replacing the intervals in the lattice with
    # their midpoints should be good.
    compute_fft_size(mid.(lattice), Ecut; kwargs...)
end

function local_potential_fourier(el::ElementCohenBergstresser, q::T) where {T <: Interval}
    lor = round(q.lo, digits=5)
    hir = round(q.hi, digits=5)
    @assert iszero(round(lor - hir, digits=3))
    T(local_potential_fourier(el, mid(q)))
end


text:
using Unitful, OrdinaryDiffEq, Test
f(du,u,p,t) = du .= 3*u"1/s"*u
prob = ODEProblem(f, [2.0u"m"], (0.0u"s", Inf*u"s"))
intg = init(prob, Tsit5())
@test_nowarn step!(intg, 0.02u"s", true)


text:
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/common_functions.jl")
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/FK_sampling_functions.jl")
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/SB_sampling_functions.jl")

using JLD, GibbsTypePriors, DataFramesMeta, RCall, Optim, OptimTestProblems, StatsBase

## functions to generate the plots

function approximation_prior_distribution(beta,sigma,N,Nt,sigma_arr)
    Pkn_numeric_ = GibbsTypePriors.Pkn_NGG.(1:N, N, beta, sigma)
    Pkn_order2_ = Pkn_NGG_pred_approx(N, beta, sigma)
    Pkn_NGGMult = Pkn_NGGM_precomp.(1:N,N,Nt,beta,sigma, [Pkn_numeric_])
    Pkn_FK_ =  Pkn_NGG_FK(N, beta, sigma, Nt; runs=2*10^2)
    #R"p_ngg_sb <- Prior_on_K_SB_NGG($beta, $sigma, $N,$(Nt), runs=2*10^2)"
    #Pkn_SB_raw  = rcopy(R"p_ngg_sb$pk")
    Pkn_SB_  = Pkn_NGG_SB_1(N,beta,sigma,Nt; runs=2*10^2)

    df= DataFrame(Pkn_numeric = Pkn_numeric_,
                  Pkn_order2 = Pkn_order2_,
                  Pkn_NGGM = Pkn_NGGMult,
                  Pkn_FK = Pkn_FK_,
                  Pkn_SB = Pkn_SB_)
       return df
end

###########

n=100
β= 1.0
ntr=250
sigma_vec= [0.25,0.75]
DF_all_100 = map(x ->approximation_prior_distribution(β,x,n,ntr,sigma_vec),sigma_vec)

n=100
β= 10.0
ntr=250
sigma_vec= [0.25,0.75]
DF_all_100_10 = map(x ->approximation_prior_distribution(β,x,n,ntr,sigma_vec),sigma_vec)

DF_100 = [DF_all_100,DF_all_100_10 ]

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_1.jld", "DF_100" ,DF_100)

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
using Test
using Cleaner: CleanTable, add_index

@testset "add_index is working as expected" begin
    testCT = CleanTable([:A], [[4, 5, 6]])

    @test add_index(testCT).row_index == [1, 2, 3]
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("colorscheme.jl"))
using Makie, DynamicalSystems, InteractiveChaos
using AbstractPlotting.MakieLayout

# %% Animation of changing ε in curve of 1D energy balance model
function fitzhugh(x, p, t)
    u, w = x
    I, a, b = p
    udot = a*u*(u-b)*(1-u) - w + I
    wdot = 0.01*(u-w)
    return SVector(udot, wdot)
end

a = 3.0
b = 0.2
I = 0.0
p = [I, a, b]
fh = ContinuousDynamicalSystem(fitzhugh, rand(2), p)

# Layout plot:
scene, layout = layoutscene(resolution = (1000, 800))
ax = layout[1,1] = LAxis(scene)
display(scene)
ax.xlabel = "u"
ax.ylabel = "w"

# plot the two nullclines
us = -0.4:0.01:1.2
w1 = us
w2 = @. a*us*(us-b)*(1-us) + I
lines!(ax, us, w1; color = COLORS[1], linewidth = 2, linestyle = :dash)
lines!(ax, us, w2; color = COLORS[2], linewidth = 2, linestyle = :dot)

spoint = select_point(ax.scene)
on(spoint) do pos
    tr = trajectory(fh, 10000, SVector(pos...))
    lines!(ax.scene, columns(tr)...; color = InteractiveChaos.randomcolor())
end

# %% orbit evolution in a 2D torus
using OrdinaryDiffEq

R = 2.0
r = 0.5

function torus(u)
    θ, φ = u
    x = (R + r*cos(θ))*cos(φ)
    y = (R + r*cos(θ))*sin(φ)
    z = r*sin(θ)
    return SVector(x, y, z)
end

function quasiperiodic_f(u, p, t)
    # here we make the frequency ratio a state variable, because the
    # trajectory_animator application doesn't allow different parameters for different
    # initial conditions
    θ, φ, ω = u
    θdot = ω
    φdot = 1.0
    return SVector(θdot, φdot, 0.0)
end

ds = ContinuousDynamicalSystem(quasiperiodic_f, [0.0, 0.0, 0.0], nothing)
u0s = [[0, 0, ω] for ω ∈ [3, sqrt(7), 5]]
diffeq = (alg = Tsit5(), dtmax = 0.01)

lims = ((-R-r, R+r), (-R-r, R+r), (-3r, 3r))

scene, main, layout, obs = interactive_evolution(
    ds, u0s; tail = 20000, diffeq, colors = COLORS, transform = torus, lims,
    plotkwargs = (linewidth = 1.0,),
)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
module ProjectAssistant

using Pkg, Dates, Logging, REPL.TerminalMenus, LibGit2
using PkgTemplates

export project_initialize, project_activate

##################################################

dirs_project = [
    "tmp",
    "logs",
    "assets",
    "data",
    "papers",
    "scripts",
    "project",
]

# DrWatson dirs, see https://juliadynamics.github.io/DrWatson.jl/dev/project/
dirs_DrWatson = [
    "_research",
    "_research/tmp",
    "data/exp_pro",
    "data/exp_raw",
    "plots",
    "notebooks",
    "papers",
]

files_project = [
  ("REQUIRE", ""),
  ("ENV", ""),
  ("project/startup.jl", ""),
  ("project.bat","@echo off\n%LOCALAPPDATA%/Julia-$(string(VERSION))/bin/julia.exe --color=yes project/project.jl\n"*"pause\n")
]

copy_files_project = [
  "project/run.jl",
  "project/project.jl",
]

##################################################

printError(ex) = @error sprint(showerror, ex, catch_backtrace())

# if VERSION.major > 1 || (VERSION.major == 1 && VERSION.minor > 3)
#     function installed() # since Julia 1.4: Pkg.installed() is deprecated
#         deps = Pkg.dependencies()
#         installs = Dict{String, VersionNumber}()
#         for (uuid, dep) in deps
#             dep.is_direct_dep || continue
#             dep.version === nothing && continue
#             installs[dep.name] = dep.version
#         end
#         return installs
#     end
# else
#     installed() = Pkg.installed()
# end

# will change when project is updated to Julia 1.4
installed() = Pkg.installed()

is_pkgs_missing(pkgs::AbstractArray) = length(filter(x->!haskey(installed(),x),pkgs))>0

  # update when some packages are not added
function install_missing_pkgs(pkgs::AbstractArray)
    missing = filter(x->!haskey(installed(),x),pkgs)
    if length(missing)>0
        for pkg in missing Pkg.add(pkg) end
        Pkg.update()
    end
end

##################################################

function project_initialize()
    PATH_OLD = pwd() # root dir (all projects)
    PROJECT = length(ARGS) > 0 ? ARGS[1] : ""
    USERNAME = length(ARGS) > 1 ? ARGS[2] : ""
    PATH_PROJECT = nothing
    SKIP = true
    TEMPLATE = nothing
    DrWatsonDirs = false

    LOGGER = ConsoleLogger(stdout, Logging.Debug)

    with_logger(LOGGER) do
        while true
            try
                if isempty(PROJECT)
                    print("New project name: "); PROJECT = readline()
                    if isempty(PROJECT) throw(ArgumentError("Project name is required")) end

                    print("Skip questions and use default settings? [yes]: "); input = readline()
                    if input == "n" || input == "no" SKIP=false end
                end

                if !SKIP
                    TEMPLATE = interactive_template()

                    print("Create additional DrWatson dirs? [no]: "); input = readline()
                    if input == "y" || input == "yes" DrWatsonDirs=true end
                else
                    if isempty(USERNAME)
                        USERNAME = LibGit2.getconfig("github.user", "")
                        if isempty(USERNAME) print("Username: "); USERNAME = readline() end
                        if isempty(USERNAME) throw(ArgumentError("Username is required")) end
                    end
                end

            catch ex
                printError(ex)
                print("Retry? [yes]: "); input = readline()
                if input == "n" || input == "no" break end
                @info "Retry..."
                continue
            end

            break
        end
    end

    if isnothing(TEMPLATE) # default
        TEMPLATE = Template(;
            user=USERNAME,
            #host="https://github.com",
            #license="MIT",
            authors=[USERNAME],
            #dir=Pkg.devdir() #abspath(pwd())
            #julia_version=VERSION,
            #ssh=false,
            dev=false, # true
            #manifest=false,
            plugins=[
                GitHubPages(),
                AppVeyor(),
                Codecov(), #; config_file=nothing
                #Coveralls(),
                TravisCI()
            ],
            #git=true
        )
    end

    generate(TEMPLATE, PROJECT)
    PATH_PROJECT = abspath(TEMPLATE.dir,PROJECT)

    if !isdir(PATH_PROJECT) throw(ArgumentError("Cannot write dir $(PATH_PROJECT)!")) end

    cd(PATH_PROJECT) # go to project dir

    path_project_files = abspath(PATH_PROJECT,"project")

    if !isdir(path_project_files)
        @info "create dir project."
        mkdir(path_project_files)
    end

    # write dirs
    dir_lists = [dirs_project]
    if DrWatsonDirs push!(dir_lists,dirs_DrWatson) end

    for dirs in dir_lists
        for name in dirs
            path_dir = abspath(name)
            if !isdir(path_dir)
                @info "create dir $(basename(name))."
                mkdir(path_dir)
            end
        end
    end

    git_files = String[]

    # write files
    for file in files_project
        name = file[1]
        content = file[2]
        path_file = abspath(name)
        if !isfile(path_file)
            @info "create file $(basename(name))."
            write(path_file, content)
            push!(git_files,path_file)
        end
    end

    # copy files
    for name in copy_files_project
        path_file = abspath(name)
        if !isfile(path_file)
            @info "copy file $(basename(name))."
            cp(abspath(@__DIR__,"../",name), path_file)
            push!(git_files,path_file)
        end
    end

    # TODO: add files to existing git repo...
    # repo = LibGit2.init(PATH_PROJECT) # wrong?
    # LibGit2.add!(repo,git_files...)
end

function project_activate(dir=pwd())

    cd(dir) # go to project dir
    Pkg.activate(pwd()) # activate this project
    ENV_backup = ENV # backup ENV

    ##################################################

    function show_menu()
        options = ["info", "start", "test", "make docs", "update packages", "exit"]
        menu = RadioMenu(options) #, pagesize=6
        selected = 1

        while true
            write(stdin.buffer, repeat("",max(selected-1,0))) # simulate input

            choice = request("Choose option:", menu)

            if choice != -1
                selected = choice
                option = options[choice]
                println("Choosed: ", option)

                try
                    if choice == length(options) break
                    elseif choice == 1 get_project_info()
                    elseif choice == 2 runProject()
                    elseif choice == 3 Pkg.test()
                    elseif choice == 4 createDocs()
                    elseif choice == 5 update()
                    end
                catch ex
                    printError(ex)
                finally
                    resetENV()
                end
            else
                println("Menu canceled.")
                break
            end
        end
    end

    ##################################################

    function get_project_info()
        file_path = abspath(pwd(),"Project.toml")
            packages = join(sort((x->Pair(x.first,isnothing(x.second) ? "" : string(x.second))).(collect(installed()))),"
        ")

        println()
        @info "Project.toml\n" * read(file_path, String)
        println()
        @info "Installed packages\n" * packages
        println()
    end

    ##################################################

    function resetENV()
        for (k,v) in ENV
            if haskey(ENV_backup, k)
                ENV[k] = ENV_backup[k]
            end
        end
    end

     # set environment variables
    function readsetENV(file::String = abspath(pwd(),"ENV"))
        if !isfile(file) write(file,"") end
        for line in readlines(file)
            len = length(line)
            if iszero(len) || line[1] == '#' continue end
            env=split(replace(line,"\""=>""),'=')
            ENV[env[1]] = length(env)>1 ? env[2] : ""
        end
    end

    ##################################################

    function runProject()
        project = Symbol(basename(pwd()))
        path_project = abspath(pwd(),"project")
        readsetENV()
        include(abspath(path_project,"startup.jl"))
        file_run = abspath(path_project,"run.jl")
        @info "start project \x1b[33m$project\x1b[39m..."
        run(`$(Sys.BINDIR)/julia.exe --color=yes --startup-file=no --project $file_run $ARGS`)
        @info "end project \x1b[33m$project\x1b[39m."
    end

    ##################################################

    function createDocs()
        file_docs_make = abspath(pwd(),"docs","make.jl")
        !isfile(file_docs_make) && error("docs/make.jl does not exist!")
        include(file_docs_make)
    end

    ##################################################

    function update()
        file_path = abspath(pwd(),"REQUIRE")

        pkgs_require = filter(x->isnothing(match(r"^#",x)),readlines(file_path))
        @info "Check update..." pkgs_require

        pkgs_list = [x.first for x in installed()]
        pkgs_install = filter(x->!in(x,pkgs_list),pkgs_require)
        pkgs_remove = filter(x->!in(x,pkgs_require),pkgs_list)

        for pkg_name in pkgs_install
            @info "Install $pkg_name..."

            try
                Pkg.add(pkg_name)
                @info "$pkg_name installed"
            catch e
                dump(e)
            end
        end

        for pkg_name in pkgs_remove
            if in(pkg_name,pkgs_project_essential) continue end # skip if package is marked as "essential"

            @info "Remove $pkg_name..."

            try
                Pkg.rm(pkg_name)
                @info "$pkg_name removed"
            catch e
                dump(e)
            end
        end

        # update when some packages are not added
        #pkgs_essential_missing = filter(x->!haskey(installed(),x),pkgs_project_essential)
        #if length(pkgs_essential_missing)>0
        #   for pkg in pkgs_essential_missing Pkg.add(pkg) end
        #end

        Pkg.update()
    end

    ##################################################

    LOGGER = ConsoleLogger(stdout, Logging.Debug)

    with_logger(LOGGER) do
        show_menu()
    end
end

end #module


text:
module Solitaire

const Card = Int8
const A = Card(53) # large joker
const B = Card(54) # small joker

const Deck = Vector{Card}

"Full deck of cards (with 2 jokers), sorted in increasing order."
sorted_deck() = Deck(1:54)

"Move A one position, with wrap-around."
function move_A(deck::Deck)
    pos = findfirst(isequal(A), deck)
    if pos == length(deck)
        return vcat(deck[1], A, deck[2:end-1])
    else
        return vcat(deck[1:pos-1], deck[pos+1], A, deck[pos+2:end])
    end
end

"Move B two positions, with wrap-around."
function move_B(deck::Deck)
    pos = findfirst(isequal(B), deck)
    if pos == length(deck)
        return vcat(deck[1:2], B, deck[3:end-1])
    elseif pos == length(deck) - 1
        return vcat(deck[1], B, deck[2:pos-1], deck[pos+1:end])
    else
        return vcat(deck[1:pos-1], deck[pos+1:pos+2], B, deck[pos+3:end])
    end
end

"Triple cut around the positions of A,B (or B,A)."
function cut_AB(deck::Deck)
    a = findfirst(isequal(A), deck)
    b = findfirst(isequal(B), deck)
    fst, snd = extrema([a, b])
    left, middle, right = deck[1:fst-1], deck[fst:snd], deck[snd+1:end]
    return vcat(right, middle, left)
end

"Cut after count to value of bottom card."
function cut_bottom(deck::Deck)
    bottom = deck[end]
    if bottom in [A, B]
        bottom = 53 # either joker counts the same
    end
    left, middle, right = deck[1:bottom], deck[bottom+1:end-1], deck[end:end]
    return vcat(middle, left, right)
end

"Extract value for key stream out of current state."
function value(deck::Deck)
    top = deck[1]
    if top in [A, B]
        top = 53 # either joker counts the same
    end
    location = (top + 1)
    val = deck[location]
    if val in [A, B]
        val = 53
    end
    return val
end

"Iterates over deck states."
struct DeckIterator
    state::Deck
end

function Base.iterate(iter::DeckIterator, state::Deck=iter.state)
    # do the steps to advance the state
    state = move_A(state)
    state = move_B(state)
    state = cut_AB(state)
    state = cut_bottom(state)

    # extract value
    v = value(state)

    # never stop iterating!
    return v, state
end

"Cut deck according to given value"
function cut_given(deck::Deck, value::Int)
    @assert 1 ≤ value ≤ length(deck) - 1
    left, middle, right = deck[1:value], deck[value+1:end-1], deck[end:end]
    return vcat(middle, left, right)
end

"Generate deck from pass phrase"
function initialize(passphrase::String)
    # convert passphrase to numerical array
    values = [letter - 'A' + 1 for letter in uppercase(passphrase)]
    if !all(1 .≤ values .≤ 26)
        error("Invalid passphrase: Only letters (A-Z) are supported.")
    end

    # start with sorted deck
    deck = sorted_deck()

    # now go through Solitaire steps, but interleave another cut
    iter = DeckIterator(deck)
    for value in values
        _, deck = iterate(iter, deck)
        deck = cut_given(deck, value)
    end

    return deck
end

"Card is not a joker."
not_joker(card) = card != A && card !=B

"Filter out jokers from key stream."
skip_jokers(cards) = Iterators.filter(not_joker, cards)

"Encrypt cleartext with given deck."
function encrypt(cleartext::String, deck::Deck)
    keygen = skip_jokers(DeckIterator(deck))
    clearnum = [letter - 'A' + 1 for letter in uppercase(cleartext)]
    ciphernum = Int[]
    for t in zip(clearnum, keygen)
        c = t[1] + t[2]
        while c > 26
            c -= 26
        end
        push!(ciphernum, c)
    end
    return String('A' .+ (ciphernum .- 1))
end

"Encrypt cleartext with given passphrase."
function encrypt(cleartext::String, passphrase::String)
    return encrypt(cleartext, initialize(passphrase))
end

"Decrypt ciphertext with given deck."
function decrypt(ciphertext::String, deck::Deck)
    keygen = skip_jokers(DeckIterator(deck))
    ciphernum = [letter - 'A' + 1 for letter in uppercase(ciphertext)]
    clearnum = Int[]
    for t in zip(ciphernum, keygen)
        c = t[1] - t[2] + 52
        while c > 26
            c -= 26
        end
        push!(clearnum, c)
    end
    return String('A' .+ (clearnum .- 1))
end

"Decrypt ciphertext with given passphrase."
function decrypt(ciphertext::String, passphrase::String)
    return decrypt(ciphertext, initialize(passphrase))
end

end


text:
using ApproxOperator
using Test

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

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D)

  # set the freestream flow parameters
  rho_freestream = 1.0
  u_freestream = 0.3
  p_freestream = inv(equations.gamma)

  theta = pi / 90.0 # analogous with a two degree angle of attack
  si, co = sincos(theta)
  v1 = u_freestream * co
  v2 = u_freestream * si

  prim = SVector(rho_freestream, v1, v2, p_freestream)
  return prim2cons(prim, equations)
end

initial_condition = uniform_flow_state

boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
boundary_conditions = Dict( :Bottom => boundary_condition_uniform_flow,
                            :Top    => boundary_condition_uniform_flow,
                            :Right  => boundary_condition_uniform_flow,
                            :Left   => boundary_condition_uniform_flow,
                            :Circle => boundary_condition_slip_wall )

###############################################################################
# Get the DG approximation space

solver = DGSEM(polydeg=4, surface_flux=flux_hll)

###############################################################################
# Get the curved quad mesh from a file

default_mesh_file = joinpath(@__DIR__, "mesh_box_around_circle.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/8b9b11a1eedfa54b215c122c3d17b271/raw/0d2b5d98c87e67a6f384693a8b8e54b4c9fcbf3d/mesh_box_around_circle.mesh",
                                       default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=10,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# based on the contents of align.jl
# consider changing terminology from align to calibrate.
# assumes κs have been updated into As, Bs, Es.

function alignquantificationregion(y_cost::Vector{Complex{T}},
    U_cost,
    P_cost,
    As,
    Bs,
    Es,
    fs,
    SW,
    Δsys_cs,
    a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
    shift_lb::Vector{T},
    shift_ub::Vector{T};
    LS_inds = 1:length(U_cost),
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-3,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    w_lb_default = 1e-3,
    w_ub_default = 1e2,
    β_max_iters = 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T <: Real

    # prepare.
    N_d = sum( getNd(Bs[n]) for n = 1:length(Bs) )
    @assert length(shift_ub) == length(shift_lb) == N_d

    U_rad_cost = U_cost .* (2*π)

    # setup inner optim over β.
    q, updatedfunc, getshiftfunc, N_vars_set,
    run_optim, obj_func_β, E_BLS, w_BLS, b_BLS, updateβfunc, updatewfunc,
    q_β = setupcostnesteddwarpw(Es, Bs, As, fs, SW, LS_inds, U_rad_cost,
        y_cost, Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings;
        β_optim_algorithm = β_optim_algorithm,
        w_lb_default = w_lb_default,
        w_ub_default = w_ub_default,
        β_max_iters = β_max_iters,
        β_xtol_rel = β_xtol_rel,
        β_ftol_rel = β_ftol_rel,
        β_maxtime = β_maxtime)

    # set up outer optim over shifts.
    #N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    #N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

    p_β = zeros(T, N_β) # persistant buffer.

    obj_func = pp->costnesteddw(U_rad_cost, y_cost, updatedfunc, updatewfunc,
    pp, Es, Bs, κs_β_orderings, κs_β_DOFs, run_optim, E_BLS, w_BLS, b_BLS, p_β)

    # optim.
    prob = MultistartOptimization.MinimizationProblem(obj_func, shift_lb, shift_ub)

    local_method = MultistartOptimization.NLoptLocalMethod(; algorithm = local_optim_algorithm,
    xtol_rel = xtol_rel,
    maxeval = maxeval,
    maxtime = maxtime)

    println("N_starts = ", N_starts)

    multistart_method = MultistartOptimization.TikTak(N_starts)
    ret_mo = MultistartOptimization.multistart_minimization(multistart_method,
        local_method, prob)
    #
    println("ret_mo = ", ret_mo)
    minf = ret_mo.value
    minx = ret_mo.location
    ret = :None
    if hasproperty(ret_mo, :ret)
        ret = ret_mo.ret
    end

    # force w_BLS to update.
    obj_func(minx)

    return obj_func, minf, minx, ret, w_BLS
end

"""
Requires user to supply a surrogate.
"""
function alignquantificationcompound(y::Vector{Complex{T}}, U_y, P_y, As, Bs, Es, fs, SW,
    Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings::Vector{Vector{Vector{Int}}},
    shift_lb::Vector{T},
    shift_ub::Vector{T},
    cost_inds_set::Vector{Vector{Int}};
    loop_range = 1:length(cost_inds_set),
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-3,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    w_lb_default = 1e-2,
    w_ub_default = 1e2,
    β_max_iters = 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T <: Real

    N_regions = length(loop_range)
    minxs = Vector{Vector{T}}(undef, N_regions)
    rets = Vector{Symbol}(undef, N_regions)
    minfs = Vector{T}(undef, N_regions)
    obj_funcs = Vector{Function}(undef, N_regions)
    ws = Vector{Vector{T}}(undef, N_regions)

    for r in loop_range
        println("Working on region $(r)")
        obj_funcs[r], minfs[r], minxs[r], rets[r],
        ws[r] = alignquantificationregion(y[cost_inds_set[r]],
            U_y[cost_inds_set[r]],
            P_y[cost_inds_set[r]],
            As,
            Bs,
            Es,
            fs,
            SW,
            Δsys_cs,
            a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
            shift_lb,
            shift_ub;
            N_starts = N_starts,
            local_optim_algorithm = local_optim_algorithm,
            xtol_rel = xtol_rel,
            maxeval = maxeval,
            maxtime = maxtime,
            β_optim_algorithm = β_optim_algorithm,
            w_lb_default = w_lb_default,
            w_ub_default = w_ub_default,
            β_max_iters = β_max_iters,
            β_xtol_rel = β_xtol_rel,
            β_ftol_rel = β_ftol_rel,
            β_maxtime = β_maxtime)
    end

    return obj_funcs, minfs, minxs, rets, ws
end


text:
# TODO PR to UnsafeArrays.jl
const UnsafeVector{T} = UnsafeArray{T,1}
const UnsafeMatrix{T} = UnsafeArray{T,2}

function modelkind(path::AbstractString)
    ext = last(split(last(splitext(path)), '.'))
    if ext in ("xml", "XML", "mjcf", "MJCF")
        return :MJCF
    elseif ext in ("mjb", "MJB")
        return :MJB
    else
        throw(ArgumentError("model path must end with one of [xml, XML, mjcf, MJCF, mjb, MJB]. Got: $ext"))
    end
end

function check_modelpath(path::AbstractString)
    if isfile(path)
        return modelkind(path)
    else
        throw(ArgumentError("model path $path does not exist"))
    end
end

tostring(error::String) = error
tostring(error::Cstring) = unsafe_string(error)
tostring(error::AbstractVector{UInt8}) = String(readuntil(IOBuffer(error), 0x00))

struct MuJoCoException <: Exception
    error::String
    extra::Union{Nothing, String}
    function MuJoCoException(error, extra = nothing)
        new(tostring(error), extra === nothing ? nothing : tostring(extra))
    end
end

Base.show(io::IO, ::MIME"text/plain", ex::MuJoCoException) = show(io, ex)
function Base.show(io::IO, ex::MuJoCoException)
    lines = readlines(IOBuffer(ex.error))
    if ex.extra !== nothing
        println(io, "MuJoCoException(", ex.extra, "):")
        foreach(l->println(io, l), lines)
    elseif length(lines) > 1
        println(io, "MuJoCoException:")
        foreach(l->println(io, l), lines)
    else
        print(io, "MuJoCoException(", first(lines), ')')
    end
end

function MuJoCoException(
    errstr::Union{SVector{N,UInt8},MVector{N,UInt8},Vector{UInt8}};
    checknullterminated::Bool = true,
) where {N}
    errstr = String(errstr)
    if checknullterminated
        nullidx = findfirst('\0', errstr)
        endidx = nullidx isa Integer && 0 <= nullidx <= length(errstr) ? nullidx :
                 length(errstr)
        errstr = errstr[1:endidx]
    end
    MuJoCoException(errstr)
end

macro mjerror(errstr)
    :(throw(MuJoCoException($(esc(errstr)))))
end

macro mjerror(msg, errstr)
    :(throw(MuJoCoException($(esc(msg)), $(esc(errstr)))))
end

macro mjerrormsg(s)
    quote
        local error = MuJoCoException($(esc(s)))
        Base.with_output_color(get(stderr, :color, false) ? Base.error_color() : :nothing, stderr) do io
            showerror(io, error)
        end
    end
end

union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type) = (x,)

macro uninitable(ex)
    sdef = MacroTools.splitstructdef(ex)
    esc(create_constructors(sdef))
end

function create_constructors(sdef)
    length(sdef[:params]) > 0 && error("Parametric types not supported")

    uninitialized_ctor = :($(sdef[:name])() = new())

    fieldexprs, fieldnames = [], []
    for (name, T) in sdef[:fields]
        push!(fieldexprs, :($name::$T))
        push!(fieldnames, name)
    end

    posonly_innerctor_body_typed = Expr(:call, :new, fieldexprs...)
    posonly_innerctor_typed = MacroTools.combinedef(Dict(
        :params => [],
        :name => sdef[:name],
        :args => fieldexprs,
        :kwargs => [],
        :body => posonly_innerctor_body_typed,
        :whereparams => ()
    ))

    posonly_innerctor_body = Expr(:call, :new, fieldnames...)
    posonly_innerctor = MacroTools.combinedef(Dict(
        :params => [],
        :name => sdef[:name],
        :args => fieldnames,
        :kwargs => [],
        :body => posonly_innerctor_body,
        :whereparams => ()
    ))

    push!(sdef[:constructors], uninitialized_ctor)
    push!(sdef[:constructors], posonly_innerctor_typed)
    push!(sdef[:constructors], posonly_innerctor)
    MacroTools.combinestructdef(sdef)
end

struct CRef{name, T}
    CRef{name, T}() where {name, T} = new{name::Symbol, T}()
end

Base.unsafe_convert(::Type{Ptr{T}}, ::CRef{name, T}) where {name,T} = getglobal(Val(name), T)
Base.unsafe_convert(::Type{Ptr}, ::CRef{name, T}) where {name,T} = getglobal(Val(name), T)

Base.getindex(::CRef{name,T}) where {name,T} = loadglobal(Val(name), T)
Base.setindex!(::CRef{name,T}, x) where {name,T} = storeglobal!(Val(name), T, x)

function loadglobal(::Val{name}, ::Type{T}, i::Integer = 1) where {name,T}
    Base.unsafe_load(getglobal(Val(name), T))
end

getglobal(::Val{name}, ::Type{T}) where {name,T} = cglobal((name, libmujoco), T)

function storeglobal!(::Val{name}, ::Type{T}, x) where {name,T}
    ptr = getglobal(Val(name), T)
    Base.unsafe_store!(ptr, Base.unsafe_convert(T, Base.cconvert(T, x)))
end