text:
module CoreEngine

import DataFrames
import JSON

import Dance.Configuration
import Dance.CoreRenderer
import Dance.Logger
import Dance.Router
import Dance.PayloadUtils
import Dance.URIUtils
include("FlamencoEngine.jl")

const OUTPUT_DATA_FORMATS = Union{DataFrames.DataFrame, Dict, String}
const ROUTE_PARAMS_FORMATS = Union{Float64, Int64, String}

"""
    build_route_params_dict(;request_route_segments::Array{String, 1}, route_path::String) :: Dict

Build dict of route params and corresponding values, by:
- Parsing route.path for only segment sections
- Looking-up request route segments at corresponding index
"""
function build_route_params_dict(;request_route_segments::Array{String, 1}, route_path::String) :: Dict
    route_params_dict::Dict{Symbol, ROUTE_PARAMS_FORMATS} = Dict()

    for (idx, item) in enumerate(collect(eachmatch(URIUtils.ROUTE_REGEX_PARAM_ONLY, route_path)))
        index::String = lstrip(split(item.match, ">")[1], '<')
        value::ROUTE_PARAMS_FORMATS = request_route_segments[idx]
        if tryparse(Float64, value) isa Number
            if tryparse(Int64, value) !== nothing
                value = parse(Int64, value)
            else
                value = parse(Float64, value)
            end
        end
        route_params_dict[Symbol(index)] = value
    end

    return route_params_dict
end

"""
    map_route_function_output()

Output from route.action can contain optional HTTP Headers params dict
"""
function map_route_function_output(output::Union{Tuple, OUTPUT_DATA_FORMATS})
    data::OUTPUT_DATA_FORMATS = ""
    headers::Dict{String, String} = Dict()

    if isa(output, Tuple)
        data = output[1]
        headers = output[2]
    else
        data = output
    end

    return headers, data
end

"""
    process_backend_function(;route::Router.Route, route_segments::Array{String, 1}, payload::String, headers::Dict{String, String}) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

2 main cases

- JSON POST request with payload
- GET request

If error during process, render 500 response

Render 400 if badly supplied payload data
"""
function process_backend_function(;route::Router.Route, route_segments::Array{String, 1}, payload::String, headers::Dict{String, String}) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    data::OUTPUT_DATA_FORMATS = ""
    output_headers::Dict{String, String} = Dict()
    output::Union{Tuple, OUTPUT_DATA_FORMATS} = ""
    received_data::Union{DataFrames.DataFrame, Union{Array{Any,1}, Dict}} = Dict()
    rendered_dict::Dict{Symbol, Union{Dict{String, String}, Int16, String}} = Dict()
    route_params_dict::Dict{Symbol, ROUTE_PARAMS_FORMATS} = Dict()

    if route.endpoint==Router.EP_JSON && length(payload)>0
        can_proceed::Bool = false
        try
            json_decoded_data::Union{Array{Any,1}, Dict} = JSON.parse(payload)
            if isa(json_decoded_data, Array{Any,1})
                received_data = PayloadUtils.convert_array_to_dataframe(json_decoded_data)
            else
                received_data = json_decoded_data
            end
            can_proceed = true
        catch e
            rendered_dict = render_400(;endpoint=route.endpoint)
        end
        if can_proceed
            try
                if route.has_regex
                    route_params_dict = build_route_params_dict(;request_route_segments=route_segments, route_path=route.path)
                    output = route.action(route_params_dict, received_data, headers)
                else
                    output = route.action(received_data, headers)
                end
                output_headers, data = map_route_function_output(output)
                rendered_dict = render_200(;headers=output_headers, endpoint=route.endpoint, data=data, html_file=route.html_file)
            catch e
                rendered_dict = render_500(;endpoint=route.endpoint, data=e, html_file=route.html_file, request_path=route.path)
            end
        end
    elseif route.endpoint==Router.EP_JSON || route.endpoint==Router.EP_HTML
        try
            if route.has_regex
                route_params_dict = build_route_params_dict(;request_route_segments=route_segments, route_path=route.path)
                output = route.action(route_params_dict, headers)
            else
                output = route.action(headers)
            end
            output_headers, data = map_route_function_output(output)
            rendered_dict = render_200(;headers=output_headers, endpoint=route.endpoint, data=data, html_file=route.html_file)
        catch e
            rendered_dict = render_500(;endpoint=route.endpoint, data=e, html_file=route.html_file, request_path=route.path)
        end
    # Static case (no headers param passed)
    else
        try
            output = route.action(route.path)
            output_headers, data = map_route_function_output(output)
            rendered_dict = render_200(;headers=output_headers, endpoint=route.endpoint, data=data, html_file=route.html_file)
        catch e
            rendered_dict = render_500(;endpoint=route.endpoint, data=e, html_file=route.html_file, request_path=route.path)
        end
    end

    return rendered_dict
end

"""
    render_200(;headers::Dict{String, String}, endpoint::String, data::OUTPUT_DATA_FORMATS, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Render HTTP 200 response
"""
function render_200(;headers::Dict{String, String}, endpoint::String, data::OUTPUT_DATA_FORMATS, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    return CoreRenderer.render(;headers=headers, status_code=Int16(200), endpoint=endpoint, data=data, html_file=html_file)
end

"""
    render_400(;endpoint::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Render HTTP 400 response

No HTML output file required here, as limited to JSON request
"""
function render_400(;endpoint::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    return CoreRenderer.render(;status_code=Int16(400), endpoint=endpoint, data="Bad Request", html_file="")
end

"""
    render_404(;endpoint::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Render HTTP 404 response
"""
function render_404(;endpoint::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    return CoreRenderer.render(;status_code=Int16(404), endpoint=endpoint, data="Not Found", html_file=Configuration.Settings[:html_base_filename]*".html")
end

"""
    render_405(;endpoint::String, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Render HTTP 405 response
"""
function render_405(;endpoint::String, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    return CoreRenderer.render(;status_code=Int16(405), endpoint=endpoint, data="Method Not Allowed", html_file=html_file)
end

"""
    render_500(;endpoint::String, data::String, html_file::String, request_path::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Render HTTP 500 response

Only if NOT Prod: output error details to page
"""
function render_500(;endpoint::String, data::String, html_file::String, request_path::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    Logger.log("500 Internal Server Error when rendering url: $request_path")

    ## Only output error details if NOT Prod ##
    if Configuration.is_prod()
        data = "Internal Server Error"
    else
        data = "Internal Server Error: " * data
    end

    return CoreRenderer.render(;status_code=Int16(500), endpoint=endpoint, data=data, html_file=html_file)
end

"""
    render(;request_headers::Array, request_method::String, request_path::String, request_payload::String) :: String

Main entry point of HTTP Rendering

- Invalid method will render 405 response
- If error during process, render 500 response

Check Headers if was JSON request, to return 404 as JSON output
"""
function render(;request_headers::Dict{String, String}, request_method::String, request_path::String, request_payload::String) ::String
    rendered_dict::Dict{Symbol, Union{Dict{String, String}, Int16, String}} = Dict()
    request_route_segments::Array{String, 1} = []

    function _render_404_from_content_type(request_headers::Dict{String, String}) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
        ## Check Headers if was JSON request, to return 404 in JSON format ##
        content_type::String = ""
        endpoint::String = Router.EP_HTML

        try
            if haskey(request_headers, "Content-Type")
                content_type = request_headers["Content-Type"]
            end
            if content_type=="application/json"
                endpoint = Router.EP_JSON
            end
        catch e
            nothing
        end

        return render_404(;endpoint=endpoint)
    end

    route::Union{Router.Route, Nothing} = Router.get_route(request_path)
    if isa(route, Router.Route)
        try
            ## Check if method allowed ##
            if request_method == Router.OPTIONS
                headers_dict::Dict{String, String} = Dict(
                    "Allow" => route.method,
                    "Access-Control-Allow-Methods" => route.method,
                    "Access-Control-Allow-Headers" => "X-PINGOTHER, Content-Type"
                )
                rendered_dict = render_200(;headers=headers_dict, endpoint=route.endpoint, data=Dict(), html_file="")
            else
                if request_method==route.method
                    if route.has_regex
                        request_route_segments = URIUtils.get_path_param_segments(;request_path=request_path, route_path=route.path)
                    end
                    rendered_dict = process_backend_function(;route=route, route_segments=request_route_segments, payload=request_payload, headers=request_headers)
                else
                    rendered_dict = render_405(;endpoint=route.endpoint, html_file=route.html_file)
                end
            end
        catch e
            rendered_dict = render_500(;endpoint=route.endpoint, data=string(e), html_file=route.html_file, request_path=route.path)
        end
    else
        rendered_dict = _render_404_from_content_type(request_headers)
    end

    # Remove empty pair from Header dict (CoreRenderer default)
    if haskey(rendered_dict[:headers], "")
        delete!(rendered_dict[:headers], "")
    end

    ## Automaticlly set `Content-Type` Header ##
    rendered_dict[:headers]["Content-Type"] = rendered_dict[:content_type]

    return respond(;
        headers=rendered_dict[:headers], status_code=rendered_dict[:status_code], body=rendered_dict[:body]
    )
end

end


text:
# Practical Machine learning
# Decision Trees example
# Chapter 5

using Base.Test
using RDatasets
using DecisionTree

iris = dataset("datasets", "iris")
features = array(iris[:, 1:4]);
labels = array(iris[:, 5]);

# train full-tree classifier
model = build_tree(labels, features)
# prune tree: merge leaves having >= 90% combined purity (default: 100%)
model = prune_tree(model, 0.9)
# pretty print of the tree, to a depth of 5 nodes (optional)
print_tree(model, 5)
# apply learned model
apply_tree(model, [5.9,3.0,5.1,1.9])
# run n-fold cross validation for pruned tree, using 90% purity threshold purning, and 3 CV folds
println("\n##### nfoldCV Classification Tree #####")
accuracy = nfoldCV_tree(labels, features, 0.9, 3)
@test mean(accuracy) > 0.8

# train random forest classifier, using 2 random features, 10 trees and 0.5 of samples per tree (optional, defaults to 0.7)
model = build_forest(labels, features, 2, 10, 0.5)
# apply learned model
apply_forest(model, [5.9,3.0,5.1,1.9])
# run n-fold cross validation for forests, using 2 random features, 10 trees, 3 folds, 0.5 of samples per tree (optional, defaults to 0.7)
println("\n##### nfoldCV Classification Forest #####")
accuracy = nfoldCV_forest(labels, features, 2, 10, 3, 0.5)
@test mean(accuracy) > 0.8

# train adaptive-boosted decision stumps, using 7 iterations
model, coeffs = build_adaboost_stumps(labels, features, 7);
# apply learned model
apply_adaboost_stumps(model, coeffs, [5.9,3.0,5.1,1.9])
# run n-fold cross validation for boosted stumps, using 7 iterations and 3 folds
println("\n##### nfoldCV Classification Adaboosted Stumps #####")
nfoldCV_stumps(labels, features, 7, 3)

ext: jl


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter
using GeoGreensFunctions

makedocs(
    doctest=false,
    modules = [GeoGreensFunctions],
    sitename = "GeoGreensFunctions",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "APIs" => "APIs.md",
    ],
)

deploydocs(
  repo = "github.com/shipengcheng1230/GeoGreensFunctions.jl.git",
  target = "build",
)


text:
using ADCME

CONDA_ROOT = abspath(joinpath(ADCME.LIBDIR, ".."))
change_directory(ADCME.PREFIXDIR)
http_file("https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-4.0.4.tar.gz", "openmpi.tar.gz")
uncompress("openmpi.tar.gz", "openmpi-4.0.4")
change_directory(joinpath("openmpi-4.0.4", "build"))
require_file("Makefile") do
    run_with_env(`../configure CC=$(ADCME.CC) CXX=$(ADCME.CXX) --enable-mpi-thread-multiple --prefix=$(CONDA_ROOT)
                    --enable-mpirun-prefix-by-default --enable-mpi-fortran=no --with-mpi-param-check=always`, Dict("LDFLAGS"=>"-L"*ADCME.LIBDIR))
end
require_library("mpi") do
    run_with_env(`make -j all`)
end
require_library(joinpath(ADCME.LIBDIR, "libmpi")) do
    run_with_env(`make install`)
end


text:
mutable struct CartesianProductIt{T, U}
  ranges::Vector{T}
  inplace::Bool
  value::Vector{U}

  function CartesianProductIt{T, U}() where {T, U}
    return new{T, U}()
  end
end

function cartesian_product_iterator(x::T, n::Int; inplace::Bool = true) where T
  v = fill(x, n)
  return cartesian_product_iterator(v, inplace = inplace)
end

function cartesian_product_iterator(v::Vector{T}; inplace::Bool = true) where T
  it = CartesianProductIt{T, eltype(T)}()
  it.ranges = v
  it.inplace = inplace
  it.value = Vector{eltype(T)}(undef, length(v))
  return it
end

function Base.iterate(F::CartesianProductIt{T, U}) where {T, U}
  r = length(F.ranges)
  if r == 0
    return nothing
  end
  x = iterate(F.ranges[1])
  if x isa Nothing
    return nothing
  end
  a, b = x
  st = Vector{typeof(b)}(undef, r)
  st[1] = b
  F.value[1] = a
  for i = 2:r
    a, b = iterate(F.ranges[i])
    st[i] = b
    F.value[i] = a
  end

  if !F.inplace
    ret = deepcopy(F.value)
  else
    ret = F.value
  end

  return ret, st
end

function Base.iterate(F::CartesianProductIt, st::Vector{S}) where {S}
  i = 1
  r = length(F.ranges) + 1
  while i != r
    it = iterate(F.ranges[i], st[i])
    if it === nothing
      a, b = iterate(F.ranges[i])
      F.value[i] = a
      st[i] = b
      #st[i] = first(F.ranges[i])
      i += 1
      continue
    end
    F.value[i] = it[1]
    st[i] = it[2]

    if !F.inplace
      ret = deepcopy(F.value)
    else
      ret = F.value
    end

    return ret, st
  end
  return nothing
end

function Base.length(F::CartesianProductIt)
  return prod(length(x) for x in F.ranges)
end

Base.IteratorSize(::Type{CartesianProductIt{T}}) where T = Base.HasLength()

Base.eltype(::CartesianProductIt{T}) where T = Vector{eltype(T)}

function Base.getindex(F::CartesianProductIt{T}, i::Int) where T
  v = Vector{eltype{T}}(undef, length(F.ranges))
  ind = i
  for s = length(v):-1:1
    d = prod(Int[length(F.ranges[j]) for j = 1:s-1])
    q, ind = divrem(ind, d)
    if iszero(ind)
      v[s] = F.ranges[s][q]
      for t = s-1:-1:1
        v[t] = last(F.ranges[t])
      end
      return v
    end
    v[s] = F.ranges[s][q+1]
  end
  return v
end

function Base.getindex(F::CartesianProductIt, coordinates::Vector{Int})
  @assert length(coordinates) == length(F.ranges)
  return Int[F.ranges[i][coordinates[i]] for i = 1:length(F.ranges)]
end


text:
using ClassicalOrthogonalPolynomials, Plots

##
#  Ideal fluid flow consists of level sets of the imagainary part of a function
# that is asymptotic to c*z and whose imaginary part vanishes on Γ
#
#
# On the unit interval, -2*H*u gives the imaginary part of cauchy(u,x)
#  So if we want to find u defined on Γ so that hilbert(u,x) = imag(c*x)
#  then c*z + 2cauchy(u,z) vanishes on Γ
##

T = ChebyshevT()
U = ChebyshevU()
x = axes(U,1)
H = inv.(x .- x')

c = exp(0.5im)

u = Weighted(U) * ((H * Weighted(U)) \ imag(c * x))

ε  = eps(); (inv.(0.1+ε*im .- x') * u + inv.(0.1-ε*im .- x') * u)/2 ≈ imag(c*0.1)
ε  = eps(); real(inv.(0.1+ε*im .- x') * u ) ≈ imag(c*0.1)

v = (s,t) -> (z = (s + im*t); imag(c*z) - real(inv.(z .- x') * u))

xx = range(-3,3; length=100)
yy = range(-1,1; length=100)
plot([-1,1],[0,0]; color=:black)
contour!(xx,yy,v.(xx',yy))


text:
#
# TODO: add isEmpty test
#
module TestConfig
  using Config
  using FactCheck

  facts("Read default value") do
    @fact Config.create().worldIpsPeriod --> 1.0
  end

  facts("Save configuration into the file") do
    local cfgFile = "config.data"
    local conf    = Config.create()

    try rm(cfgFile) end
    @fact Config.save(conf, cfgFile) --> true
    @fact isfile(cfgFile)      --> true
    rm(cfgFile)
  end
  facts("Save/load configuration into/from file with correct data") do
    local cfg::Config.ConfigData = Config.create()
    local cfgFile::String = "config.data"

    try rm(cfgFile) end
    cfg.worldWidth = 666
    @fact Config.save(cfg, cfgFile)   --> true
    cfg.worldWidth = 777
    cfg = Config.load(cfgFile)
    @fact cfg.worldWidth --> 666
    rm(cfgFile)
  end
  facts("Load configuration from incorrect file") do
    local cfgFile::String = "config.data"
    local conf = Config.create()

    try rm(cfgFile) end
    io = open(cfgFile, "w")
    write(io, "Hello!")
    close(io)
    println("val: ", typeof(Config.load(cfgFile)))
    @fact Config.isEmpty(Config.load(cfgFile)) --> true
    rm(cfgFile)
  end
  facts("Load configuration from file, which doesn't exist") do
    local cfgFile::String = "config.data"
    try rm(cfgFile) end
    @fact Config.isEmpty(Config.load(cfgFile)) --> true
  end

  facts("Combined test of all functions") do
    local cfg::Config.ConfigData = Config.create()
    local cfgFile1::String = "config1.data"
    local cfgFile2::String = "config2.data"

    try rm(cfgFile1) end
    try rm(cfgFile2) end
    cfg.worldWidth = 666
    @fact Config.save(cfg, cfgFile1)       --> true # IPS === 666
    @fact cfg.worldWidth --> 666
    cfg.worldWidth = 777
    @fact cfg.worldWidth                    --> 777
    cfg = Config.load(cfgFile1)
    @fact Config.isEmpty(cfg)              --> false
    @fact cfg.worldWidth                    --> 666
    cfg.worldWidth = 777
    @fact cfg.worldWidth                    --> 777

    @fact Config.save(cfg, cfgFile2)       --> true      # IPS === 777
    cfg.worldWidth = 666
    @fact cfg.worldWidth                    --> 666
    cfg = Config.load(cfgFile2)
    @fact Config.isEmpty(cfg)              --> false
    @fact cfg.worldWidth                    --> 777

    cfg = Config.load(cfgFile1)
    @fact Config.isEmpty(cfg)              --> false     # IPS === 666
    @fact cfg.worldWidth                    --> 666
    cfg = Config.load(cfgFile2)
    @fact Config.isEmpty(cfg)              --> false     # IPS === 777
    @fact cfg.worldWidth                    --> 777

    rm(cfgFile1)
    rm(cfgFile2)
  end
 end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
export CoreGeometry
export BHloop, SpecificPowerLossData
export FerriteProperties
export PCB_Specification

include("ConverterTypes.jl")
"""
    CoreGeometry

**Fields**
- `name`                      -- string identifying the core
  `winding_aperture_height`  -- m
- `winding_aperture`          -- m
- `half_center_width`         -- m
- `center_length`             -- m
- `effective_volume`          -- m^3
- `effective_area`            -- m^2
- `effective_length`          -- m
- `mass`                      -- Kg
"""
struct CoreGeometry
  # all values are in MKS units
  name :: String
  winding_aperture_height :: Float64
  winding_aperture :: Float64
  half_center_width :: Float64  # center_radius for ER
  center_length :: Float64
  effective_volume :: Float64
  effective_area :: Float64
  effective_length :: Float64
  mass :: Float64
end

"""
    BHloop(hc,bs,br)

Define the BH loop of a magnetic material.

**Fields**
- `hc`      -- Coercive force (A-turn/m)
- `bs`      -- Remnant flux density (Tesla)
- `br`      -- Saturation flux density (Tesla)
"""
struct BHloop
  hc :: Float64
  bs :: Float64
  br :: Float64
end
hc(x::BHloop) = x.hc
bs(x::BHloop) = x.bs
br(x::BHloop) = x.br

"""
    SpecificPowerLossData(frequency::Tuple, mb::Tuple)
    SpecificPowerLossData(input::SplInput)

Capture specific power loss data from datasheet.

Data is stored as a series of linear approximations on a log log plot, one for
each frequency.

**Fields**
- `frequency`   -- frequency of linear approximation (Hz)
- `mb`          -- Tuple (slope, offset) defining the linear approximation

In order to simplify manual data entry, data my also be passed as a SplInput
  object.  All data is in MKS units (Hz, Tesla, W/m^3).
"""
SpecificPowerLossData

struct SpecificPowerLossData
  frequency :: Tuple  # in Hz
  mb :: Tuple  # y = mx + b
  # y in W/m^3
  # x in Tesla
  function SpecificPowerLossData(frequency,mb)
    if length(frequency) == 0
      throw(ArgumentError("no frequency list"))
    end
    if length(frequency) == 1
      throw(ArgumentError("single frequency not allowed"))
    end
    if length(frequency)!=length(mb)
      throw(ArgumentError("lengths do not match"))
    end
    new(frequency,mb)
  end
end

"""
    FerriteProperties(frequency_range, troom, thot, bh_room, bh_hot, spl_hot)

Store material data.

**Fields**
- `name`              -- string identifying the ferrite
- `fmin`              -- minimum recommended operating frequency
- `fmax`              -- maximum recommended operating frequency
- `troom`             -- typicaly 25C
- `thot`              -- typicaly 100C
- `bh_room`           -- BH loop at room temperature
- `bh_hot`            -- BH loop at hot temperature
- `spl_hot`           -- specfic power loss data at hot temperature
"""
struct FerriteProperties
  name :: String
  fmin :: Float64
  fmax :: Float64
  troom :: Float64
  thot :: Float64
  bh_room :: BHloop
  bh_hot :: BHloop
  spl_hot :: SpecificPowerLossData
end

include("Stackup.jl")

"""
    PCB_Specification(trace_edge_gap  :: Float64,
                      trace_trace_gap :: Float64,
                      stackup         :: Stackup)

Store PCB data.

**Fields**
- `trace_edge_gap`          -- minimum distance from trace to PCB edge (m)
- `trace_trace_gap`         -- minimum distance between traces (m)
- `stackup`                 -- defines material and thickness of the layers
"""
struct PCB_Specification
  trace_edge_gap :: Float64
  trace_trace_gap :: Float64
  stackup :: Stackup
end

include("Winding.jl")
include("Transformer.jl")

struct ChanInductor
  hc :: Float64
  bs :: Float64
  br :: Float64
  a :: Float64
  lm :: Float64
  lg :: Float64
  n :: Float64
end
function Base.show(io::IO,ci::ChanInductor)
  println(io,"Hc=",ci.hc,", Bs=",ci.bs,", Br=",ci.br,", A=",
          ci.a,", Lm=",ci.lm,", Lg=",ci.lg,", N=",ci.n)
end


text:
@testset "Test Distributed Order Relaxation" begin
    h = 0.5; t = collect(0:h:1)
    fun(t)=0
    prob = DODEProblem([1, 0.1], [x->6*x*(1-x), 0], (0, 1), fun, 1, t)
    result = solve(prob, h, DOMatrixDiscrete())
    @test isapprox(result.u, [0.934688468068206, 0.9020578767132177, 0.879670766076252]; atol=1e-3)
end

@testset "Test Distributed Order Oscillator" begin
    h = 0.5; t = collect(0:h:1)
    dooscillator(t) = 0 ≤ t ≤ 1 ? (return 8) : (return 0)
    prob = DODEProblem([1, 1, 1], [2, x->6*x*(1-x), 0], (0, 1), dooscillator, [0; 0], t)
    sol = solve(prob, h, DOMatrixDiscrete())
    @test isapprox(sol.u, [0.0, 0.0, 1.243950673139022]; atol=1e-3)
end


text:
test = """
    6,10
    0,14
    9,10
    0,3
    10,4
    4,11
    6,0
    6,12
    4,1
    0,13
    10,12
    3,4
    3,0
    8,4
    1,10
    2,14
    8,10
    9,0

    fold along y=7
    fold along x=5
    """

function parse_input(input::AbstractString)
    dots, folds = split.(split(input, "\n\n"), r"\n"; keepempty = false)
    dots = map(p -> parse.(Int, split(p, ',')), dots)
    folds = [let m = match(r"^fold along (.)=([0-9]+)$", fold)
                 (Symbol(m.captures[1]), parse(Int, m.captures[2]))
             end for fold in folds]
    card = fill(false, maximum(last.(dots)) + 1, maximum(first.(dots)) + 1)
    for (x, y) in dots
        card[y + 1, x + 1] = true
    end
    return (card, folds)
end

fold(card::AbstractMatrix, (axis, v)) = fold(Val(axis), card, v)

function fold(::Val{:x}, card::AbstractMatrix, v::Int)
    return card[:, 1:v] .| card[:, 2v + 1:-1:v + 2]
end

function  fold(::Val{:y}, card::AbstractMatrix, v::Int)
    return card[1:v, :] .| card[2v + 1:-1:v + 2, :]
end

q1(card, folds) = sum(fold(card, folds[1]))
q2(card, folds) = reduce(fold, folds; init = card) |> pretty

function pretty(card::AbstractMatrix)::String
    join([join([v ? "##" : "  " for v in row], "")
          for row in eachrow(card)], "\n")
end

let (card, folds) = parse_input(test)
    @assert q1(card, folds) == 17
end

(card, folds) = parse_input(read("day13.in", String))
println("Q1: ", q1(card, folds))
println("Q2:\n", q2(card, folds))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using Plots
using YAML
using Configurations
using Distributions

include("gen_clients.jl")
include("load_data.jl")
include("save_results.jl")
include("build_instance.jl")

@option struct Config
    map_size::Int64
    line_id::String
    day_id::String
    stat_id::String
    min_var::Float64
    max_var::Float64
    scale::Int64
    data_file::String
    result_file::String
    result_key::String
    plot::Bool
    csr_mean::Float64
    csr_var::Float64
    client_mean::Float64
    client_var::Float64
end

function get_log_normal_params(μ, σ)
    x = log(μ^2 / √(μ^2 + σ^2))
    y = log(1 + σ^2 / μ^2)
    x, y
end

function log_normal_factory(μ, σ)
    x, y = get_log_normal_params(μ, σ)
    LogNormal(x, y)
end

function factory_distributions(config::Config)
    d₁ = log_normal_factory(config.csr_mean, config.csr_var)
    d₂ = log_normal_factory(config.client_mean, config.client_var)
    d₁, d₂
end

function main(config::Config)

    data = load_data(config.data_file)

    counts = get_data(
        data, config.line_id, config.day_id, config.stat_id)
    clients = generate_clients(
        counts, config.map_size, config.min_var,
        config.max_var, config.scale)

    close_data(data)

    d₁, d₂ = factory_distributions(config)
    instance = build_instance(clients, d₁, d₂)
    save_results(config.result_file, config.result_key, instance)
    if config.plot
        x, y = size(clients)
        contour(1:x, 1:y, clients)
    end
end

conf_dict = YAML.load_file("config.yml"; dicttype=Dict{String, Any})
config = from_dict(Config, conf_dict)
main(config)


text:
"""
    XCFunctional

An enumerated type corresponding to different kinds of exchange-correlation functionals. The
enumeration corresponds to the values used by abinit.

The enumeration of functionals is identical to the one given in the abinit documentation:
https://docs.abinit.org/variables/basic/#ixc

Functionals are described as `LDA` for local density approximations (without support for spin),
`LSDA` for local spin-density approximations (which support systems with and without spin), and
`GGA` for generalized gradient approximations.
"""
@enum XCFunctional begin
    NOXC = 0
    LSDA_TETER = 1
    LDA_PZCA = 2
    LDA_TETER_OLD = 3
    LDA_WIGNER = 4
    LDA_HEDIN_LUNDQVIST = 5
    LDA_X_ALPHA = 6
    LSDA_PERDEW_WANG = 7
    LSDA_PERDEW_WANG_X = 8
    LSDA_PERDEW_WANG_X_RPA = 9
    GGA_PBE = 11
    GGA_PBE_X = 12
    GGA_VAN_LEEUWEN_BAERENDS = 13
    GGA_REVPBE = 14
    GGA_RPBE = 15
    GGA_HTCH93 = 16
    GGA_HTCH120 = 17
    FERMI_AMALDI = 20
    FERMI_AMALDI_LDA = 21
    FERMI_AMALDI_BPG = 22
    GGA_WU2006 = 23
    GGA_C09X = 24
    GGA_HTCH147 = 26
    GGA_HTCH407 = 27
    HF = 40
    PBE0 = 41
    PBE013 = 42
end

"""
    HGHPseudopotential

Hartwigsen-Goedecker-Hutter pseudopotentials, as defined in "Relativistic separable dual-space
Gaussian pseudopotentials from H to Rn" (Hartwigsen, Goedecker, Hutter, 1990) and used in
abinit norm-conserving pseudopotential calculations.

HGH pseudopotentials are relativistic and separable, and consist of a local (real space) portion
and a nonlocal (reciprocal space) portion.

Although HGH pseudopotential data is freely available, and in principle we could load in all
existing HGH pseudopotentials into this package, this data structure allows alteration of the
pseudopotential parameters.
"""
struct HGHPseudopotential <: AbstractPotential
    # Atomic number
    zatom::Int8
    # Total charge
    zion::Int8
    # Exchange-correlation functional to be used
    pspxc::XCFunctional
    # Maximum value of l (for spherical harmonics)
    lmax::Int8
    # Not sure what this particular value is...
    rloc::Float64
    # Some coefficients
    c::SVector{4,Float64}
    # r coefficients
    r::Vector{Float64}
    # h and k coefficients: dimensions 3*3*(lmax+1)
    # Note that values at l need to be accessed using the index [i,j,l+1]
    h::Array{Float64,3}
    k::Array{Float64,3}
end

"""
    zatom(p::AbstractPotential)

Gets the atomic number associated with a potential.
"""
zatom(p::AbstractPotential) = p.zatom

"""
    zion(psp::AbstractPseudopotential)

Gets the effective charge associated with a pseudopotential.
"""
zion(psp::AbstractPseudopotential) = psp.zion
zion(p::AbstractPotential) = zatom(p)

"""
    AlchemicalPotential{T<:AbstractPotential}

An alchemical potential, which can be used to describe a crystal site with mixed occupancy.

The `pot` field consists of the potentials used to construct the alchemical potential, and the
`coeff` field is a `Vector{Float64}` describing the relative contribution of each listed potential.
"""
struct AlchemicalPotential{T<:AbstractPotential} <: AbstractPotential
    pot::Vector{T}
    coeff::Vector{Float64}
    function AlchemicalPotential(pot::AbstractVector{T}, coeff::AbstractVector{<:Real}) where T
        @assert length(pot) == length(coeff) string(
            "Potential and coefficient vectors have different lengths"
        )
        # Remove duplicate values
    end
end

max_stars_repo_name: brainandforce/Xtal

max_issues_repo_name: brainandforce/Xtal

max_forks_repo_name: brainandforce/Xtal


text:
using GenericLinearAlgebra
using LinearAlgebra
using SpecialFunctions
using Dates
using PolynomialRoots
using Printf
using FFTW
using PyPlot
using MAT

function GSLW(P, p_P, stop_criteria, R_init, R_max, if_polish = false)
    #------------------------------------------------------------------------------------------------------------
    # Input:
    #     P: Chebyshev coefficients of polynomial P, only need to provide non-zero coefficient.
    #        P should satisfy parity constraint and |P|^2 \le 1
    #     p_P: parity of P, 0 -- even, 1 -- odd
    #     stop_eps: algorithm will stop if it find factors that approximate polynomials with error less than
    #               stop_eps on Chebyshev points
    #     R_init: number of bits used at the beginning
    #     R_max: number of bits available
    #     if_polish: whether polish the roots, it cost more time but may improve performance
    #
    # Output:
    #    Phase factor Phi such that real part of (0,0)-component of U_\Phi(x) approximates P(x),
    #    where U_\Phi(x) = e^{\I \phi_0 \sigma_z} \prod_{j=1}^{\qspdeg} \left[ W(x) e^{\I \phi_j \sigma_z} \right]
    #
    #    Besides, the algorithm will return the L^∞ error of such approximation on Chebyshev points
    #
    #------------------------------------------------------------------------------------------------------------
    #
    # Reference:
    #     A. Gily ́en, Y. Su, G. H. Low, and N. Wiebe.
    #     Quantum singular value transformation and beyond: exponentialimprovements for quantum matrix arithmetics.
    #
    # Author: X.Meng, Y. Dong
    # Version 1.0 .... 02/2020
    #
    #------------------------------------------------------------------------------------------------------------

    # Step 1: Find all the roots of 1-P^2

    R = R_init
    while(true)
        if(R>=R_max)
            return Inf,[]
        end
        setprecision(BigFloat, R)
        P = big.(P)
        degree = length(P)

        if(p_P==0)
            coef_PQ = zeros(BigFloat,degree*2)
            for j=1:degree
                for k=1:degree
                    coef_PQ[j+k-1] -= P[j]*P[k]
                end
            end
            coef_PQ[1] += 1
            coef_PQ1 = zeros(BigFloat,degree)
            for j=1:degree
                coef_PQ1[j] -= P[j]
            end
            coef_PQ1[1] += 1
            coef_PQ2 = zeros(BigFloat,degree)
            for j=1:degree
                coef_PQ2[j] += P[j]
            end
            coef_PQ2[1] += 1
            Proot1 = roots(coef_PQ1, polish = if_polish, epsilon = big.(0.0))
            Proot2 = roots(coef_PQ2, polish = if_polish, epsilon = big.(0.0))
            Proot = [Proot1;Proot2]
        else
            coef_PQ = zeros(BigFloat,degree*2)
            for j=1:degree
                for k=1:degree
                    coef_PQ[j+k] -= P[j]*P[k]
                end
            end
            coef_PQ[1] += 1
            Proot = roots(coef_PQ, polish = if_polish, epsilon = big.(0.0))
        end

    # Step 2: Find root of 1-P^2, construct full P and Q by FFT

        # recover full root list

        all_root = zeros(Complex{BigFloat},length(Proot)*2)
        for i=1:length(Proot)
            tmpnorm = norm(Proot[i])
            tmpangle = angle(Proot[i])
            all_root[2*i-1] = sqrt(tmpnorm)*exp(1im*tmpangle/2)
            all_root[2*i] = -sqrt(tmpnorm)*exp(1im*tmpangle/2)
        end

        # Construct W such that W(x)W(x)^*=1-P^2(x)

        eps = 1e-16
        S_0 = 0
        S_1 = 0
        S_2 = 0
        S_3 = 0
        S_4 = 0
        S1_list = zeros(Complex{BigFloat},length(all_root))
        S2_list = zeros(Complex{BigFloat},length(all_root))
        S3_list = zeros(Complex{BigFloat},length(all_root))
        S4_list = zeros(Complex{BigFloat},length(all_root))
        for i=1:length(all_root)
            if(abs(all_root[i])<eps)
                S_0 += 1
                continue
            end
            if(abs(imag(all_root[i]))<eps&&real(all_root[i])>0)
                if(real(all_root[i])<1-eps)
                    S_1 += 1
                    S1_list[S_1] = real(all_root[i])
                else
                    S_2 += 1
                    S2_list[S_2] = findmax([real(all_root[i]),1])[1]
                end
                continue
            end
            if(abs(real(all_root[i]))<eps&&imag(all_root[i])>0)
                S_3 += 1
                S3_list[S_3] = all_root[i]
                continue
            end
            if(imag(all_root[i])>0&&real(all_root[i])>0)
                S_4 += 1
                S4_list[S_4] = all_root[i]
            end
        end
        K = abs(P[end])

        function get_w(x,use_real = true) # W(x)
            x = big.(x)
            Wx = K*x^(S_0/2)
            eps3 = 1e-24
            if(x==1) # if x==\pm 1, silghtly move x such that make \sqrt{1-x^2}>0
                x -= eps3
            elseif(x==-1)
                x += eps3
            end
            for i=1:S_1
                Wx *= sqrt(x^2-S1_list[i]^2)
            end
            for i=1:S_2
                Wx *= sqrt(S2_list[i]^2-big.(1))*x+im*S2_list[i]*sqrt(big.(1)-x^2)
            end
            for i=1:S_3
                Wx *= sqrt(abs(S3_list[i])^2+big.(1))*x+im*abs(S3_list[i])*sqrt(big.(1)-x^2)
            end
            for i=1:S_4
                tmpre = real(S4_list[i])
                tmpim = imag(S4_list[i])
                tmpc = tmpre^2+tmpim^2+sqrt(2*(tmpre^2+1)*tmpim^2+(tmpre^2-1)^2+tmpim^4)
                Wx *= tmpc*x^2-(tmpre^2+tmpim^2)+im*sqrt(tmpc^2-1)*x*sqrt(big.(1)-x^2)
            end
            if(use_real)
                return real(Wx)
            else
                return imag(Wx)/sqrt(big.(1)-x^2)
            end
        end

        function get_p(x) # P(x)
            P_t = big.(0)
            for j=1:length(P)
                if(p_P==1)
                    P_t += P[j]*x^big.(2*j-1)
                else
                    P_t += P[j]*x^big.(2*j-2)
                end
            end
            return P_t
        end

        # Represent full P and Q under Chevyshev basis

        get_wr(x) = get_w(x,true)
        get_wi(x) = get_w(x,false)
        DEG = 2^ceil(Int,log2(degree)+1)
        coef_r = ChebyExpand(get_wr, DEG)
        coef_i = ChebyExpand(get_wi, DEG)
        coef_p = ChebyExpand(get_p, DEG)
        if(p_P==0)
            P_new = 1im.*coef_r[1:2*degree-1]+coef_p[1:2*degree-1]
            Q_new = coef_i[1:2*degree-2].*1im
        else
            P_new = 1im.*coef_r[1:2*degree]+coef_p[1:2*degree]
            Q_new = coef_i[1:2*degree-1].*1im
        end

    # Step 3: Get phase factors and check convergence

        phi = get_factor(P_new,Q_new)
        max_err = 0
        t = cos.(collect(1:2:(2*degree-1))*big.(pi)/big.(4)/big.(degree))
        for i=1:length(t)
            targ, ret = QSPGetUnitary(phi, t[i])
            P_t = big.(0)
            for j=1:degree
                if(p_P==1)
                    P_t += P[j]*t[i]^big.(2*j-1)
                else
                    P_t += P[j]*t[i]^big.(2*j-2)
                end
            end
            t_err = norm(real(ret[1,1])-P_t)
            if(t_err>max_err)
                max_err = t_err
            end
        end
        @printf("For degree N = %d, precision R = %d, the estimated inf norm of err is %5.4e\n",length(phi)-1,R,max_err)
        if(max_err<stop_criteria)
            return max_err,phi
        else
            @printf("Error is too big, increase R.\n")
        end
        R = R*2
    end
end

function get_factor(P,Q)
    # From polynomials P, Q generate phase factors phi such that
    # U_\Phi(x) = [P & i\sqrt{1-x^2}Q \\ i\sqrt{1-x^2}Q^* & P]
    # phase factors are generated via a reduction procedure under Chebyshev basis
    phi = zeros(BigFloat,length(P))
    lenP = length(P)
    for i=1:lenP-1
        P, Q, phit = ReducePQ(P, Q)
        phi[end+1-i] = real(phit)
    end
    phi[1] = angle(P[1])
    return phi
end

function ReducePQ(P, Q)
    # A single reduction step
    P = big.(P)
    Q = big.(Q)
    colP = length(P)
    colQ = length(Q)
    degQ = colQ-1

    tmp1 = zeros(Complex{BigFloat},colP+1)
    tmp2 = zeros(Complex{BigFloat},colP+1)
    tmp1[2:end] = big.(0.5)*P
    tmp2[1:end-2] = big.(0.5)*P[2:end]
    Px = tmp1 + tmp2
    Px[2] = Px[2] + big.(0.5)*P[1]

    if(degQ>0)
        tmp1 = zeros(Complex{BigFloat},colQ+2)
        tmp2 = zeros(Complex{BigFloat},colQ+2)
        tmp3 = zeros(Complex{BigFloat},colQ+2)
        tmp1[1:end-2] = big.(0.5)*Q
        tmp2[3:end] = -big.(1)/big.(4)*Q
        tmp3[1:end-4] = -big.(1)/big.(4)*Q[3:end]
        Q1_x2 = tmp1 + tmp2 + tmp3
        Q1_x2[2] = Q1_x2[2] - 1/big.(4)*Q[2]
        Q1_x2[3] = Q1_x2[3] - 1/big.(4)*Q[1]
    else
        Q1_x2 = zeros(Complex{BigFloat},3)
        Q1_x2[1] = big.(0.5)*Q[1]
        Q1_x2[end] = -big.(0.5)*Q[1]
    end

    tmp1 = zeros(Complex{BigFloat},colQ+1)
    tmp2 = zeros(Complex{BigFloat},colQ+1)
    tmp1[2:end] = big.(0.5)*Q
    tmp2[1:end-2] = big.(0.5)*Q[2:end]
    Qx = tmp1 + tmp2
    Qx[2] = Qx[2] + big.(0.5)*Q[1];

    if(degQ>0)
        ratio = P[end]/Q[end]*big.(2)
    else
        ratio = P[end]/Q[end]
    end
    phi = big.(0.5)*angle(ratio)
    rexp = exp(-1im*phi)
    Ptilde = rexp * (Px + ratio*Q1_x2)
    Qtilde = rexp * (ratio*Qx - P)
    Ptilde = Ptilde[1:degQ+1]
    Qtilde = Qtilde[1:degQ]

    return Ptilde,Qtilde,phi
end

function QSPGetUnitary(phase, x)
    # Given phase factors Phi and x, yield U_\Phi(x)
    phase = big.(phase)
    Wx = zeros(Complex{BigFloat},2,2)
    Wx[1,1] = x
    Wx[2,2] = x
    Wx[1,2] = sqrt(1-x^2)*1im
    Wx[2,1] = sqrt(1-x^2)*1im
    expphi = exp.(1im*phase)
    ret = zeros(Complex{BigFloat},2,2)
    ret[1,1] = expphi[1]
    ret[2,2] = conj(expphi[1])

    for k = 2:length(expphi)
        temp = zeros(Complex{BigFloat},2,2)
        temp[1,1] = expphi[k]
        temp[2,2] = conj(expphi[k])
        ret = ret * Wx * temp
    end
    targ = real(ret[1,1])
    return targ,ret
end

function ChebyExpand(func, maxorder)
    # Evaluate Chebyshev coefficients of a polynomial of degree at most maxorder
    M = maxorder
    theta = zeros(BigFloat,2*M)
    for i=1:2*M
        theta[i] = (i-1)*big.(pi)/M
    end
    f = func.(-cos.(theta))
    c = real.(BigFloatFFT(f))
    c = copy(c[1:M+1])
    c[2:end-1] = c[2:end-1]*2
    c[2:2:end] = -copy(c[2:2:end])
    c = c./(big.(2)*big.(M))
    return c
end

function Chebytonormal(coef)
    #Convert Chebyshev basis to normal basis
    coef = big.(coef)
    coef2 = zeros(BigFloat,length(coef))
    A = zeros(BigFloat,length(coef),length(coef))
    b = zeros(BigFloat,length(coef))
    t = cos.(collect(1:2:(2*length(coef)-1))*big.(pi)/big.(4)/big.(length(coef)))
    t2 = collect(1:2:(2*length(coef)-1))*big.(pi)/big.(4)/big.(length(coef))
    for i=1:length(coef)
        for j=1:length(coef)
            A[i,j] = t[i]^(j-1)
            b[i] += coef[j]*cos((j-1)*t2[i])
        end
    end
    coef2 = A\b
    #@printf("Error is %5.4e\n",norm(A*coef2-b))
    return coef2
end

function BigFloatFFT(x)
    # Perform FFT on vector x
    # This function only works for length(x) = 2^k
    N = length(x);
    xp = x[1:2:end];
    xpp = x[2:2:end];
    if(N>=8)
        Xp = BigFloatFFT(xp);
        Xpp = BigFloatFFT(xpp);
        X = zeros(Complex{BigFloat},N,1);
        Wn = exp.(big.(-2)im*big.(pi)*(big.(0:N/2-1))/big.(N));
        tmp = Wn .* Xpp;
        X = [(Xp + tmp);(Xp -tmp)];
    elseif(N==2)
        X = big.([1 1;1 -1])*x;
    elseif(N==4)
        X = big.([1 0 1 0; 0 1 0 -1im; 1 0 -1 0;0 1 0 1im]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1])*x;
    end
    return X
end

# Test case 1: Hamiltonian simulation
#
# Here we want to approxiamte e^{-i\tau x} by Jacobi-Anger expansion:
#
# e^{-i\tau x} = J_0(\tau)+2\sum_{k even} (-1)^{k/2}J_{k}(\tau)T_k(x)+2i\sum_{k odd} (-1)^{(k-1)/2}J_{k}(\tau) T_k(x)
#
# We truncate the series up to N = 1.4\tau+log(10^{14}), which gives an polynomial approximation of e^{-i\tau x} with
# accuracy 10^{-14}. Besides, we deal with real and imaginary part of the truncated series seperatly and divide them
# by a constant factor 2 to enhance stability.
#
# parameters
#     stop_eps: desired accuracy
#     tau: the duration \tau in Hamiltonian simulation
#     R_init: number of bits used at the beginning
#     R_max: number of bits available

stop_eps = 1e-12
tau = 100
R_init = 1024
R_max = 1025

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

phi1 = []
phi2 = []
for p_P=0:1
    N = ceil.(Int,tau*1.4+log(1e14))
    if(p_P==0)
        setprecision(BigFloat,4096)
        if(mod(N,2)==1)
            N -= 1
        end
        coef = zeros(BigFloat,N+1)
        for i=1:(round(Int,N/2)+1)
            coef[2*i-1] = (-1)^(i-1)*besselj(big.(2.0*(i-1)),tau)
        end
        coef[1] = coef[1]/2
        P = Chebytonormal(coef)
        P = P[1:2:end]
    else
        setprecision(BigFloat,4096)
        if(mod(N,2)==0)
            N += 1
        end
        coef = zeros(BigFloat,N+1)
        for i=1:round(Int,(N+1)/2)
            coef[2*i] = (-1)^(i-1)*besselj(big.(2*i-1),tau)
        end
        P = Chebytonormal(coef)[2:2:end]
    end

    start_time = time()
    err,phi = GSLW(P,p_P,stop_eps,R_init,R_max)
    elpased_time = time()-start_time
    @printf("Elapsed time is %4.2e s\n", elpased_time)
end

# Test case 2: Eigenstate filter
#
# Here we want to generate factors for the eigenstate filter function:
#
# f_n(x,\delta)=\frac{T_n(-1+2\frac{x^2-\delta^2}{1-\delta^2})}{T_n(-1+2\frac{-\delta^2}{1-\delta^2})}.
#
# We divide f_n by a constant factor \sqrt{2} to enhance stability.
#
# Reference: Lin Lin and Yu Tong
#            Solving quantum linear system problem with near-optimal complexity
#
# parameters
#     stop_eps: desired accuracy
#     n, \delta: parameters of f_n
#     R_init: number of bits used at the beginning
#     R_max: number of bits available
#

stop_eps = 1e-12
n = 100
delta = 0.03
R_init = 1024
R_max = 1025

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

function f_n(x,n,delta)
    val = copy(x)
    delta = big.(delta)
    fact = chebyshev(-big.(1)-big.(2)*delta^2/(big.(1)-delta^2),n)
    if(length(x)==1)
        return chebyshev(-big.(1)+big.(2)*(x^2-delta^2)/(big.(1)-delta^2),n)/fact
    else
        for i=1:length(x)
            val[i] = chebyshev(-1+2*(x[i]^2-delta^2)/(1-delta^2),n)/fact
        end
        return val
    end
end

function chebyshev(x,n) # T_n(x)
    if(abs(x)<=1)
        return cos(big.(n)*acos(x))
    elseif(x>1)
        return cosh(big.(n)*acosh(x))
    else
        return big.((-1)^n)*cosh(big.(n)*acosh(-x))
    end
end

# Obtain expansion of f_n under Chebyshev basis via FFT

setprecision(BigFloat,1024)
M = 2*n
theta = range(0, stop=2*pi, length=2*M+1)
theta = theta[1:2*M]
f = f_n(-cos.(theta),n,delta)
c = real(fft(f))
c = c[1:M+1]
c[2:end-1] = c[2:end-1]*2
c[2:2:end] = - c[2:2:end]
c = c / (2*M)

setprecision(BigFloat,4096)
P = Chebytonormal(c)[1:2:end]/sqrt(big.(2.0))
p_P = 0
start_time = time()
err,phi = GSLW(P,p_P,stop_eps,R_init,R_max)
elpased_time = time()-start_time
@printf("Elapsed time is %4.2e s\n", elpased_time)

# Test case 3: Matrix inversion
#
# We would like to approximate 1/x over [1/kappa,1] by a polynomial, such polynomial is generated
# by Remez algorithm and the approximation error is bounded by 10^{-6}
#
# parameters
#     stop_eps: desired accuracy
#     kappa: parameters of polynomial approximation
#     R_init: number of bits used at the beginning
#     R_max: number of bits available
#

stop_eps = 1e-12
kappa = 20
R_init = 2048
R_max = 2049

#------------------------------------------------------------------
# even approximation

# enter your path here
matpath2 = "../QSPPACK/Data/inversex/"

vars = matread(matpath2 * "coef_xeven_" * string(kappa)*"_6"* ".mat")
coef = vars["coef"]
setprecision(BigFloat,4096)
coef2 = zeros(2*length(coef)-1)
coef2[1:2:end] = coef
P = Chebytonormal(coef2)[1:2:end]
p_P = 0
start_time = time()
err,phi1 = GSLW(P,p_P,stop_eps,R_init,R_max)
elpased_time = time()-start_time
@printf("Elapsed time is %4.2e s\n", elpased_time)

# odd approximation

vars = matread(matpath2 * "coef_xodd_" * string(kappa)*"_6"* ".mat")
coef = vars["coef"]
setprecision(BigFloat,4096)
coef2 = zeros(2*length(coef))
coef2[2:2:end] = coef
P = Chebytonormal(coef2)[2:2:end]
start_time = time()
p_P = 1
err,phi = GSLW(P,p_P,stop_eps,R_init,R_max)
elpased_time = time()-start_time
@printf("Elapsed time is %4.2e s\n", elpased_time)

max_stars_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_issues_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_forks_repo_name: Aaron-Gregory/RE-QSVT-HHL


text:
@testset "gradient TV." begin
    M = Circle()
    N = Power(M,(4,))
    x = PowPoint(S1Point.([0.1,0.2,0.3,0.5]))
    tvTestξ = PowTVector( S1TVector.([-1.0,0.,0.,1.]) )
    @test gradTV(N,x) == tvTestξ
    @test gradTV(M, (x[1],x[1])) == (zeroTVector(M,x[1]), zeroTVector(M,x[1]))
    @test norm(N,x, gradTV(N,x,2) - tvTestξ) ≈ 0
    tv2Testξ = PowTVector(S1TVector.([0.,.5,-1.,0.5]))
    @test gradTV2(N,x) == tv2Testξ
    @test norm(N,x, forwardLogs(N,x)- PowTVector(S1TVector.([0.1, 0.1, 0.2, 0.]))) ≈ 0 atol=10^(-16)
    x2 = PowPoint(S1Point.([0.1,0.2,0.3]))
    N2 = Power(M,size(getValue(x2)))
    @test gradTV2(N2,x2) == PowTVector(S1TVector.(zeros(3)))
    @test gradTV2(N2,x2,2) == PowTVector(S1TVector.(zeros(3)))
    @test norm(N, x,
        gradIntrICTV12(N,x,x,x,1.,1.)[1]
        - PowTVector([S1TVector(-1.),S1TVector(0.), S1TVector(0.), S1TVector(1.)])#
        ) ≈ 0
    @test norm(N, x, gradIntrICTV12(N,x,x,x,1.,1.)[2] - zeroTVector(N,x)) ≈ 0
    @test gradTV(M, (S1Point(0.), S1Point(0.)),2) == (S1TVector(0.), S1TVector(0.))
    # 2d forward forwardLogs
    N3 = Power(M,(2,2))
    N3C = Power(M, (2,2,2) )
    x3 = PowPoint(S1Point.( [0.1 0.2;0.3 0.5] ))
    x3C = PowPoint( cat(getValue(x3), getValue(x3); dims=3) )
    tC = PowTVector(S1TVector.(cat([.2 .3; 0. 0.], [.1 0.; .2 0.] ; dims=3)))
    @test norm(N3C, x3C, forwardLogs(N3,x3)-tC) ≈ 0 atol=10^(-16)
end


text:
"""
    sampleplane(vp, v1, v2, lengtht, sizet)

Create a tuple of points on a plane and their normals based on a point (`vp`) and
two vectors (`v1`, `v2`).

# Arguments:
- `vp`: a point of the plane.
- `v1`, `v2`: two vectors that define a plane.
- `lengtht`: tuple containing the length of the plane along each side.
- `sizet`: tuple containing the number of samples along each side.
"""
function sampleplane(vp, v1, v2, lengtht, sizet)
    @assert size(v1) == size(v2) == size(vp) "Vector's size must be the same!"
    v1n = normalize(v1)
    v2n = normalize(v2)
    n = cross(v1n, v2n)
    @assert norm(n) != 0 "Vectors can't be collinear!"
    s1, s2 = sizet
    @assert s1 > 0 && s2 > 0 "Should sample more than 0."
    s1l, s2l = lengtht
    pn = [n for i in 1:s1*s2]
    vpn = vp - (s1l*v1n + s2l*v2n)/2
    ps = [vpn + v1n*u/s1l + v2n*v/s2l for u in 1:s1 for v in 1:s2]
    return (ps, pn)
end

function sampleplanefromcorner(vp, v1, v2, lengtht, sizet)
    @assert size(v1) == size(v2) == size(vp) "Vector's size must be the same!"
    v1n = normalize(v1)
    v2n = normalize(v2)
    n = cross(v1n, v2n)
    @assert norm(n) != 0 "Vectors can't be collinear!"
    s1, s2 = sizet
    @assert s1 > 0 && s2 > 0 "Should sample more than 0."
    s1l, s2l = lengtht
    pn = [n for i in 1:(s1+1)*(s2+1)]
    ps = [vp + v1n*s1l*u/s1 + v2n*s2l*v/s2 for u in 0:s1 for v in 0:s2]
    return (ps, pn)
end

"""
    sampleplane(vp, n, lengtht, sizet)

Create a tuple of points on a plane and their normals based on a point (`vp`) and
a surface normal (`n`).

# Arguments:
- `vp`: a point of the plane.
- `n`: surface normal.
- `lengtht`: tuple containing the length of the plane along each side.
- `sizet`: tuple containing the number of samples along each side.
"""
function sampleplane(vp, n, lengtht, sizet)
    arbnorm = arbitrary_orthogonal2(n)
    norm2 = cross(n, arbnorm)
    return sampleplane(vp, arbnorm, norm2, lengtht, sizet)
end

"""
    samplecylinder(ax, vp, R, h, sizet)

Create a tuple of points on a cylinder and their normals based on an axis (`ax`),
a nullpoint (`vp`), radius and height.

# Arguments:
- `ax`: axis of the cylinder.
- `vp`: nullpoint of the cylinder.
- `R`: radius of the cylinder.
- `h`: height of the cylinder.
- `sizet`: tuple containing the number of samples along the circumference and the height.
"""
function samplecylinder(ax, vp, R, h, sizet)
    @assert size(ax,1) == size(vp,1) == 3 "Axis and nullpoint must be 3 dimensional."
    @assert R > 0 && h > 0 "Can't construct a line or circle."
    s1, s2 = sizet
    @assert s1 > 1 && s2 > 1 "Should sample more than 1."

    axn = normalize(ax)
    aort = normalize(arbitrary_orthogonal2(axn))

    function dontMul(M, v, kit)
        if kit == 0
            return v
        else
            return M^(kit-1)*v
        end
    end

    rotMat = rodriguesrad(axn, 2*π/s1)
    ps = [vp + R.*normalize(dontMul(rotMat, aort, i)) + ((j-1)*h/(s2-1)).*axn   for i in 1:s1 for j in 1:s2]
    ns = [normalize(p - vp - axn*dot( axn, p-vp )) for p in ps]
    return (ps, ns)
end

"""
    samplesphere(cp, R, sizet)

Create a tuple of points on a sphere and their normals.

# Arguments:
- `cp`: centerpoint of sphere.
- `R`: radius of the sphere.
- `sizet`: tuple containing the number of samples along two great circles.
"""
function samplesphere(cp, R, sizet)
    @assert R>0 "Can't construct sphere with radius < 0."
    s1, s2 = sizet
    @assert s1 > 2 && s2 > 2 "Should sample more than 2."
    r1 = range(0, π, length = s1+2)[2:end-1]
    r2 = range(0, 2π, length = s2+1)[1:end-1]
    vert = [cp + R*SVector{3}([sin(i)*cos(j), sin(i)*sin(j), cos(i)]) for i in r1 for j in r2]
    rV = SVector(0,0,R)
    append!(vert, [cp-rV,cp+rV])

    ns = similar(vert)
    for i in eachindex(ns)
        ns[i] = normalize(vert[i]-cp)
    end
    return vert, ns
end

function samplecone(ap, ax, opangr, h, sizet)
    # s1 along ax
    # s2 along circle
    s1, s2 = sizet

    surfl = h/cos(opangr/2)
    axn = normalize(ax)
    v1s = collect(range(0, stop=surfl, length=s1+1))
    popfirst!(v1s)
    # points along the axis
    alongax = [axn*i for i in v1s]
    aort = normalize(arbitrary_orthogonal2(axn))
    ns = normalize(cross(aort, axn))
    rM = rodriguesrad(aort, opangr/2)

    # one line along the surface of the cone
    alongsurf = [rM*s for s in alongax]
    alongsurfn = fill!(similar(alongsurf), rM*ns)
    rotMat = rodriguesrad(axn, 2*π/s2)
    # translate with the apex too
    ps = [rotMat^j*s + ap for s in alongsurf for j in 1:s2]
    final_ns = [rotMat^j*s for s in alongsurfn for j in 1:s2]
    return ps, final_ns
end

"""
    normalsforplot(verts, norms, arrowsize = 0.5)

Create an array of pair of points for plotting the normals with Makie.
Only the direction of the normals is presented, their size not.

# Arguments:
- `arrowsize`: scaling factor for the size of the lines.
"""
function normalsforplot(verts, norms, arrowsize = 0.5)
    @assert size(verts) == size(norms) "They should have the same size."
    as = arrowsize
    return [verts[i] => verts[i] + as .*normalize(norms[i]) for i in 1:length(verts) ]
end

"""
    noisifyvertices(verts, allvs, scalef = 1)

Add gaussian noise to vertices.
Random subset or all vertices can be selected.

# Arguments:
- `allvs::Bool`: adds noise to every vertice if true.
- `scalef::Real`: scale the noise from the [-1,1] interval.
"""
function noisifyvertices(verts, allvs, scalef = 1)
    randis = similar(verts)
    for i in eachindex(randis)
        randis[i] = (2*scalef) .*SVector(rand(eltype(eltype(verts)),3)...) .-1
    end
    allvs && return verts+randis
    bools = rand(Bool, size(verts))
    verts[bools] += randis[bools]
    return verts
end

"""
    noisifyv_fixed(verts, allvs, scalef = 1)

Fixed version of `noisifyvertices`.

# Arguments:
- `allvs::Bool`: adds noise to every vertice if true.
- `scalef::Real`: scale the noise from the [-1,1] interval.
"""
function noisifyv_fixed(verts, allvs, scalef = 1)
    randis = similar(verts)
    for i in eachindex(randis)
        randis[i] = (2*scalef) .*eltype(verts)(rand(3)) .- scalef
    end
    allvs && return verts+randis
    bools = rand(Bool, size(verts))
    verts[bools] += randis[bools]
    return verts
end

"""
    noisifynormals(norms, maxrot)

Add gaussian noise to normals.
Rotates the normals around a random axis with `maxrot` degrees.

# Arguments:
- `maxrot::Real`: maximum rotation in degrees.
"""
function noisifynormals(norms, maxrot)
    retn = similar(norms)
    # an array of random axes
    randis = similar(norms)
    for i in eachindex(randis)
        randis[i] = SVector{3}(rand(eltype(eltype(norms)),3))
    end
    # random rotations: array of random numbers
    randrots = maxrot.*rand(eltype(eltype(norms)), length(norms)).-maxrot/2
    for i in 1:length(retn)
        nr = normalize(randis[i])
        while norm(normalize(cross(norms[i], nr))) < 0.1
            nr = SVector{3}(normalize(rand(eltype(norms),3)))
        end
        crv = cross(nr, norms[i])
        rM = rodriguesdeg(crv, randrots[i])
        retn[i] = rM*norms[i]
    end
    return retn
end

## Whole new world

"""
    applynormalnoise!(normals, normalnoise::URN=nothing)

Rotates the normals around a random axis with `maxrot` degrees.
`URN = Union{Real,Nothing}`, where nothing means, that noise is not applied.
Uses `randn()`, say normally-distributed random numbers with mean 0 and standard deviation 1.

# Arguments:
- `normals`: an array of `SVector{3}`.
- `vertexnoise::URN=nothing`: additional noise in the form of a rotation around a random axis, with `normalnoise*randn()` degrees.
"""
function applynormalnoise!(normals, normalnoise::URN=nothing)
    normalnoise === nothing && return normals
    ftype = eltype(normals[1])
    for i in eachindex(normals)
        nr = normalize(SVector{3}(randn(ftype, 3)))
        while norm(normalize(cross(normals[i], nr))) < 0.1
            nr = normalize(SVector{3}(randn(ftype, 3)))
        end
        crv = cross(nr, normals[i])
        rM = rodriguesdeg(crv, normalnoise*randn(ftype))
        normals[i] = rM*normals[i]
    end
    return normals
end

"""
    applyvertexnoise!(vertices, vertexnoise::URN=nothing)

Apply gaussian noise to all vertices.
`URN = Union{Real,Nothing}`, where nothing means, that noise is not applied.
Uses `randn()`, say normally-distributed random numbers with mean 0 and standard deviation 1.

# Arguments:
- `vertices`: an array of `SVector{3}`.
- `vertexnoise::URN=nothing`: additional noise in the form of `noise=vertexnoise*randn(3)`.
"""
function applyvertexnoise!(vertices, vertexnoise::URN=nothing)
    vertexnoise === nothing && return vertices
    @assert vertexnoise > 0 "Only makes sense if you scale by a positive number."
    # apply noise to all points
    for i in eachindex(vertices)
        vertices[i] = vertices[i] + vertexnoise*randn(3)
    end
    return vertices
end

"""
    addoutliers!(vertices, normals, indexes, outliercount::URN=nothing)

Add outliers.
`URN = Union{Real,Nothing}`, where nothing means, that noise is not applied.
Outliers are added in the 1.2*bounding box space.
Outliers are marked in the index array with -1.

# Arguments:
- `vertices`: an array of `SVector{3}`.
- `normals`: an array of `SVector{3}`.
- `indexes`: an array of indexes, that store, that given points belong to wich primitive. Outliers are marked with -1.
- `outliercount::URN=nothing`: add plus `outliercount`% outlier points-normals to the point cloud, uniformly distributed in the 1.1*bounding box space.
"""
function addoutliers!(vertices, normals, indexes, outliercount::URN=nothing)
    outliercount === nothing && return (vertices, normals, indexes)
    @assert 0 <= outliercount <= 100 "Percent should be <= 0 <= 100."
    s = size(vertices,1)
    c = clamp(0, floor(Int, outliercount * s/100), s)
    c == 0 && return (vertices, normals, indexes)

    vtype = eltype(vertices)
    ntype = eltype(normals)

    aa, bb = RANSAC.findAABB(vertices)

    extend = 0.2
    absaabb = (1+extend) .* (bb - aa)
    delta = aa - (extend/2)*(bb-aa)

    outlv = [vtype(rand(3) .* absaabb+delta) for i in 1:c]
    append!(vertices, outlv)
    outln = [normalize(ntype(2 .*rand(3) .-1)) for i in 1:c]
    append!(normals, outln)
    append!(indexes, fill(-1, c))

    return (vertices, normals, indexes)
end


text:
using Test, LowLevelJPEG

using Colors
using FixedPointNumbers
using StaticArrays

ambiguities_lljpeg = detect_ambiguities(LowLevelJPEG, Base, Core, Colors)
ambiguities_sarray = detect_ambiguities(StaticArrays, Base, Core)

@test isempty(setdiff(ambiguities_lljpeg, ambiguities_sarray))

const dir = joinpath(dirname(pathof(LowLevelJPEG)), "..", "test", "res")

function readtest(filename::AbstractString, t::Type, offset::Integer)
    open(joinpath(dir, filename), "r") do f
        seek(f, offset)
        return read(f, t)
    end
end
readtest(t::Type, offset::Integer) = readtest("julia.jpg", t, offset)

function shown_string(obj)
    io = IOBuffer()
    show(io, obj)
    String(take!(io))
end

function shown_plaintext(obj)
    io = IOBuffer()
    show(io, MIME("text/plain"), obj)
    seekstart(io)
    join(map(rstrip, eachline(io)), "\n") # remove trailing whitespace
end

function written_hex(obj)
    io = IOBuffer()
    write(io, obj)
    bytes2hex(take!(io))
end

@testset "Colors" begin
    include("colors.jl")
end
@testset "Markers" begin
    include("markers.jl")
end
@testset "CodeBits" begin
    include("codebits.jl")
end
@testset "Coeffs" begin
    include("coeffs.jl")
end
@testset "Blocks" begin
    include("blocks.jl")
end
@testset "QTs" begin
    include("qt.jl")
end
@testset "Huffman" begin
    include("huffman.jl")
end


text:
# This file assumes that the directories are organized as follows:
#
# ./SomeDirectory/RRTMGP.jl/
# ./SomeDirectory/rte-rrtmgp/
#

include("ConfigRRTMGP.jl")
include("CleanRRTMGP.jl")
traverse_input_params = false
run_clear_sky = false
run_all_sky = true

cd(joinpath(root_dir, "rte-rrtmgp", "build")) do
    run(`make`)
end

if run_clear_sky
    cd(rte_rrtmgp_rfmip_clear_sky) do
        run(`make`)
        if !traverse_input_params
            # stage_files.py calls urllib.request.urlretrieve, so call sparingly (can get blocked out by too many calls)
            run(`python stage_files.py`)
            run(`python run-rfmip-examples.py`)
        else
            println("Running RFMIP drivers")
            conds_file = joinpath(".", "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc")
            lw_coeffs = joinpath(rte_rrtmgp_dir, "rrtmgp", "data", "rrtmgp-data-lw-g256-2018-12-04.nc")
            sw_coeffs = joinpath(rte_rrtmgp_dir, "rrtmgp", "data", "rrtmgp-data-sw-g224-2018-12-04.nc")
            for wavelength in ["l", "s"]
                for direction in ["u", "d"]
                    for forcing_index in ["1"]
                        # for forcing_index in ["1","2","3"]
                        for physics_index in ["1"]
                            # for physics_index in ["1","2"]
                            @show wavelength, direction, forcing_index, physics_index
                            # Arguments: block size, input conditions, coefficient files, forcing index, physics index
                            run(
                                `$(joinpath(".", "rrtmgp_rfmip_lw")) 8 $(conds_file) $(lw_coeffs) $(forcing_index) $(physics_index)`,
                            )
                            run(
                                `$(joinpath(".", "rrtmgp_rfmip_sw")) 8 $(conds_file) $(sw_coeffs) $(forcing_index) $(physics_index)`,
                            )
                        end
                    end
                end
            end

        end
    end
end

if run_all_sky
    cd(rte_rrtmgp_rfmip_all_sky) do
        # subprocess.run([all_sky_exe_name, atmos_file, lw_gas_coeffs_file, lw_clouds_coeff_file, ncol_str, nloops_str])
        run(`make`)
        run(`python $(joinpath(".", "run-allsky-example.py"))`)
        # run(`python $(joinpath(".", "compare-to-reference.py"))`)
    end
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
#---------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Model of a Gibbs reactor
#*----------------------------------------------------------------------
#*
#*   Description:
#*       Thermodynamic equilibrium modeling of a reactor using Gibbs
#*	free energy minimization approach.
#*
#*   Assumptions:
#*		* single-phases involved
#*       * thermodynamic equilibrium
#*		* steady-state
#*
#*	Specify:
#*		* inlet stream
#*		* number of elements related to components
#*		* matrix of elements by components
#*		* equilibrium temperature
#*
#*----------------------------------------------------------------------
#* Author: Rodolfo Rodrigues
#* $Id$
#*--------------------------------------------------------------------
#---------------------------------------------------------------------
#*	only vapour phase
#*--------------------------------------------------------------------
type gibbs_vap
	gibbs_vap()=begin
		NElem=outers.NElem
		new(
			tank_vap(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of elements",
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Universal gas constant",
				:Unit=>"J/mol/K",
				:Default=>8.314
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Number of elements per component"
			)),(NElem,NComp)),
			fill(pressure(Dict{Symbol,Any}(
				:Brief=>"Fugacity in standard state",
				:Default=>1,
				:DisplayUnit=>"atm"
			)),(NComp)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Reference temperature",
				:Default=>298.15
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			fill(energy_mol(Dict{Symbol,Any}(
				:Brief=>"Gibbs free-energy change of formation"
			)),(NComp)),
			fill(energy_mol(Dict{Symbol,Any}(
				:Brief=>"Lagrangian multiplier",
				:Symbol=>"\\lambda"
			)),(NElem)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Activity",
				:Symbol=>"\\hat{a}",
				:Lower=>0
			)),(NComp)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Overall component rate of reaction"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fractional conversion of component",
				:Symbol=>"X",
				:Default=>0
			)),(NComp)),
			fill(flow_mol(Dict{Symbol,Any}(
				:Brief=>"Component molar flow rate"
			)),(NComp)),
			[
				:(Outlet.F*Outlet.z = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z + rate*_base_1._base_1.Tank.V),
				:(Outlet.P = _base_1._base_1.Outletm.P),
				:(Outlet.F = sum(Fi)),
				:(Fi = Outlet.F*Outlet.z),
				:(Outlet.F*Outlet.h = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.h),
				:(sumt(Fi*na) = sumt(_base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z*na)),
				:(G = _base_1._base_1.PP.IdealGasGibbsOfFormation(Outlet.T)),
				:(G([1:_base_1._base_1.NComp]) + sumt(lambda*na(:,[1:_base_1._base_1.NComp])) = -Rg*Outlet.T*ln(activ([1:_base_1._base_1.NComp]))),
				:(Fi([1:_base_1._base_1.NComp]) = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z([1:_base_1._base_1.NComp])*(1 - conv([1:_base_1._base_1.NComp]))),
				:(conv([1:_base_1._base_1.NComp]) = 1),
				:(conv([1:_base_1._base_1.NComp]) = 0),
				:(activ = _base_1._base_1.PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z) *Outlet.P*Outlet.z/fs),
			],
			[
				"Outlet stream","Mechanical equilibrium","Steady-state","Component molar flow rate","Energy balance","Element balance","Gibbs free-energy of formation","Lagrangian multiplier","Molar conversion","Molar conversion","Molar conversion","Activity",
			],
			[:NElem,:Rg,:na,:fs,:To,],
			[:Outlet,:G,:lambda,:activ,:rate,:conv,:Fi,]
		)
	end
	_base_1::tank_vap
	NElem::DanaInteger
	Rg::DanaReal
	na::Array{DanaReal}
	fs::Array{pressure}
	To::temperature
	Outlet::vapour_stream
	G::Array{energy_mol}
	lambda::Array{energy_mol}
	activ::Array{DanaReal}
	rate::Array{reaction_mol}
	conv::Array{DanaReal}
	Fi::Array{flow_mol}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export gibbs_vap
function setEquationFlow(in::gibbs_vap)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	#	"Gibbs free-energy of formation without Cp correction"
	#	G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
	#		+ PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
	addEquation(8)
	if (Outletm.z(i) > 1e-16)
		addEquation(9)
	else
		if (Outlet.z(i) > 0)
			addEquation(10)
			# ?

		else
			addEquation(11)
			# ?

		end
	end
	addEquation(12)
end
function atributes(in::gibbs_vap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/cstr"
	fields[:Brief]="Model of a generic vapour-phase Gibbs CSTR"
	fields[:Info]="
== Assumptions ==
* thermodynamic equilibrium
* steady-state

== Specify ==
* inlet stream
* number of elements related to components
* matrix of elements by components
* equilibrium temperature
"
	drive!(fields,_)
	return fields
end
gibbs_vap(_::Dict{Symbol,Any})=begin
	newModel=gibbs_vap()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Combinatorics

export expval_as_corrs, corr_as_expvals

function term2corr(A::QuTerm,C::Tuple)
    corrs = noaliascopy(A.corrs)
    sizehint!(corrs,length(corrs)+length(C))
    for grp in C
        cc = Corr(BaseOpProduct(collect(map(i -> A.bares.v[i], grp))))
        push!(corrs,cc)
    end
    sort!(corrs)
    QuTerm(A.nsuminds,A.δs,A.params,A.expvals,corrs,BaseOpProduct())
end

"""
    expval_as_corrs(expr::QuExpr)

Take an expression `expr=A B C + D E...` and write its expectation value in
terms of correlations ``⟨A⟩_c, ⟨B⟩_c, ⟨AB⟩_c, ⟨ABC⟩_c, \\ldots``. Note that
``⟨A⟩_c = ⟨A⟩``.

E.g., `expval_as_corrs(adag(:n)*a(:n))` returns ``⟨a^\\dagger_n a_n⟩_c +
⟨a^\\dagger_n⟩_c ⟨a_n⟩_c`` (which is equal to ``⟨a^\\dagger_n a_n⟩``), while
`expval_as_corrs(adag(:n)*a(:m)*a(:n))` returns ``\\langle a_{n}^\\dagger a_{m}
a_{n} \\rangle_{c} + \\langle a_{n}^\\dagger \\rangle_{c} \\langle a_{m}
\\rangle_{c} \\langle a_{n} \\rangle_{c} + \\langle a_{n}^\\dagger \\rangle_{c}
\\langle a_{m} a_{n} \\rangle_{c} + \\langle a_{m} \\rangle_{c} \\langle
a_{n}^\\dagger a_{n} \\rangle_{c} + \\langle a_{n} \\rangle_{c} \\langle
a_{n}^\\dagger a_{m} \\rangle_{c}``.

See also: [`expval`](@ref), [`corr`](@ref)"""
function expval_as_corrs(A::QuExpr)
    newA = QuExpr()
    corrterms = QuExpr()
    for (t,s) in A.terms
        # first calculate the correlation for the term in the sum with the "bare" indices, which means that the sum index
        # is assumed to be distinct from the indices of the expressions
        # then, calculate the correction expression for the term where the sum index is the same as any of the non-summed indices,
        # and add the correction between that and the "incorrect" term where index identity is not taken into account
        _add_corrs!(newA,t,s)

        if t.nsuminds>0 && !isempty(t.bares)
            extinds = unique(extindices(t.bares))
            for ii in 1:t.nsuminds
                sumind = sumindex(ii)
                for ind in extinds
                    f = replace_inds(sumind=>ind,(sumindex(jj)=>sumindex(jj-1) for jj=ii+1:t.nsuminds)...)
                    tb = f(t.bares)
                    tbc = noaliascopy(tb)
                    normal_order!(tb,corrterms)
                    @assert isempty(corrterms)
                    if tb != tbc
                        tright = _normalize_without_commutation(QuTerm(t.nsuminds-1,f.(t.δs),f.(t.params),f.(t.expvals),f.(t.corrs),tb))
                        _add_corrs!(newA,tright,s)
                        twrong = QuTerm(tright.nsuminds,tright.δs,tright.params,tright.expvals,tright.corrs,tbc)
                        _add_corrs!(newA,twrong,-s,_add_with_normal_order!)
                    end
                end
            end
        end
    end
    newA
end

function _add_corrs!(A,t,s,addfun! = _add_sum_term!)
    if isempty(t.bares)
        addfun!(A,t,s)
    else
        for C in expval2corrs_inds(length(t.bares))
            addfun!(A,term2corr(t,C),s)
        end
    end
end

function _add_expvs!(A,t,s,addfun! = _add_sum_term!)
    if isempty(t.bares)
        addfun!(A,t,s)
    else
        for (C,sC) in corr2expvals_inds(length(t.bares))
            addfun!(A,term2expvals(t,C),s*sC)
        end
    end
end

"""
    CorrTup_isless(a,b)

isless for Tuples of integers that represent Corr of sorted operators (with n representing An such that n<m == An<Am).
(n,m,...) ≡ Corr(An*Am*...). Defined in such a way that the same order is obtained as with `BaseOperator` objects
"""
CorrTup_isless(a::NTuple{N,Int},b::NTuple{M,Int}) where {N,M} = N==M ? a<b : N<M

"""
    CorrPerm_isless(a,b)

isless for Tuples of Tuples representing products of Corr (see above) of a permutation of operators.
E.g., ((1,3),(2,)) represents <A1 A3>_C <A2>_C.
It is assumed that the total number of operators in a and b is equal, i.e., that `sum(length.(a)) == sum(length.(b))`."""
function CorrPerm_isless(a::Tuple,b::Tuple)
    for (ca,cb) in zip(a,b)
        ca == cb || return CorrTup_isless(ca,cb)
    end
    # if we reach here, they are equal
    return false
end

"""
    expval2corrs_inds(N::Int)

for N operators, create an array of tuples of tuples that represents the terms
in a sum of products of correlators. Each tuple corresponds to a sum term, see
[`CorrTup_isless`](@ref) and [`CorrPerm_isless`](@ref) for details of the
format. The returned array and terms are sorted such that if the N operators are
sorted, the represented expression is also sorted with the conventions of the
QuantumAlgebra package. This allows to directly return a normal-ordered form.
"""
function expval2corrs_inds(N::Int)
    # get return value for N from cache if present, otherwise calculate and cache it
    get!(_EXPVAL2CORRS_CACHE,N) do
        terms = Any[Tuple(tuple.(1:N))]
        for n = 2:N-1
            append!(terms,ncomb_inds(n,1:N))
        end
        [(tuple(1:N...),),sort!(terms,lt=CorrPerm_isless)...]
    end
end
# preload the cache for N=0 and N=1 so we do not have to special-case those above
const _EXPVAL2CORRS_CACHE = Dict{Int,Vector{Any}}(0 => [((),)], 1 => [((1,),)])

function ncomb_inds(n,inds,used_combs=Set())
    terms = []
    for c in combinations(inds,n)
        C = Tuple(c)
        C in used_combs && continue
        push!(used_combs,C)
        notCinds = setdiff(inds,c)
        S = tuple.(notCinds)
        push!(terms,Tuple(sort!([C,S...],lt=CorrTup_isless)))
        if n<length(inds)
            # pass a copy of used_combs here so getting, e.g., <jk>*C here does not prevent a term <jk> <a><b> later
            for CC in ncomb_inds(n,notCinds,copy(used_combs))
                push!(terms,Tuple(sort!([C,CC...],lt=CorrTup_isless)))
            end
        end
    end
    terms
end

function term2expvals(A::QuTerm,C::Tuple)
    expvals = noaliascopy(A.expvals)
    sizehint!(expvals,length(expvals)+length(C))
    for grp in C
        cc = ExpVal(BaseOpProduct(collect(map(i -> A.bares.v[i], grp))))
        push!(expvals,cc)
    end
    sort!(expvals)
    QuTerm(A.nsuminds,A.δs,A.params,expvals,A.corrs,BaseOpProduct())
end

_canon_op_to_num(A::BaseOperator) = _canon_ind_to_num(only(A.inds))
_expval2tuple(A::ExpVal) = Tuple(_canon_op_to_num.(A.ops.v))
const _CORR2EXPVALS_CACHE = Dict{Int,Vector{Any}}()
function corr2expvals_inds(N)
    get!(_CORR2EXPVALS_CACHE,N) do
        A = QuExpr(QuTerm(BaseOpProduct(BosonDestroy.(:a,_canon_ind.(1:N)))))
        # expval_as_corrs(As...) = <As...> in terms of correlators <A1 A2...>c
        # <As...>c = <As...> - <As...> + <As...>c, where RHS expression only contains <As...> and lower-order correlators
        cA = expval(A) - expval_as_corrs(A) + corr(A)
        # rewrite the lower-order correlators in terms of expectation values
        ncA = QuExpr()
        for (t,s) in cA.terms
            nt = noaliascopy(t)
            empty!(nt.corrs)
            nts = QuExpr(nt)
            for co in t.corrs
                nts *= QuExpr((term2expvals(QuTerm(co.ops),C)=>sC for (C,sC) in corr2expvals_inds(length(co))))
            end
            for (ntt,ns) in nts.terms
                _add_with_normal_order!(ncA,ntt,ns*s)
            end
        end
        [Tuple(_expval2tuple.(t.expvals))=>s for (t,s) in sort!(collect(ncA.terms),by=x->x[1])]
    end
end

function corr_as_expvals(A::QuExpr)
    newA = QuExpr()
    corrterms = QuExpr()
    for (t,s) in A.terms
        # first calculate the correlation for the term in the sum with the "bare" indices, which means that the sum index
        # is assumed to be distinct from the indices of the expressions
        # then, calculate the correction expression for the term where the sum index is the same as any of the non-summed indices,
        # and add the correction between that and the "incorrect" term where index identity is not taken into account
        _add_expvs!(newA,t,s)

        if t.nsuminds>0 && !isempty(t.bares)
            extinds = unique(extindices(t.bares))
            for ii in 1:t.nsuminds
                sumind = sumindex(ii)
                for ind in extinds
                    f = replace_inds(sumind=>ind,(sumindex(jj)=>sumindex(jj-1) for jj=ii+1:t.nsuminds)...)
                    tb = f(t.bares)
                    tbc = noaliascopy(tb)
                    normal_order!(tb,corrterms)
                    @assert isempty(corrterms)
                    if tb != tbc
                        tright = _normalize_without_commutation(QuTerm(t.nsuminds-1,f.(t.δs),f.(t.params),f.(t.expvals),f.(t.corrs),tb))
                        _add_expvs!(newA,tright,s)
                        twrong = QuTerm(tright.nsuminds,tright.δs,tright.params,tright.expvals,tright.corrs,tbc)
                        _add_expvs!(newA,twrong,-s,_add_with_normal_order!)
                    end
                end
            end
        end
    end
    newA
end


text:
using ESDLZarr
using ESDL
using Zarr
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

c1 = Cube("/BGI/scratch/DataCube/v1.0.0/low-res/")
d1 = getCubeData(c1,variable="air_temperature_2m",longitude=(30,31),latitude=(50,51),time=(Date("2002-01-01"),Date("2008-12-31")))

c2=Cube(zopen("/home/fgans/zarr/esdc-8d-0.25deg-1x720x1440-1.0.1_1_zarr"))
d2 = getCubeData(c2,variable="air_temperature_2m",longitude=(30,31),latitude=(50,51),time=(Date("2002-01-01"),Date("2008-12-31")))

a1 = readCubeData(d1)

a2 = readCubeData(d2)


text:
using DashBootstrapComponents

nav = dbc_nav([
    dbc_navitem(dbc_navlink("Active", active = true, href = "#")),
    dbc_navitem(dbc_navlink("A link", href = "#")),
    dbc_navitem(dbc_navlink("Another link", href = "#")),
    dbc_navitem(dbc_navlink("Disabled", disabled = true, href = "#")),
    dbc_dropdownmenu(
        [dbc_dropdownmenuitem("Item 1"), dbc_dropdownmenuitem("Item 2")],
        label = "Dropdown",
        nav = true,
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
#VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Playground

using Compat
using ArgParse

include("constants.jl")
include("config.jl")
include("utils.jl")
include("parsing.jl")
include("install.jl")
include("create.jl")
include("activate.jl")
include("execute.jl")
include("list.jl")
include("clean.jl")

export
    # methods
    main,
    argparse,
    load_config,
    install,
    dirinstall,
    #gitinstall,
    create,
    activate,
    execute,
    list,
    clean,
    config_path,

    # Constants
    DEFAULT_CONFIG

function main(cmd_args=ARGS, config="", root="")
    if config == "" && root == ""
        config = joinpath(config_path(), "config.yml")
        root = config_path()
    end

    args = argparse(cmd_args)

    cmd = args["%COMMAND%"]
    args = args[cmd]

    config = load_config(config, root)

    if cmd == "install"
        install_cmd = args["%COMMAND%"]
        args = args[install_cmd]

        if install_cmd == "download"
            install(
                config,
                args["version"];
                labels=args["labels"],
            )
        elseif install_cmd == "link"
            dirinstall(
                config,
                abspath(args["dir"]);
                labels=args["labels"],
            )
        # elseif install_cmd == "build"
        #     error("Building from source isn't supported yet.")
        end
    elseif cmd == "create"
        create(
            config;
            dir=args["dir"],
            name=args["name"],
            julia=args["julia-version"],
            reqs_file=args["requirements"],
            reqs_type=args["req-type"],
        )
    elseif cmd == "activate"
        activate(
            config;
            dir=args["dir"],
            name=args["name"],
        )
    elseif cmd == "exec"
        execute(
            config,
            `$(Base.shell_split(args["cmd"]))`;
            dir=args["dir"],
            name=args["name"],
        )
    elseif cmd == "list"
        list(
            config;
            show_links=args["show-links"],
        )
    elseif cmd == "clean"
        clean(config)
    elseif cmd == "rm"
        rm(
            config;
            name=args["name"],
            dir=args["dir"],
        )
    end
end

end


text:
using HorizonSideRobots
include("../basics.jl")

md"""
        task_1, но с перегородками
"""

 """Решение задачи 14"""
 function solve_problem_14!(r)
    visited = Set()
    (M, N), (I, J) = get_field_size_and_position!(r)
    cells = Set()

    for i in 1:M
        push!(cells, (i, J))
    end
    for j in 1:N
        push!(cells, (I, j))
    end

    mark_some_cells_with_dfs!(r, visited, (I, J), false, cells)

    show(r)
end


text:
using MinimaxCalibratedEBayes
using ECOS

const MCEB = MinimaxCalibratedEBayes

@testset "Check Proposition 1" begin

prior_dbn =  MixtureModel(Normal, [(-2,0.4), (2, 0.4), (0, 1.0)])
target = PriorDensityTarget(1.0)
prior_classes = [GaussianMixturePriorClass(0.4, -3:0.1:3, ECOS.Optimizer),
                 HermitePriorClass(qmax=40, sobolev_order=2.0, sobolev_bound= 1.0, solver=ECOS.Optimizer)]

for prior_class in prior_classes
	Zs_discr = DiscretizedStandardNormalSamples(-4:0.1:4)
	C∞ = 0.02
	Zs_discr_nbhood = set_neighborhood(Zs_discr, prior_dbn; C∞ = C∞)
	delta_tuner= MCEB.FixedDelta(0.001)

	steinminimax_fit = SteinMinimaxEstimator(Zs_discr, prior_class,
								         target, delta_tuner)

    steinminimax_fit_nbhood = SteinMinimaxEstimator(Zs_discr_nbhood, prior_class,
								         target, delta_tuner)

    # Check if Delta is correct
    @test steinminimax_fit_nbhood.δ == delta_tuner.δ

    @test MCEB.get_δ(steinminimax_fit_nbhood.model, recalculate_δ=false) ≈ delta_tuner.δ atol=1e-6
    @test MCEB.get_δ(steinminimax_fit_nbhood.model; recalculate_δ=true ) ≈ delta_tuner.δ atol=1e-6

    # Check statements about bias
   bias_at_g2 = worst_case_bias(steinminimax_fit_nbhood.Q, Zs_discr_nbhood,
                             steinminimax_fit_nbhood.g2, target)[:bias]

   bias_at_g1 = worst_case_bias(steinminimax_fit_nbhood.Q, Zs_discr_nbhood,
							steinminimax_fit_nbhood.g1, target)[:bias]

   worst_case_bias_stein = worst_case_bias(steinminimax_fit_nbhood.Q, Zs_discr_nbhood,
                                        prior_class, target)

   @test 	worst_case_bias_stein.max_bias ≈ -worst_case_bias_stein.min_bias atol=0.001
   @test 	worst_case_bias_stein.max_bias ≈ worst_case_bias_stein.max_abs_bias atol=0.001
   @test 	abs2(worst_case_bias_stein.max_bias) ≈ worst_case_bias_stein.max_squared_bias atol=0.001

   @test 	worst_case_bias_stein.max_bias ≈ bias_at_g1 atol=0.001
   @test 	worst_case_bias_stein.min_bias ≈ bias_at_g2 atol=0.001
   @test   worst_case_bias_stein.max_bias ≈ steinminimax_fit_nbhood.max_bias atol=0.001
   @test   2*worst_case_bias_stein.max_bias ≈ steinminimax_fit_nbhood.ω_δ - steinminimax_fit_nbhood.δ*steinminimax_fit_nbhood.ω_δ_prime atol=0.001
end
end


text:
struct Stepper{A<:Random.AbstractRNG,M<:AbstractModel,S<:AbstractSampler,K}
    rng::A
    model::M
    sampler::S
    kwargs::K
end

Base.iterate(stp::Stepper) = step(stp.rng, stp.model, stp.sampler; stp.kwargs...)
function Base.iterate(stp::Stepper, state)
    return step(stp.rng, stp.model, stp.sampler, state; stp.kwargs...)
end

Base.IteratorSize(::Type{<:Stepper}) = Base.IsInfinite()
Base.IteratorEltype(::Type{<:Stepper}) = Base.EltypeUnknown()

function steps(
    model::AbstractModel,
    sampler::AbstractSampler;
    kwargs...
)
    return steps(Random.GLOBAL_RNG, model, sampler; kwargs...)
end

"""
    steps([rng, ]model, sampler; kwargs...)

Create an iterator that returns samples from the `model` with the Markov chain Monte Carlo
`sampler`.

# Examples

```jldoctest; setup=:(using AbstractMCMC: steps)
julia> struct MyModel <: AbstractMCMC.AbstractModel end

julia> struct MySampler <: AbstractMCMC.AbstractSampler end

julia> function AbstractMCMC.step(rng, ::MyModel, ::MySampler, state=nothing; kwargs...)
           # all samples are zero
           return 0.0, state
       end

julia> iterator = steps(MyModel(), MySampler());

julia> collect(Iterators.take(iterator, 10)) == zeros(10)
true
```
"""
function steps(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler;
    kwargs...
)
    return Stepper(rng, model, sampler, kwargs)
end


text:
module DEMC
    using LinearAlgebra
    using Statistics
    # using Plots
    using Distributed
    using SharedArrays
    using ParallelDataTransfer
    #gr()

    struct MC
        chain::Array{Float64, 3} # parameter population for all generations
        log_obj::Array{Float64, 2} # log obj along the chain
        Xcurrent::Array{Float64, 2} # population
        log_objcurrent::Array{Float64, 1} # log obj values
    end

    struct MCShared
        chain::SharedArray{Float64, 3} # parameter population for all generations
        log_obj::SharedArray{Float64, 2} # log obj along the chain
        Xcurrent::SharedArray{Float64, 2} # population
        log_objcurrent::SharedArray{Float64, 1} # log obj values
    end

    mutable struct DEMCopt
        N::Int
        K::Int
        Ngeneration::Int
        Nblocks::Int
        blockindex::Array
        eps_scale::Array{Float64, 1}
        γ::Float64
        verbose::Bool
        print_step::Int
        T0::Float64
        TN::Float64
        autostop::Symbol
        autostop_every::Int
        autostop_Rhat
    end

    function demcopt(Npar; N=4, K=10, Ngeneration=5000, Nblocks=1, blockindex=[1:Npar], eps_scale=1e-4*ones(Npar), γ=2.38, verbose = true, print_step=100, T0 = 3, TN = 1e-3, autostop = :Rhat, autostop_every = 1000, autostop_Rhat=1.05)
        return DEMCopt(N, K, Ngeneration, Nblocks, blockindex, eps_scale, γ, verbose, print_step, T0, TN, autostop, autostop_every, autostop_Rhat)
    end

    include("demcz.jl")
    include("utils.jl")
    include("demcz_anneal.jl")

end # module


text:
module box

import ..segment: Segment
import ..style: apply_style
import Term: int, chars, remove_markup_open, get_last_valid_str_idx, get_next_valid_str_idx

export get_row, get_title_row
export ASCII,
    ASCII2, ASCII_DOUBLE_HEAD, SQUARE, SQUARE_DOUBLE_HEAD, MINIMAL, MINIMAL_HEAVY_HEAD
export MINIMAL_DOUBLE_HEAD, SIMPLE, SIMPLE_HEAD, SIMPLE_HEAVY, HORIZONTALS, ROUNDED, HEAVY
export HEAVY_EDGE, HEAVY_HEAD, DOUBLE, DOUBLE_EDGE

"""
  loop_last(v::Vector)

  Returns an iterable yielding tuples (is_last, value).
"""
function loop_last(v::Vector)
    is_last = [i == length(v) for i in 1:length(v)]
    return zip(is_last, v)
end

# ---------------------------------------------------------------------------- #
#                                      BOX                                     #
# ---------------------------------------------------------------------------- #
"""
  BoxLine

Stores the characters for a line of a `Box` object.
"""
struct BoxLine
    left::Char
    mid::Char
    vertical::Char
    right::Char
end

"""
  Box

Defines characters to render boxes.

Row names:

┌─┬┐ top
│ ││ head
├─┼┤ head_row
│ ││ mid
├─┼┤ row
├─┼┤ foot_row
│ ││ foot
└─┴┘ bottom

each row is an instance of `BoxLine`
"""
struct Box
    name::String
    top::BoxLine
    head::BoxLine
    head_row::BoxLine
    mid::BoxLine
    row::BoxLine
    foot_row::BoxLine
    foot::BoxLine
    bottom::BoxLine
end

"""
    Box(string)

Construct a `Box` objet out of a box string.
"""
function Box(box_name::String, box::String)
    top, head, head_row, mid, row, foot_row, foot, bottom = split(box, "\n")

    return Box(
        box_name,
        BoxLine(chars(top)...),
        BoxLine(chars(head)...),
        BoxLine(chars(head_row)...),
        BoxLine(chars(mid)...),
        BoxLine(chars(row)...),
        BoxLine(chars(foot_row)...),
        BoxLine(chars(foot)...),
        BoxLine(chars(bottom)...),
    )
end

function Base.show(io::IO, box::Box)
    if io == stdout
        print(io, "Box ($(box.name))\n$(fit(box, [1, 3, 1]))")
    elseif io == stderr
        print(io, "err")
    else
        print(io, "Box\e[2m($(box.name))\e[0m")
    end
end

"""
    get_row(box, [1, 2, 3], :row)

Gets characters for a row of a Box object.

The level Symbol can be used to specify the box level (:top, :footer...)
"""
function get_row(box::Box, widths::Vector{Int}, level::Symbol)::String
    # get the correct level of the box
    level = getfield(box, level)

    parts = [string(level.left)]
    for (last, w) in loop_last(widths)
        segment = last ? level.mid^w * level.right : level.mid^w * level.vertical
        push!(parts, segment)
    end
    return join(parts)
end

"""
  get_title_row(row::Symbol, box::Box, title::Union{Nothing, AbstractString}; <keyword arguments>)

Create a box row with a title string.

Can create both titles in the top and bottom row to produce subtitles.

#Arguments:
- width::Int: width of line
- style::Union{Nothing:  String}: style of line
- title_style::Union{Nothing:  AbstractString}: style of title string
- justify::Symbol=:left: position of title string

See also [`get_row`](@ref).

"""
function get_title_row(
    row::Symbol,
    box::Box,
    title::Union{Nothing,AbstractString};
    width::Int,
    style::Union{Nothing,String},
    title_style::Union{Nothing,AbstractString},
    justify::Symbol = :left,
)
    initial_line = remove_markup_open(Segment(get_row(box, [width], row), style).plain)
    # @info "Getting title row" initial_line

    if isnothing(title)
        return Segment(initial_line, style)
    else

        # get title
        title = apply_style("[$title_style]" * title)
        @assert Segment(title).measure.w < width - 4 "Title too long for panel of width $width: $title, ($(Segment(title).measure.w))"

        # compose title line
        line = getfield(box, row)

        if justify == :left
            cut_start = get_last_valid_str_idx(initial_line, 4)
            pre = Segment(
                Segment(initial_line[1:cut_start], style) * "\e[0m" * " " * title * " "
            )

            post = line.mid^(length(initial_line) - pre.measure.w - 1) * line.right

        elseif justify == :right
            cut_start = get_next_valid_str_idx(initial_line, ncodeunits(initial_line) - 8)
            post = Segment(
                "\e[0m" * " " * title * " " * Segment(initial_line[cut_start:end], style)
            )

            pre = Segment(
                line.left * line.mid^(length(initial_line) - post.measure.w - 1), style
            )

        else  # justify :center
            cutval = int(ncodeunits(initial_line) / 2 - ncodeunits(title) / 2 - 15)

            cut_start = get_last_valid_str_idx(initial_line, cutval)
            # @info width cutval cut_start length(initial_line) length(initial_line) ncodeunits(initial_line)

            pre = Segment(
                Segment(initial_line[1:cut_start], style) * "\e[0m" * " " * title * " "
            )

            post = line.mid^(length(initial_line) - pre.measure.w - 1) * line.right
        end

        return pre * Segment(post, style)
        # return Segment(post, style)
    end
end

"""
  fit(box::Box, widths::Vector{Int})::String

Creates a box.

The box has one of each level type with columns
widths specified by a vector of widhts.
"""
function fit(box::Box, widths::Vector{Int})::String
    strings = [
        get_row(box, widths, :top),
        get_row(box, widths, :head),
        get_row(box, widths, :head_row),
        get_row(box, widths, :mid),
        get_row(box, widths, :row),
        get_row(box, widths, :foot_row),
        get_row(box, widths, :foot),
        get_row(box, widths, :bottom),
    ]
    return join_lines(strings)
end

# ---------------------------------------------------------------------------- #
#                                   Box types                                  #
# ---------------------------------------------------------------------------- #

ASCII = Box(
    "ASCII",
    """
    +--+
    | ||
    |-+|
    | ||
    |-+|
    |-+|
    | ||
    +--+
    """,
)

ASCII2 = Box(
    "ASCII2",
    """
    +-++
    | ||
    +-++
    | ||
    +-++
    +-++
    | ||
    +-++
    """,
)

ASCII_DOUBLE_HEAD = Box(
    "ASCII_DOUBLE_HEAD",
    """
    +-++
    | ||
    +=++
    | ||
    +-++
    +-++
    | ||
    +-++
    """,
)

SQUARE = Box(
    "SQUARE",
    """
    ┌─┬┐
    │ ││
    ├─┼┤
    │ ││
    ├─┼┤
    ├─┼┤
    │ ││
    └─┴┘
    """,
)

SQUARE_DOUBLE_HEAD = Box(
    "SQUARE_DOUBLE_HEAD",
    """
    ┌─┬┐
    │ ││
    ╞═╪╡
    │ ││
    ├─┼┤
    ├─┼┤
    │ ││
    └─┴┘
    """,
)

MINIMAL = Box(
    "MINIMAL",
    """
      ╷
      │
    ╶─┼╴
      │
    ╶─┼╴
    ╶─┼╴
      │
      ╵
    """,
)

MINIMAL_HEAVY_HEAD = Box(
    "MINIMAL_HEAVY_HEAD",
    """
      ╷
      │
    ╺━┿╸
      │
    ╶─┼╴
    ╶─┼╴
      │
      ╵
    """,
)

MINIMAL_DOUBLE_HEAD = Box(
    "MINIMAL_DOUBLE_HEAD",
    """
      ╷
      │
     ═╪
      │
     ─┼
     ─┼
      │
      ╵
    """,
)

SIMPLE = Box(
    "SIMPLE",
    """

     ──

     ──

    """,
)

SIMPLE_HEAD = Box(
    "SIMPLE_HEAD",
    """

     ──

    """,
)

SIMPLE_HEAVY = Box(
    "SIMPLE_HEAVY",
    """

     ━━

     ━━

    """,
)

HORIZONTALS = Box(
    "HORIZONTALS",
    """
     ──

     ──

     ──
     ──

     ──
    """,
)

ROUNDED = Box(
    "ROUNDED",
    """
    ╭─┬╮
    │ ││
    ├─┼┤
    │ ││
    ├─┼┤
    ├─┼┤
    │ ││
    ╰─┴╯
    """,
)

HEAVY = Box(
    "HEAVY",
    """
    ┏━┳┓
    ┃ ┃┃
    ┣━╋┫
    ┃ ┃┃
    ┣━╋┫
    ┣━╋┫
    ┃ ┃┃
    ┗━┻┛
    """,
)

HEAVY_EDGE = Box(
    "HEAVY_EDGE",
    """
    ┏━┯┓
    ┃ │┃
    ┠─┼┨
    ┃ │┃
    ┠─┼┨
    ┠─┼┨
    ┃ │┃
    ┗━┷┛
    """,
)

HEAVY_HEAD = Box(
    "HEAVY_HEAD",
    """
    ┏━┳┓
    ┃ ┃┃
    ┡━╇┩
    │ ││
    ├─┼┤
    ├─┼┤
    │ ││
    └─┴┘
    """,
)

DOUBLE = Box(
    "DOUBLE",
    """
    ╔═╦╗
    ║ ║║
    ╠═╬╣
    ║ ║║
    ╠═╬╣
    ╠═╬╣
    ║ ║║
    ╚═╩╝
    """,
)

DOUBLE_EDGE = Box(
    "DOUBLE_EDGE",
    """
    ╔═╤╗
    ║ │║
    ╟─┼╢
    ║ │║
    ╟─┼╢
    ╟─┼╢
    ║ │║
    ╚═╧╝
    """,
)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: personalize_events
using AWS.Compat
using AWS.UUIDs
"""
    PutEvents()

Records user interaction event data. For more information see event-record-api.

# Required Parameters
- `eventList`: A list of event data from the session.
- `sessionId`: The session ID associated with the user's visit. Your application generates the sessionId when a user first visits your website or uses your application. Amazon Personalize uses the sessionId to associate events with the user before they log in. For more information see event-record-api.
- `trackingId`: The tracking ID for the event. The ID is generated by a call to the CreateEventTracker API.

# Optional Parameters
- `userId`: The user associated with the event.
"""

put_events(eventList, sessionId, trackingId; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/events", Dict{String, Any}("eventList"=>eventList, "sessionId"=>sessionId, "trackingId"=>trackingId); aws_config=aws_config)
put_events(eventList, sessionId, trackingId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/events", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("eventList"=>eventList, "sessionId"=>sessionId, "trackingId"=>trackingId), args)); aws_config=aws_config)

"""
    PutItems()

Adds one or more items to an Items dataset. For more information see importing-items.

# Required Parameters
- `datasetArn`: The Amazon Resource Number (ARN) of the Items dataset you are adding the item or items to.
- `items`: A list of item data.

"""

put_items(datasetArn, items; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/items", Dict{String, Any}("datasetArn"=>datasetArn, "items"=>items); aws_config=aws_config)
put_items(datasetArn, items, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/items", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("datasetArn"=>datasetArn, "items"=>items), args)); aws_config=aws_config)

"""
    PutUsers()

Adds one or more users to a Users dataset. For more information see importing-users.

# Required Parameters
- `datasetArn`: The Amazon Resource Number (ARN) of the Users dataset you are adding the user or users to.
- `users`: A list of user data.

"""

put_users(datasetArn, users; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/users", Dict{String, Any}("datasetArn"=>datasetArn, "users"=>users); aws_config=aws_config)
put_users(datasetArn, users, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = personalize_events("POST", "/users", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("datasetArn"=>datasetArn, "users"=>users), args)); aws_config=aws_config)