text:
function jump_max_area_perimeter_Qn(n)
    if mod(log2(n),1)==0 && n>=4
        d = pi/4-asin(cos(pi/n)/sqrt(2));
        Ll = 2*n*sin(pi/(2*n))*cos(d/2);
        Lu = 2*n*sin(pi/(2*n));
        model = Model(Ipopt.Optimizer)
        #model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe))
        # Variables
        @variable(model, a[i = 1:Int(n/2)] >= 0, start = pi/n+(-1)^i*d)
        # Perimeter
        @NLobjective(model, Max, sum(4*sin(a[i]/2) for i in 1:Int(n/2)))
        # Optimal diameter graph
        @constraint(model, sum(a[i] for i in 1:Int(n/2)) == pi/2)
        @NLconstraint(model, sum((-1)^(i-1)*sin(sum(a[j] for j in 1:i)) for i in 1:Int(n/2-1)) == 1/2)
        @constraint(model, a[1] <= pi/6)
        @constraint(model, [i = 2:Int(n/2)], a[i] <= pi/3)
        # Bounds
        @NLconstraint(model, Ll <= sum(4*sin(a[i]/2) for i in 1:Int(n/2)) <= Lu)
        # Solving
        optimize!(model)
        L = objective_value(model)
        a = value.(a)
        ct = solve_time(model)
        return(L,a,ct)
    end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
module RayTracingWeekend

# Write your package code here.
using Images

end

ext: jl


text:
using AdaptiveResonance
using Test
using Logging
using DelimitedFiles

# Set the log level
LogLevel(Logging.Info)

# Auxiliary generic functions for loading data, etc.
include("test_utils.jl")

# Load the data and test across all supervised modules
data = load_iris("../data/Iris.csv")

@testset "common.jl" begin
    @info "------- Common Code Tests -------"
    # Example arrays
    three_by_two = [1 2; 3 4; 5 6]

    # Test DataConfig constructors
    @info "Testing DataConfig..."
    dc1 = DataConfig()                  # Default constructor
    dc2 = DataConfig(0, 1, 2)           # When min and max are same across all features
    dc3 = DataConfig([0, 1], [2, 3])    # When min and max differ across features
    dc4 = DataConfig(three_by_two)      # When a data matrix is provided

    # Test get_n_samples
    @info "Testing get_n_samples..."
    @test get_n_samples([1,2,3]) == 1           # 1-D array case
    @test get_n_samples(three_by_two) == 2      # 2-D array case

    # Test data_setup!
    @info "Testing data_setup!..."
    data_setup!(DDVFA(), three_by_two)
    data_setup!(DDVFA().config, three_by_two)

    # Test breaking situations
    @info "Testing common code error handling..."
    @test_throws ErrorException performance([1,2],[1,2,3])
    @test_logs (:warn,) AdaptiveResonance.data_setup!(dc3, three_by_two)
    bad_config =  DataConfig(1, 0, 3)
    @test_throws ErrorException linear_normalization(three_by_two, config=bad_config)
end # @testset "common.jl"

@testset "constants.jl" begin
    @info "------- Constants Tests -------"
    ddvfa_methods = [
        "single",
        "average",
        "complete",
        "median",
        "weighted",
        "centroid"
    ]
    @test AdaptiveResonance.DDVFA_METHODS == ddvfa_methods
end # @testset "constants.jl"

@testset "AdaptiveResonance.jl" begin
    # Module loading
    include("modules.jl")
end # @testset "AdaptiveResonance.jl"

@testset "Train Test" begin
    # All ART modules
    arts = [
        FuzzyART,
        DVFA,
        DDVFA,
        SFAM,
        DAM,
    ]
    n_arts = length(arts)

    # All common ART options
    art_opts = [
        (display = true,),
        # (display = false,),
    ]

    # Specific ART options
    art_specifics = Dict(
        DDVFA => [
            (gamma_normalization=true,),
            (gamma_normalization=false,),
        ],
        FuzzyART => [
            (gamma_normalization=true,),
            (gamma_normalization=false,),
        ],
    )

    # All test option permutations
    test_opts = [
        (get_bmu = true,),
        (get_bmu = false,)
    ]
    n_test_opts = length(test_opts)

    @info "-------------- BEGIN TRAIN TEST --------------"
    # Performance baseline for all algorithms
    perf_baseline = 0.7

    # Iterate over all ART modules
    for ix = 1:n_arts
        # Iterate over all test options
        for jx = 1:n_test_opts
            # If we are testing a module with different options, merge
            if haskey(art_specifics, arts[ix])
                local_art_opts = vcat(art_opts, art_specifics[arts[ix]])
            else
                local_art_opts = art_opts
            end
            # Iterate over all options
            for kx = 1:length(local_art_opts)
                # Only do the unsupervised method if we have an ART module (not ARTMAP)
                if arts[ix] isa ART
                    # Unsupervised
                    train_test_art(arts[ix](;local_art_opts[kx]...), data; test_opts=test_opts[jx])
                end
                # Supervised
                @test train_test_art(arts[ix](;local_art_opts[kx]...), data; supervised=true, test_opts=test_opts[jx]) >= perf_baseline
            end
        end
    end

    @info "-------------- END TRAIN TEST --------------"
end # @testset "Train Test"

@testset "kwargs" begin
    @info "--------- KWARGS TEST ---------"

    arts = [
        FuzzyART,
        DVFA,
        DDVFA,
        SFAM,
        DAM
    ]

    for art in arts
        art_module = art(alpha=1e-3, display=false)
    end

    @info "--------- END KWARGS TEST ---------"
end # @testset "kwargs"

@testset "FuzzyART" begin
    @info "------- FuzzyART Testing -------"

    # FuzzyART train and test
    my_FuzzyART = FuzzyART()
    # local_complement_code = AdaptiveResonance.complement_code(data.train_x)
    # train!(my_FuzzyART, local_complement_code, preprocessed=true)
    train!(my_FuzzyART, data.train_x)

    # Similarity methods
    methods = [
        "single",
        "average",
        "complete",
        "median",
        "weighted",
        "centroid"
    ]

    # Both field names
    field_names = ["T", "M"]

    # Compute a local sample for FuzzyART similarity method testing
    # local_sample = local_complement_code[:, 1]
    # local_complement_code = AdaptiveResonance.complement_code(data.train_x)
    # local_sample = data.train_x[:, 1]
    local_sample = AdaptiveResonance.complement_code(data.train_x[:, 1], config=my_FuzzyART.config)

    # Compute the local activation and match
    # AdaptiveResonance.activation_match!(my_FuzzyART, local_sample)

    # # Declare the true activation and match magnitudes
    # truth = Dict(
    #     "single" => Dict(
    #         "T" => 0.9988714513100155,
    #         "M" => 2.6532834139109758
    #     ),
    #     "average" => Dict(
    #         "T" => 0.33761483787933894,
    #         "M" => 1.1148764060015297
    #     ),
    #     "complete" => Dict(
    #         "T" => 0.018234409874338647,
    #         "M" => 0.07293763949735459
    #     ),
    #     "median" => Dict(
    #         "T" => 0.2089217851518073,
    #         "M" => 0.835687140607229
    #     ),
    #     "weighted" => Dict(
    #         "T" => 0.5374562506748786,
    #         "M" => 1.4396083090159748
    #     ),
    #     "centroid" => Dict(
    #         "T" => 0.0,
    #         "M" => 0.0
    #     )
    # )

    # # Test every method and field name
    # for method in methods
    #     results = Dict()
    #     for field_name in field_names
    #         results[field_name] = AdaptiveResonance.similarity(method, my_FuzzyART, field_name, local_sample, my_FuzzyART.opts.gamma_ref)
    #         @test isapprox(truth[method][field_name], results[field_name])
    #     end
    #     @info "Method: $method" results
    # end

    # Check the error handling of the similarity function
    # Access the wrong similarity metric keyword ("asdf")
    @test_throws ErrorException AdaptiveResonance.similarity("asdf", my_FuzzyART, "T", local_sample, my_FuzzyART.opts.gamma_ref)
    # Access the wrong output function ("A")
    @test_throws ErrorException AdaptiveResonance.similarity("centroid", my_FuzzyART, "A", local_sample, my_FuzzyART.opts.gamma_ref)

end # @testset "FuzzyART"

@testset "ARTSCENE.jl" begin
    # ARTSCENE training and testing
    include("test_artscene.jl")
end # @testset "ARTSCENE.jl"


text:
struct ErrorBarGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    tag::Symbol
end
ErrorBarGeometry(default_statistic=Gadfly.Stat.identity(); tag=empty_tag) =
    ErrorBarGeometry(default_statistic, tag)

struct XErrorBarGeometry <: Gadfly.GeometryElement
    tag::Symbol
end
XErrorBarGeometry(; tag=empty_tag) = XErrorBarGeometry(tag)

struct YErrorBarGeometry <: Gadfly.GeometryElement
    tag::Symbol
end
YErrorBarGeometry(; tag=empty_tag) = YErrorBarGeometry(tag)

"""
    Geom.errorbar

Draw vertical error bars if the `x`, `ymin`, and `ymax` aesthetics
are specified and/or horizontal error bars for `y`, `xmin`, and `xmax`.
Optionally color them with `color`.

See also [`Geom.xerrorbar`](@ref) and [`Geom.yerrorbar`](@ref).

"""
const errorbar = ErrorBarGeometry

"""
    Geom.xerrorbar

Draw horizontal error bars at `y` from `xmin` to `xmax`.  Optionally
color them with `color`.
"""
const xerrorbar = XErrorBarGeometry

"""
    Geom.yerrorbar

Draw vertical error bars at `x` from `ymin` to `ymax`.  Optionally
color them with `color`.
"""
const yerrorbar = YErrorBarGeometry

element_aesthetics(::ErrorBarGeometry) = [:x, :y, :xmin, :xmax, :ymin, :ymax, :color]
element_aesthetics(::XErrorBarGeometry) = [:y, :xmin, :xmax, :color]
element_aesthetics(::YErrorBarGeometry) = [:x, :ymin, :ymax, :color]
default_statistic(geom::ErrorBarGeometry) = geom.default_statistic

# Generate a form for the errorbar geometry.
#
# Args:
#   geom: errorbar geometry
#   theme: the plot's theme
#   aes: aesthetics
#
# Returns:
#   A compose Form.
#
function render(geom::ErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    # check for X and Y error bar aesthetics
    if isempty(Gadfly.undefined_aesthetics(aes, :y, :xmin, :xmax))
        xctx = render(xerrorbar(), theme, aes)
    else
        xctx = nothing
    end
    if isempty(Gadfly.undefined_aesthetics(aes, :x, :ymin, :ymax))
        yctx = render(yerrorbar(), theme, aes)
    else
        yctx = nothing
    end
    compose(context(order=3), xctx, yctx, svgclass("geometry"))
end

function render(geom::YErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.errorbar", aes, :x, :ymin, :ymax)
    Gadfly.assert_aesthetics_equal_length("Geom.errorbar", aes, :x, :ymin, :ymax)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = RGB{Float32}[theme.default_color]
    aes = inherit(aes, default_aes)
    caplen = theme.errorbar_cap_length/2
    ttc, teb, tbc = subtags(geom.tag, :top_cap, :error_bar, :bottom_cap)
    cs = stroke(theme.stroke_color.(aes.color))

    ctx = compose!(
        context(order=3; tag=geom.tag),

        # top cap
        (context(),
            Compose.line([[(x*cx - caplen, ymax), (x*cx + caplen, ymax)]
                      for (x, ymax) in zip(aes.x, aes.ymax)], ttc), cs,
            svgclass("marker")),

        # error bar
        Compose.line([[(x*cx, ymax), (x*cx, ymin)]
                      for (x, ymin, ymax) in zip(aes.x, aes.ymin, aes.ymax)], teb), cs,

        # bottom cap
        (context(),
            Compose.line([[(x*cx - caplen, ymin), (x*cx + caplen, ymin)]
                      for (x, ymin) in zip(aes.x, aes.ymin)], tbc), cs,
            svgclass("marker")),

        linewidth(theme.line_width))

    (aes.color_key_continuous == true || aes.color == nothing) || compose!(ctx,
            svgclass([svg_color_class_from_label(aes.color_label([c])[1]) for c in aes.color]))

    return ctx
end

function render(geom::XErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.errorbar", aes, :y, :xmin, :xmax)
    Gadfly.assert_aesthetics_equal_length("Geom.errorbar", aes, :y, :xmin, :xmax)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = RGB{Float32}[theme.default_color]
    aes = inherit(aes, default_aes)
    caplen = theme.errorbar_cap_length/2
    tlc, teb, trc = subtags(geom.tag, :left_cap, :error_bar, :right_cap)
    cs = stroke(theme.stroke_color.(aes.color))

    ctx = compose!(
        context(order=3, tag=geom.tag),

        # left cap
        (context(),
            Compose.line([[(xmin, y*cy - caplen), (xmin, y*cy + caplen)]
                      for (xmin, y) in zip(aes.xmin, aes.y)], tlc), cs,
            svgclass("marker")),

        # error bar
        Compose.line([[(xmin, y*cy), (xmax, y*cy)]
                      for (xmin, xmax, y) in zip(aes.xmin, aes.xmax, aes.y)], teb), cs,

        # right cap
        (context(),
            Compose.line([[(xmax, y*cy - caplen), (xmax, y*cy + caplen)]
                      for (xmax, y) in zip(aes.xmax, aes.y)], trc), cs,
            svgclass("marker")),

        linewidth(theme.line_width))

    (aes.color_key_continuous == true || aes.color == nothing) || compose!(ctx,
            svgclass([svg_color_class_from_label(aes.color_label([c])[1])
                      for c in aes.color]))

    return ctx
end


text:
"""
`module  JAC.PhotoIonization`
    ... a submodel of JAC that contains all methods for computing photoionization properties between some initial
        and final-state multiplets.
"""
module PhotoIonization

    using Printf, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..Radial, ..Nuclear, ..ManyElectron, ..PhotoEmission,
                  ..PlasmaShift, ..TableStrings

    """
    `struct  PhotoIonization.Settings  <:  AbstractProcessSettings`  ... defines a type for the details and parameters of computing photoionization lines.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + photonEnergies          ::Array{Float64,1}             ... List of photon energies.
        + calcAnisotropy          ::Bool                         ... True, if the beta anisotropy parameters are to be calculated and false otherwise (o/w).
        + calcPartialCs           ::Bool                         ... True, if partial cross sections are to be calculated and false otherwise.
        + calcTensors             ::Bool                         ... True, if statistical tensors of the excited atom are to be calculated and false o/w.
        + printBefore             ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection           ::LineSelection                ... Specifies the selected levels, if any.
        + stokes                  ::ExpStokes                    ... Stokes parameters of the incident radiation.
    """
    struct Settings  <:  AbstractProcessSettings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        calcAnisotropy            ::Bool
        calcPartialCs             ::Bool
        calcTensors               ::Bool
        printBefore               ::Bool
        lineSelection             ::LineSelection
        stokes                    ::ExpStokes
    end

    """
    `PhotoIonization.Settings()`  ... constructor for the default values of photoionization line computations
    """
    function Settings()
        Settings(Basics.EmMultipole[E1], Basics.UseGauge[Basics.UseCoulomb, Basics.UseBabushkin], Float64[], false, false, false, false,
                 LineSelection(), Basics.ExpStokes())
    end

    """
    `PhotoIonization.Settings(set::PhotoIonization.Settings;`

            multipoles=..,          gauges=..,                  photonEnergies=..,          calcAnisotropy=..,
            calcPartialCs0..,       calcTensors=..,             printBefore=..,             lineSelection=..,
            stokes=..)

        ... constructor for modifying the given PhotoIonization.Settings by 'overwriting' the previously selected parameters.
    """
    function Settings(set::PhotoIonization.Settings;
        multipoles::Union{Nothing,Array{EmMultipole,1}}=nothing,                gauges::Union{Nothing,Array{UseGauge,1}}=nothing,
        photonEnergies::Union{Nothing,Array{Float64,1}}=nothing,                calcAnisotropy::Union{Nothing,Bool}=nothing,
        calcPartialCs::Union{Nothing,Bool}=nothing,                             calcTensors::Union{Nothing,Bool}=nothing,
        printBefore::Union{Nothing,Bool}=nothing,                               lineSelection::Union{Nothing,LineSelection}=nothing,
        stokes::Union{Nothing,ExpStokes}=nothing)

        if  multipoles      == nothing   multipolesx      = set.multipoles        else  multipolesx      = multipoles       end
        if  gauges          == nothing   gaugesx          = set.gauges            else  gaugesx          = gauges           end
        if  photonEnergies  == nothing   photonEnergiesx  = set.photonEnergies    else  photonEnergiesx  = photonEnergies   end
        if  calcAnisotropy  == nothing   calcAnisotropyx  = set.calcAnisotropy    else  calcAnisotropyx  = calcAnisotropy   end
        if  calcPartialCs   == nothing   calcPartialCsx   = set.calcPartialCs     else  calcPartialCsx   = calcPartialCs    end
        if  calcTensors     == nothing   calcTensorsx     = set.calcTensors       else  calcTensorsx     = calcTensors      end
        if  printBefore     == nothing   printBeforex     = set.printBefore       else  printBeforex     = printBefore      end
        if  lineSelection   == nothing   lineSelectionx   = set.lineSelection     else  lineSelectionx   = lineSelection    end
        if  stokes          == nothing   stokesx          = set.stokes            else  stokesx          = stokes           end

        Settings( multipolesx, gaugesx, photonEnergiesx, calcAnisotropyx, calcPartialCsx, calcTensorsx, printBeforex,
                  lineSelectionx, stokesx)
    end

    # `Base.show(io::IO, settings::PhotoIonization.Settings)`  ... prepares a proper printout of the variable settings::PhotoIonization.Settings.
    function Base.show(io::IO, settings::PhotoIonization.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "calcAnisotropy:           $(settings.calcAnisotropy)  ")
        println(io, "calcPartialCs:            $(settings.calcPartialCs)  ")
        println(io, "calcTensors:              $(settings.calcTensors)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "lineSelection:            $(settings.lineSelection)  ")
        println(io, "stokes:                   $(settings.stokes)  ")
    end

    """
    `struct  PhotoIonization.Channel`
        ... defines a type for a photoionization channel to help characterize a single multipole and scattering (continuum) state
            of many electron-states with a single free electron.

        + multipole      ::EmMultipole          ... Multipole of the photon absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... Photoionization amplitude associated with the given channel.
   """
    struct  Channel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Line`  ... defines a type for a photoionization line that may include the definition of channels.

        + initialLevel   ::Level                  ... initial-(state) level
        + finalLevel     ::Level                  ... final-(state) level
        + electronEnergy ::Float64                ... Energy of the (outgoing free) electron.
        + photonEnergy   ::Float64                ... Energy of the absorbed photon.
        + crossSection   ::EmProperty             ... Cross section for this photoionization.
        + hasChannels    ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                      free-electron energy, kappa, multipole, etc., or not.
        + channels       ::Array{PhotoIonization.Channel,1}  ... List of PhotoIonization.Channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        electronEnergy   ::Float64
        photonEnergy     ::Float64
        crossSection     ::EmProperty
        hasChannels      ::Bool
        channels         ::Array{PhotoIonization.Channel,1}
    end

    """
    `PhotoIonization.Line(initialLevel::Level, finalLevel::Level, crossSection::Float64)`
        ... constructor for an photoionization line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, crossSection::EmProperty)
        Line(initialLevel, finalLevel, totalRate, 0., 0., crossSection, false, PhotoChannel[] )
    end

    # `Base.show(io::IO, line::PhotoIonization.Line)`  ... prepares a proper printout of the variable line::PhotoIonization.Line.
    function Base.show(io::IO, line::PhotoIonization.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "electronEnergy:    $(line.electronEnergy)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
    end

    """
    `PhotoIonization.amplitude(kind::String, channel::PhotoIonization.Channel, energy::Float64, continuumLevel::Level,
                                   initialLevel::Level, grid::Radial.Grid)`
        ... to compute the kind = (photoionization) amplitude  <(alpha_f J_f, epsilon kappa) J_t || O^(photoionization) || alpha_i J_i>
            due to the electron-photon interaction for the given final and initial level, the partial wave of the outgoing
            electron as well as the given multipole and gauge. A value::ComplexF64 is returned.
    """
    function amplitude(kind::String, channel::PhotoIonization.Channel, energy::Float64, continuumLevel::Level, initialLevel::Level, grid::Radial.Grid)
        if      kind in [ "photoionization"]
        #-----------------------------------
            amplitude = PhotoEmission.amplitude("absorption", channel.multipole, channel.gauge, energy, continuumLevel, initialLevel, grid,
                                                display=false, printout=false)
            amplitude = im^Basics.subshell_l(Subshell(101, channel.kappa)) * exp( -im*channel.phase ) * amplitude
        else    error("stop b")
        end

        return( amplitude )
    end

    """
    `PhotoIonization.computeAmplitudesProperties(line::PhotoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                                     settings::PhotoIonization.Settings; printout::Bool=false)`
        ... to compute all amplitudes and properties of the given line; a line::PhotoIonization.Line is returned for which the amplitudes and
            properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::PhotoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                          settings::PhotoIonization.Settings; printout::Bool=false)
        newChannels = PhotoIonization.Channel[];;   contSettings = Continuum.Settings(false, nrContinuum);    csC = 0.;    csB = 0.
        for channel in line.channels
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel, line.initialLevel.basis.subshells)
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel, newiLevel.basis.subshells)
            newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newiLevel)
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newfLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
            newChannel = PhotoIonization.Channel(channel.multipole, channel.gauge, channel.kappa, channel.symmetry, phase, 0.)
            amplitude  = PhotoIonization.amplitude("photoionization", newChannel, line.photonEnergy, newcLevel, newiLevel, grid)
            push!( newChannels, PhotoIonization.Channel(newChannel.multipole, newChannel.gauge, newChannel.kappa, newChannel.symmetry,
                                                        newChannel.phase, amplitude) )
            if       channel.gauge == Basics.Coulomb     csC = csC + abs(amplitude)^2
            elseif   channel.gauge == Basics.Babushkin   csB = csB + abs(amplitude)^2
            elseif   channel.gauge == Basics.Magnetic    csB = csB + abs(amplitude)^2;   csC = csC + abs(amplitude)^2
            end
        end
        Ji2 = Basics.twice(line.initialLevel.J)
        ##x csFactor     = 4 * pi^2 * Defaults.getDefaults("alpha") * line.photonEnergy / (2*(Ji2 + 1))
        ##x csFactor     = 4 * pi^2 * Defaults.getDefaults("alpha") / line.photonEnergy / (Ji2 + 1)
        ##x csFactor     = 4 * pi^2 / Defaults.getDefaults("alpha") / line.photonEnergy / (Ji2 + 1)
        csFactor     = 8 * pi^3 / Defaults.getDefaults("alpha") / line.photonEnergy
        crossSection = EmProperty(csFactor * csC, csFactor * csB)
        newLine = PhotoIonization.Line( line.initialLevel, line.finalLevel, line.electronEnergy, line.photonEnergy,
                                        crossSection, true, newChannels)
        return( newLine )
    end

    """
    `PhotoIonization.computeAmplitudesPropertiesPlasma(line::PhotoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid,
                                                           settings::PlasmaShift.PhotoSettings)`
        ... to compute all amplitudes and properties of the given line but for the given plasma model;
            a line::PhotoIonization.Line is returned for which the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesPropertiesPlasma(line::PhotoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, settings::PlasmaShift.PhotoSettings)
        newChannels = PhotoIonization.Channel[];;   contSettings = Continuum.Settings(false, grid.NoPoints-50);    csC = 0.;    csB = 0.
        for channel in line.channels
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel)
            newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newiLevel)
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel)
            @warn "Adapt a proper continuum orbital for the plasma potential"
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newfLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
            newChannel = PhotoIonization.Channel(channel.multipole, channel.gauge, channel.kappa, channel.symmetry, phase, 0.)
            @warn "Adapt a proper Auger amplitude for the plasma e-e interaction"
            amplitude = 1.0
            # amplitude  = PhotoIonization.amplitude("photoionization", channel, line.photonEnergy, newcLevel, newiLevel, grid)
            push!( newChannels, PhotoIonization.Channel(newChannel.multipole, newChannel.gauge, newChannel.kappa, newChannel.symmetry,
                                                        newChannel.phase, amplitude) )
            if       channel.gauge == Basics.Coulomb     csC = csC + abs(amplitude)^2
            elseif   channel.gauge == Basics.Babushkin   csB = csB + abs(amplitude)^2
            elseif   channel.gauge == Basics.Magnetic    csB = csB + abs(amplitude)^2;   csC = csC + abs(amplitude)^2
            end
        end
        Ji2 = Basics.twice(line.initialLevel.J)
        csFactor     = 4 * pi^2 * Defaults.getDefaults("alpha") * line.photonEnergy / (2*(Ji2 + 1))
        crossSection = EmProperty(csFactor * csC, csFactor * csB)
        println("plasma-photo cs = $crossSection")
        newline = PhotoIonization.Line( line.initialLevel, line.finalLevel, line.electronEnergy, line.photonEnergy, crossSection, true, newChannels)
        return( newline )
    end

    """
    `PhotoIonization.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                      settings::PhotoIonization.Settings; output::Bool=true)`
        ... to compute the photoIonization transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{PhotoIonization.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                           settings::PhotoIonization.Settings; output::Bool=true)
        println("")
        printstyled("PhotoIonization.computeLines(): The computation of photo-ionization and properties starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        lines = PhotoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    PhotoIonization.displayLines(lines)    end
        # Determine maximum energy and check for consistency of the grid
        maxEnergy = 0.;   for  line in lines   maxEnergy = max(maxEnergy, line.electronEnergy)   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newLines = PhotoIonization.Line[]
        for  line in lines
            println("\n>> Calculate photoionization amplitudes and properties for line: $(line.initialLevel.index) - $(line.finalLevel.index) " *
                    "for the photon energy $(Defaults.convertUnits("energy: from atomic", line.photonEnergy)) " * Defaults.GBL_ENERGY_UNIT)
            newLine = PhotoIonization.computeAmplitudesProperties(line, nm, grid, nrContinuum, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        PhotoIonization.displayResults(stdout, newLines, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoIonization.displayResults(iostream, newLines, settings)     end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `PhotoIonization.computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                         settings::PhotoIonization.Settings; output=true, printout::Bool=true)`
        ... to compute the photoionization transition amplitudes and all properties as requested by the given settings. The computations
            and printout is adapted for large cascade computations by including only lines with at least one channel and by sending
            all printout to a summary file only. A list of lines::Array{PhotoIonization.Lines} is returned.
    """
    function  computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                  settings::PhotoIonization.Settings; output=true, printout::Bool=true)

        lines = PhotoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        # if  settings.printBefore    PhotoIonization.displayLines(lines)    end
        # Determine maximum energy and check for consistency of the grid
        maxEnergy = 0.;   for  en in settings.photonEnergies   maxEnergy = max(maxEnergy, Defaults.convertUnits("energy: to atomic", en))   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newLines = PhotoIonization.Line[]
        for  (i,line)  in  enumerate(lines)
            if  rem(i,10) == 0    println("> Photo line $i:  ... calculated ")    end
            newLine = PhotoIonization.computeAmplitudesProperties(line, nm, grid, nrContinuum, settings, printout=printout)
            ##x if  rem(i,10) == 0    println("> Photo line $i:  ... not calculated ")    end
            ##x newLine = PhotoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, line.photonEnergy, Basics.EmProperty(1.0, 0.),
            ##x                                false, PhotoIonization.Channel[] )
            push!( newLines, newLine)
        end
        # Print all results to a summary file, if requested
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoIonization.displayResults(iostream, newLines, settings)   end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `PhotoIonization.computeLinesPlasma(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                            settings::PlasmaShift.PhotoSettings; output::Bool=true)`
        ... to compute the photoIonization transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{PhotoIonization.Lines} is returned.
    """
    function  computeLinesPlasma(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                 settings::PlasmaShift.PhotoSettings; output::Bool=true)
        println("")
        printstyled("PhotoIonization.computeLinesPlasma(): The computation of photo-ionization cross sections starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        photoSettings = PhotoIonization.Settings(settings.multipoles, settings.gauges, settings.photonEnergies, false, false, false,
                                                 settings.printBefore, settings.selectLines, settings.selectedLines, Basics.ExpStokes() )

        lines = PhotoIonization.determineLines(finalMultiplet, initialMultiplet, photoSettings)
        # Display all selected lines before the computations start
        if  settings.printBefore    PhotoIonization.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = PhotoIonization.Line[]
        for  line in lines
            newLine = PhotoIonization.computeAmplitudesPropertiesPlasma(line, nm, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        PhotoIonization.displayResults(stdout, lines, photoSettings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoIonization.displayResults(iostream, lines, photoSettings)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `PhotoIonization.computePartialCrossSectionUnpolarized(gauge::EmGauge, Mf::AngularM64, line::PhotoIonization.Line)`
        ... to compute the partial photoionization cross section for initially unpolarized atoms by unpolarized plane-wave photons.
            A value::Float64 is returned.
    """
    function  computePartialCrossSectionUnpolarized(gauge::EmGauge, Mf::AngularM64, line::PhotoIonization.Line)
        # Define an internal Racah expression for the summation over t, lambda
        function Racahexpr(kappa::Int64, Ji::AngularJ64, Jf::AngularJ64, Mf::AngularM64, J::AngularJ64, Jp::AngularJ64,
                           L::Int64, Lp::Int64, p::Int64, pp::Int64)
            # Determine the allowed values of t
            t1 = Basics.oplus( AngularJ64(Lp), Jf);    t2 = Basics.oplus( AngularJ64(L), Jf);    tList = intersect(t1, t2)
            wb = 0.
            for  t  in tList
                for  lambda = -1:2:1
                    j = AngularMomentum.kappa_j(kappa);   Mf_lambda = Basics.add(AngularM64(lambda), Mf)
                    wb = wb + (1.0im * lambda)^p * (-1.0im * lambda)^pp *
                         AngularMomentum.ClebschGordan( AngularJ64(Lp), AngularM64(lambda), Jf, Mf, t, Mf_lambda) *
                         AngularMomentum.ClebschGordan( AngularJ64(L),  AngularM64(lambda), Jf, Mf, t, Mf_lambda) *
                         AngularMomentum.Wigner_9j(j, Jp, Jf, J, Ji, AngularJ64(L), Jf, AngularJ64(Lp), t)
                end
            end

            return( wb )
        end

        if  !line.hasChannels   error("No channels are defined for the given PhotoIonization.line.")         end
        wa = 0.0im;    Ji = line.initialLevel.J;    Jf = line.finalLevel.J;
        kappaList = PhotoIonization.getLineKappas(line)
        for  kappa in kappaList
            for  cha  in line.channels
                if  kappa != cha.kappa  ||  (gauge != cha.gauge  &&  gauge != Basics.Magnetic)   continue    end
                J = cha.symmetry.J;    L = cha.multipole.L;    if  cha.multipole.electric   p = 1   else    p = 0   end
                for  chp  in line.channels
                    if  kappa != chp.kappa  ||  (gauge != cha.gauge  &&  gauge != Basics.Magnetic)    continue    end
                    Jp = chp.symmetry.J;    Lp = chp.multipole.L;    if  chp.multipole.electric   pp = 1   else    pp = 0   end
                    wa = wa + 1.0im^(L - Lp) * (-1)^(L + Lp) * AngularMomentum.bracket([AngularJ64(L), AngularJ64(Lp), J, Jp]) *
                         Racahexpr(kappa, Ji, Jf, Mf, J, Jp, L, Lp, p, pp) * cha.amplitude * conj(chp.amplitude)
                end
            end
        end
        csFactor = 8 * pi^3 * Defaults.getDefaults("alpha") / (2*line.photonEnergy * (Basics.twice(Ji) + 1))
        wa       = csFactor * wa

        return( wa )
    end

    """
    `PhotoIonization.computeStatisticalTensorUnpolarized(k::Int64, q::Int64, gauge::EmGauge, line::PhotoIonization.Line,
                                                             settings::PhotoIonization.Settings)`
        ... to compute the statistical tensor of the photoion in its final level after the photoionization of initially unpolarized atoms
            by plane-wave photons with given Stokes parameters (density matrix). A value::ComplexF64 is returned.
    """
    function  computeStatisticalTensorUnpolarized(k::Int64, q::Int64, gauge::EmGauge, line::PhotoIonization.Line,
                                                  settings::PhotoIonization.Settings)
        if  !line.hasChannels   error("No channels are defined for the given PhotoIonization.line.")         end
        wa = 0.0im;    Ji = line.initialLevel.J;    Jf = line.finalLevel.J
        kappaList = PhotoIonization.getLineKappas(line);    P1 = settings.stokes.P1;   P2 = settings.stokes.P2;   P3 = settings.stokes.P3
        for  kappa in kappaList
            j = AngularMomentum.kappa_j(kappa)
            for  cha  in line.channels
                if  kappa != cha.kappa  ||  (gauge != cha.gauge  &&  gauge != Basics.Magnetic)   continue    end
                J = cha.symmetry.J;    L = cha.multipole.L;    if  cha.multipole.electric   p = 1   else    p = 0   end
                for  chp  in line.channels
                    if  kappa != chp.kappa  ||  (gauge != cha.gauge  &&  gauge != Basics.Magnetic)    continue    end
                    Jp = chp.symmetry.J;    Lp = chp.multipole.L;    if  chp.multipole.electric   pp = 1   else    pp = 0   end
                    #
                    for  lambda = -1:2:1
                        for  lambdap = -1:2:1
                            if  lambda == lambdap   wb = (1.0 + 0.0im +lambda*P3)    else    wb = P1 - lambda * P2 * im    end
                            wa = wa + wb * 1.0im^(L - Lp + p - pp) * lambda^p * lambdap^pp *
                                 sqrt( AngularMomentum.bracket([AngularJ64(L), AngularJ64(Lp), J, Jp]) ) *
                                 AngularMomentum.phaseFactor([J, +1, Jp, +1, Jf, +1, Ji, +1, j, +1, AngularJ64(1)]) *
                                 AngularMomentum.ClebschGordan( AngularJ64(L),  AngularM64(lambda), AngularJ64(Lp),  AngularM64(-lambda),
                                                                    AngularJ64(k),  AngularM64(q)) *
                                 AngularMomentum.Wigner_6j(Jf, j, Jp, J, AngularJ64(k), Jf) *
                                 AngularMomentum.Wigner_6j(Jp, Ji, AngularJ64(Lp), AngularJ64(L), AngularJ64(k), J) *
                                 cha.amplitude * conj(chp.amplitude)
                        end
                    end
                end
            end
        end

        wa = pi / (Basics.twice(Ji) + 1) * wa
        return( wa )
    end

    """
    `PhotoIonization.determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoIonization.Settings)`
        ... to determine a list of photoionization Channel for a transitions from the initial to final level and by taking into account
            the particular settings of for this computation; an Array{PhotoIonization.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoIonization.Settings)
        channels = PhotoIonization.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        if  Basics.UseCoulomb  in  settings.gauges   gaugeM = Basics.UseCoulomb    else   gaugeM = Basics.UseBabushkin    end
        for  mp in settings.multipoles
            for  gauge in settings.gauges
                symList = AngularMomentum.allowedMultipoleSymmetries(symi, mp)
                ##x println("mp = $mp   symi = $symi   symList = $symList")
                for  symt in symList
                    kappaList = AngularMomentum.allowedKappaSymmetries(symt, symf)
                    for  kappa in kappaList
                        # Include further restrictions if appropriate
                        if     string(mp)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                            push!(channels, PhotoIonization.Channel(mp, Basics.Coulomb,   kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                            push!(channels, PhotoIonization.Channel(mp, Basics.Babushkin, kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'M'  &&   gauge == gaugeM
                            push!(channels, PhotoIonization.Channel(mp, Basics.Magnetic,  kappa, symt, 0., Complex(0.)) )
                        end
                    end
                end
            end
        end
        return( channels )
    end

    """
    `PhotoIonization.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoIonization.Settings)`
        ... to determine a list of PhotoIonization.Line's for transitions between levels from the initial- and final-state multiplets,
            and  by taking into account the particular selections and settings for this computation; an Array{PhotoIonization.Line,1}
            is returned. Apart from the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoIonization.Settings)
        lines = PhotoIonization.Line[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    for  omega in settings.photonEnergies
                        # Photon energies are still in 'pre-defined' units; convert to Hartree
                        omega_au = Defaults.convertUnits("energy: to atomic", omega)
                        energy   = omega_au - (fLevel.energy - iLevel.energy)
                        if  energy < 0    continue   end
                        channels = PhotoIonization.determineChannels(fLevel, iLevel, settings)
                        push!( lines, PhotoIonization.Line(iLevel, fLevel, energy, omega_au, EmProperty(0., 0.), true, channels) )
                    end
                end
            end
        end
        return( lines )
    end

    """
    `PhotoIonization.displayLines(lines::Array{PhotoIonization.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all selected
            transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{PhotoIonization.Line,1})
        nx = 175
        println(" ")
        println("  Selected photoionization lines:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"   ; na=2);                       sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(10, "Energy_fi"; na=3);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(10, "omega"; na=3);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(12, "Energy e_p"; na=3);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(57, "List of multipoles, gauges, kappas and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(57, "partial (multipole, gauge, total J^P)                  "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        nchannels = 0
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=3)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))              * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
            kappaMultipoleSymmetryList = Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry}[]
            for  i in 1:length(line.channels)
                push!( kappaMultipoleSymmetryList, (line.channels[i].kappa, line.channels[i].multipole, line.channels[i].gauge,
                                                    line.channels[i].symmetry) )
                nchannels = nchannels + 1
            end
            ##x println("PhotoIonization-diplayLines-ad: kappaMultipoleSymmetryList = ", kappaMultipoleSymmetryList)
            wa = TableStrings.kappaMultipoleSymmetryTupels(85, kappaMultipoleSymmetryList)
            sb = sa * wa[1];    println( sb )
            for  i = 2:length(wa)
                sb = TableStrings.hBlank( length(sa) ) * wa[i];    println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx), "\n")
        println("  A total of $nchannels channels need to be calculated. \n")
        #
        return( nothing )
    end

    """
    `PhotoIonization.displayResults(stream::IO, lines::Array{PhotoIonization.Line,1}, settings::PhotoIonization.Settings)`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayResults(stream::IO, lines::Array{PhotoIonization.Line,1}, settings::PhotoIonization.Settings)
        nx = 130
        println(stream, " ")
        println(stream, "  Total photoionization cross sections for initially unpolarized atoms by unpolarized plane-wave photons:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"   ; na=2);                       sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "f--Energy--i"; na=4)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(12, "omega"     ; na=4)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(12, "Energy e_p"; na=3)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(10, "Multipoles"; na=1);                              sb = sb * TableStrings.hBlank(13)
        sa = sa * TableStrings.center(30, "Cou -- Cross section -- Bab"; na=3)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section") * "          " *
                                              TableStrings.inUnits("cross section"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=3)
            en = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
            multipoles = EmMultipole[]
            for  ch in line.channels
                multipoles = push!( multipoles, ch.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=2)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Coulomb))     * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Babushkin))   * "                 "
            ##x sa = sa * @sprintf("%.6e", line.crossSection.Coulomb)     * "    "
            ##x sa = sa * @sprintf("%.6e", line.crossSection.Babushkin)   * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        #
        if  settings.calcPartialCs
            nx = 144
            println(stream, " ")
            println(stream, "  Partial cross sections for initially unpolarized atoms by unpolarized plane-wave photons:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
                sa = "  ";   sb = "  "
            sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(18)
            sa = sa * TableStrings.center(18, "i--J^P--f"   ; na=2);                       sb = sb * TableStrings.hBlank(22)
            sa = sa * TableStrings.center(12, "f--Energy--i"; na=4)
            sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "omega"     ; na=4)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "Energy e_p"; na=3)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
            sa = sa * TableStrings.center(10, "Multipoles"; na=1);                              sb = sb * TableStrings.hBlank(13)
            sa = sa * TableStrings.center( 7, "M_f"; na=1);                                     sb = sb * TableStrings.hBlank(11)
            sa = sa * TableStrings.center(30, "Cou -- Partial cross section -- Bab"; na=3)
            sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section") * "          " *
                                              TableStrings.inUnits("cross section"); na=3)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
            #
            for  line in lines
                sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                             fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
                sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
                sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=3)
                en = line.finalLevel.energy - line.initialLevel.energy
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                  * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
                multipoles = EmMultipole[]
                for  ch in line.channels
                    multipoles = push!( multipoles, ch.multipole)
                end
                multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
                sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=2)
                println(stream, sa)
                MfList = Basics.projections(line.finalLevel.J)
                for  Mf in MfList
                    sb  = TableStrings.hBlank(97)
                    wac = PhotoIonization.computePartialCrossSectionUnpolarized(Basics.Coulomb, Mf, line)
                    wab = PhotoIonization.computePartialCrossSectionUnpolarized(Basics.Babushkin, Mf, line)
                    sb  = sb * TableStrings.flushright( 8, string(Mf))                             * "       "
                    sb  = sb * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", wac.re))     * "    "
                    sb  = sb * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", wab.re))     * "    "
                    println(stream, sb)
                end
            end
            println(stream, "  ", TableStrings.hLine(nx))
        end
        #
        #
        if  settings.calcTensors
            nx = 144
            println(stream, " ")
            println(stream, "  Reduced statistical tensors of the photoion in its final level after the photoionization ")
            println(stream, "  of initially unpolarized atoms by plane-wave photons with given Stokes parameters (density matrix):")
            println(stream, "\n     + tensors are printed for k = 0, 1, 2 and if non-zero only.")
            println(stream,   "     + Stokes parameters are:  P1 = $(settings.stokes.P1),  P2 = $(settings.stokes.P2),  P3 = $(settings.stokes.P3) ")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
                sa = "  ";   sb = "  "
            sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(18)
            sa = sa * TableStrings.center(18, "i--J^P--f"   ; na=2);                       sb = sb * TableStrings.hBlank(22)
            sa = sa * TableStrings.center(12, "f--Energy--i"; na=4)
            sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "omega"     ; na=4)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "Energy e_p"; na=3)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
            sa = sa * TableStrings.center(10, "Multipoles"; na=1);                              sb = sb * TableStrings.hBlank(13)
            sa = sa * TableStrings.center(10, "k    q"; na=4);                                  sb = sb * TableStrings.hBlank(11)
            sa = sa * TableStrings.center(30, "Cou --  rho_kq (J_f)  -- Bab"; na=3)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
            #
            for  line in lines
                sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                             fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
                sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
                sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=3)
                en = line.finalLevel.energy - line.initialLevel.energy
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                  * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
                multipoles = EmMultipole[]
                for  ch in line.channels
                    multipoles = push!( multipoles, ch.multipole)
                end
                multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
                sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=2)
                println(stream, sa)
                for  k = 0:2
                    for q = -k:k
                        sb   = TableStrings.hBlank(102)
                        rhoc = PhotoIonization.computeStatisticalTensorUnpolarized(k, q, Basics.Coulomb,   line, settings)
                        rhob = PhotoIonization.computeStatisticalTensorUnpolarized(k, q, Basics.Babushkin, line, settings)
                        sb   = sb * string(k) * " " * TableStrings.flushright( 4, string(q))             * "       "
                        sb   = sb * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", rhoc.re))     * "    "
                        sb   = sb * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", rhob.re))     * "    "
                        println(stream, sb)
                    end
                end
            end
            println(stream, "  ", TableStrings.hLine(nx))
        end
        #
        return( nothing )
    end

    """
    `PhotoIonization.getLineKappas(line::PhotoIonization.Line)`
        ... returns a list of kappa-values (partial waves) which contribute to the given line, to which one or several channels are
            assigned. An kappaList::Array{Int64,1} is returned.
    """
    function getLineKappas(line::PhotoIonization.Line)
        kappaList = Int64[]
        for  ch  in line.channels
            kappaList = union(kappaList, ch.kappa)
        end

        return( kappaList )
    end

end # module


text:
# Fallback constructor for diffusivity types without precomputed diffusivities.
DiffusivityFields(arch::AbstractArchitecture, grid::AbstractGrid, args...) = nothing

# For a tuple of closures
DiffusivityFields(
    arch::AbstractArchitecture, grid::AbstractGrid, tracers,
    bcs::NamedTuple, closure_tuple::Tuple
) = Tuple(DiffusivityFields(arch, grid, tracers, bcs, closure) for closure in closure_tuple)

#####
##### For closures that only require an eddy viscosity νₑ field.
#####

const ViscosityClosures = Union{AbstractSmagorinsky, AbstractLeith}

DiffusivityFields(
    arch::AbstractArchitecture, grid::AbstractGrid, tracers, ::ViscosityClosures;
    νₑ = CellField(arch, grid, DiffusivityBoundaryConditions(grid), zeros(arch, grid))
) = (νₑ=νₑ,)

function DiffusivityFields(arch::AbstractArchitecture, grid::AbstractGrid, tracers,
                           bcs::NamedTuple, ::ViscosityClosures)
    νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid)
    νₑ = CellField(arch, grid, νₑ_bcs, zeros(arch, grid))
    return (νₑ=νₑ,)
end

#####
##### For closures that also require tracer diffusivity fields κₑ on each tracer.
#####

const ViscosityDiffusivityClosures = Union{VAMD, RAMD}

function TracerDiffusivityFields(arch, grid, tracer_names; kwargs...)
    κ_fields =
        Tuple(c ∈ keys(kwargs) ?
              kwargs[c] :
              CellField(arch, grid, DiffusivityBoundaryConditions(grid), zeros(arch, grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(κ_fields)
end

function TracerDiffusivityFields(arch, grid, tracer_names, bcs::NamedTuple)
    κ_fields =
        Tuple(c ∈ keys(bcs) ?
              CellField(arch, grid, bcs[c],                              zeros(arch, grid)) :
              CellField(arch, grid, DiffusivityBoundaryConditions(grid), zeros(arch, grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(κ_fields)
end

function DiffusivityFields(
    arch::AbstractArchitecture, grid::AbstractGrid, tracers, ::ViscosityDiffusivityClosures;
    νₑ = CellField(arch, grid, DiffusivityBoundaryConditions(grid), zeros(arch, grid)), kwargs...)
    κₑ = TracerDiffusivityFields(arch, grid, tracers; kwargs...)
    return (νₑ=νₑ, κₑ=κₑ)
end

function DiffusivityFields(arch::AbstractArchitecture, grid::AbstractGrid, tracers,
                           bcs::NamedTuple, ::ViscosityDiffusivityClosures)

    νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid)
    νₑ = CellField(arch, grid, νₑ_bcs, zeros(arch, grid))

    κₑ = :κₑ ∈ keys(bcs) ?
        TracerDiffusivityFields(arch, grid, tracers, bcs[:κₑ]) :
        TracerDiffusivityFields(arch, grid, tracers)

    return (νₑ=νₑ, κₑ=κₑ)
end


text:
using XGBoost

# load file from text file, also binary buffer generated by xgboost

dtrain = DMatrix("../data/agaricus.txt.train")
dtest = DMatrix("../data/agaricus.txt.test")

#Defining parameters for xgboost

param = ["max_depth" => 2,
         "eta" => 1,
         "silent" => 1,
         "objective" => "binary:logistic"]
num_round = 2
nfold = 5

print ("running cross validation\n")
# do cross validation, this will print result out as
# [iteration]  metric_name:mean_value+std_value
# std_value is standard deviation of the metric
nfold_cv(dtrain, num_round, nfold, param = param, metrics = ["error"], seed = 0)

print ("running cross validation, disable standard deviation display\n")
# do cross validation, this will print result out as
# [iteration]  metric_name:mean_value+std_value
# std_value is standard deviation of the metric
nfold_cv(dtrain, num_round, nfold, param = param, metrics = ["error"], seed = 0, show_stdv = false)

print ("running cross validation, with preprocessing function\n")
# define the preprocessing function
# used to return the preprocessed training, test data, and parameter
# we can use this to do weight rescale, etc.
# as a example, we try to set scale_pos_weight
function fpreproc(dtrain::DMatrix, dtest::DMatrix, param)
    label = get_info(dtrain, "label")
    ratio = sum(label == 0) / sum(label == 1)
    param["scale_pos_weight"] = ratio
    return dtrain, dtest, param
end
# do cross validation, for each fold
# the dtrain, dtest, param will be passed into fpreproc
# then the return value of fpreproc will be used to generate
# results of that fold
nfold_cv(dtrain, num_round, nfold, param = param, metrics = ["auc"],
         seed = 0, show_stdv = false, fpreproc = fpreproc)

print("running cross validation, with customized loss function\n")
###
# you can also do cross validation with cutomized loss function
# See custom_objective.py
##

function logregobj(preds::Vector{Float32}, dtrain::DMatrix)
        labels = get_info(dtrain, "label")
        preds = 1.0 ./ (1.0 + exp(-preds))
        grad = preds - labels
        hess = preds .* (1. - preds)
        return grad, hess
end

function evalerror(preds::Vector{Float32}, dtrain::DMatrix)
    labels = get_info(dtrain, "label")
    # return a pair metric_name, result
    # since preds are margin(before logistic transformation, cutoff at 0)
    return "self-error", sum((preds .> 0.0) .!= labels) / float(size(preds)[1])
end

# train with customized objective
nfold_cv(dtrain, num_round, nfold, metrics = [], seed = 0, obj = logregobj, feval = evalerror,
         max_depth = 2, eta = 1, silent = 1)


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

"""

    MultiSelectMenu

A menu that allows a user to select a multiple options from a list.

# Sample Output

```julia
julia> request(MultiSelectMenu(options))
Select the fruits you like:
[press: d=done, a=all, n=none]
   [ ] apple
 > [X] orange
   [X] grape
   [ ] strawberry
   [ ] blueberry
   [X] peach
   [ ] lemon
   [ ] lime
You like the following fruits:
  - orange
  - grape
  - peach
```

"""
mutable struct MultiSelectMenu{C} <: _ConfiguredMenu{C}
    options::Array{String,1}
    pagesize::Int
    pageoffset::Int
    selected::Set{Int}
    config::C
end

"""

    MultiSelectMenu(options::Array{String,1}; pagesize::Int=10, selected=[], kwargs...)

Create a MultiSelectMenu object. Use `request(menu::MultiSelectMenu)` to get
user input. It returns a `Set` containing the indices of options that
were selected by the user.

# Arguments

  - `options::Array{String, 1}`: Options to be displayed
  - `pagesize::Int=10`: The number of options to be displayed at one time, the menu will scroll if length(options) > pagesize
  - `selected=[]`: pre-selected items. `i ∈ selected` means that `options[i]` is preselected.

Any additional keyword arguments will be passed to [`TerminalMenus.MultiSelectConfig`](@ref).
"""
function MultiSelectMenu(options::Array{String,1}; pagesize::Int=10, selected=Int[], warn::Bool=true, kwargs...)
    length(options) < 2 && error("MultiSelectMenu must have at least two options")

    # if pagesize is -1, use automatic paging
    pagesize = pagesize == -1 ? length(options) : pagesize
    # pagesize shouldn't be bigger than options
    pagesize = min(length(options), pagesize)
    # after other checks, pagesize must be greater than 2
    pagesize < 2 && error("pagesize must be >= 2")

    pageoffset = 0
    _selected = Set{Int}()
    for item in selected
        push!(_selected, item)
    end

    if !isempty(kwargs)
        MultiSelectMenu(options, pagesize, pageoffset, _selected, MultiSelectConfig(; kwargs...))
    else
        warn && Base.depwarn("Legacy `MultiSelectMenu` interface is deprecated, set a configuration option such as `MultiSelectMenu(options; charset=:ascii)` to trigger the new interface.", :MultiSelectMenu)
        MultiSelectMenu(options, pagesize, pageoffset, _selected, CONFIG)
    end

end

# AbstractMenu implementation functions
# See AbstractMenu.jl
#######################################

header(m::MultiSelectMenu) = "[press: d=done, a=all, n=none]"

options(m::MultiSelectMenu) = m.options

cancel(m::MultiSelectMenu) = m.selected = Set{Int}()

# Do not exit menu when a user selects one of the options
function pick(menu::MultiSelectMenu, cursor::Int)
    if cursor in menu.selected
        delete!(menu.selected, cursor)
    else
        push!(menu.selected, cursor)
    end

    return false #break out of the menu
end

function writeline(buf::IOBuffer, menu::MultiSelectMenu{MultiSelectConfig}, idx::Int, iscursor::Bool)
    if idx in menu.selected
        print(buf, menu.config.checked, " ")
    else
        print(buf, menu.config.unchecked, " ")
    end

    print(buf, replace(menu.options[idx], "\n" => "\\n"))
end

# d: Done, return from request
# a: Select all
# n: Deselect all
function keypress(menu::MultiSelectMenu, key::UInt32)
    if key == UInt32('d') || key == UInt32('D')
        return true # break
    elseif key == UInt32('a') || key == UInt32('A')
        menu.selected = Set(1:length(menu.options))
    elseif key == UInt32('n') || key == UInt32('N')
        menu.selected = Set{Int}()
    end
    false # don't break
end

## Legacy interface
function TerminalMenus.writeLine(buf::IOBuffer, menu::MultiSelectMenu{<:Dict}, idx::Int, cursor::Bool)
    # print a ">" on the selected entry
    cursor ? print(buf, menu.config[:cursor]," ") : print(buf, "  ")
    if idx in menu.selected
        print(buf, menu.config[:checked], " ")
    else
        print(buf, menu.config[:unchecked], " ")
    end

    print(buf, replace(menu.options[idx], "\n" => "\\n"))
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
"""
    struct ElectricPotential{T, N, S, AT} <: AbstractArray{T, N}

Electric potential of the simulation in units of volt (V).

## Parametric types
* `T`: Element type of `data`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{T, N}`: Array containing the values of the electric potential at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points for which the electric potential is determined.
"""
struct ElectricPotential{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S, AT}
end

@inline size(epot::ElectricPotential{T, N, S}) where {T, N, S} = size(epot.data)
@inline length(epot::ElectricPotential{T, N, S}) where {T, N, S} = length(epot.data)
@inline getindex(epot::ElectricPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(epot.data, I...)
@inline getindex(epot::ElectricPotential{T, N, S}, i::Int) where {T, N, S} = getindex(epot.data, i)
@inline getindex(epot::ElectricPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(epot.grid, s)

function NamedTuple(epot::ElectricPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(epot.grid),
        values = epot.data * u"V",
    )
end
Base.convert(T::Type{NamedTuple}, x::ElectricPotential) = T(x)

function ElectricPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    ElectricPotential{T, N, S, typeof(grid.axes)}( ustrip.(uconvert.(u"V", nt.values)), grid)
end
Base.convert(T::Type{ElectricPotential}, x::NamedTuple) = T(x)

"""
    get_ticks_at_positions_of_large_gradient(epot::ElectricPotential)

The electric potential is analyzed in order to find and return ticks where the gradient (electric field) is strong (relativ to its maximum).

In the 1D case of a pn-junction, the electric field strength is the largest at the position of the pn-junction.
Thus, this function is likely to return ticks which are located close to the pn-junction of a semiconductor.
"""
function get_ticks_at_positions_of_large_gradient(epot::ElectricPotential{T, 3}) where T
    dims = (1, 2, 3)
    mps::NTuple{3, Vector{T}} = broadcast(idim -> StatsBase.midpoints(epot.grid[idim]), dims)
    extended = length.(mps) .> 0
    max_diffs = broadcast(idim -> extended[idim] ? map(i->maximum(abs.(selectdim(epot.data, idim, i+1) .- selectdim(epot.data, idim, i))), eachindex(mps[idim]))::Vector{T} : T[], dims);
    max_diff = broadcast(idim -> extended[idim] ? maximum(max_diffs[idim]) : T(0), dims)
    normalized_max_diffs = broadcast(idim -> extended[idim] ? max_diffs[idim] ./ max_diff[idim] : T[], dims)
    inds = broadcast(idim -> extended[idim] ? findall(Δ -> Δ > 0.5, normalized_max_diffs[idim]) : Int[], dims) # 0.5 seems to be a good limit as we don't want to add too many ticks.
    return broadcast(idim -> extended[idim] ? mps[idim][inds[idim]] : T[], dims)
end


text:
function formatinfo(s::Spec, x::AbstractChar)
    c = Char(x)
    size = ncodeunits(c)
    return size + paddingsize(s, 1), c
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, ::AbstractChar, c::Char)
    align = default(s.align, ALIGN_LEFT)
    pw = paddingwidth(s, 1)
    p = padleft(data, p, s.fill, align, pw)
    p = char(data, p, c)
    p = padright(data, p, s.fill, align, pw)
    return p
end


text:
struct MountainCar{S, A, R} <: DiscEnv{S, A, R}
    current::Vector{S}
    observationspace::ObservationSpace{Vector{S}}
    actionspace::ActionSpace{A}
    function MountainCar{T}() where T <: AbstractFloat
        observationspace = ObservationSpace{Vector{T}}(2)
        actionspace = ActionSpace(collect(T, -1:1))
        new{T, T, T}(randpos(), observationspace, actionspace)
    end
end

MountainCar() = MountainCar{Float32}()

struct MountainCarContinuous{S, A, R} <: ContEnv{S, A, R}
    current::Vector{S}
    observationspace::ObservationSpace{Vector{S}}
    actionspace::ActionSpace{Tuple{A, A}}
    function MountainCarContinuous{T}() where T <: AbstractFloat
        observationspace = ObservationSpace{Vector{T}}(2)
        actionspace = ActionSpace([(-one(T), one(T))])
        new{T, T, T}(randpos(), observationspace, actionspace)
    end
end

MountainCarContinuous() = MountainCarContinuous{Float32}()

randpos() = [rand(Uniform(-0.6, -0.4)), 0.0]

step!(env::MountainCar, action::Integer) = step!(env, env.actionspace.actions[action])
function step!(env::MountainCar{S, A, R}, action::A) where {S, A, R}
    p, v = env.current
    v += action * 0.001 - 0.0025 * cos(3.0*p)
    v = clamp(v, -0.07, 0.07)
    p = clamp(p + v, -1.2, 0.6)
    if p == -1.2 && v < 0.0
        v = 0.0
    end
    d = p >= 0.5
    r = d ? zero(R) : -one(R)
    env.current .= p, v
    return env.current, r, d
end

step!(env::MountainCarContinuous{S, A, R}, action::Vector{A}) where {S, A, R} = step!(env, action[1])
function step!(env::MountainCarContinuous{S, A, R}, action::A) where {S, A, R}
    p, v = env.current
    v += clamp(action, -1.0, 1.0) * 0.0015 - 0.0025 * cos(3.0*p)
    v = clamp(v, -0.07, 0.07)
    p = clamp(p + v, -1.2, 0.6)
    if p == -1.2 && v < 0.0
        v = 0.0
    end
    d = p >= 0.45
    r = (d ? 100.0 : 0.0) - (action ^ 2.0) * 0.1
    env.current .= p, v
    return env.current, R(r), d
end

reset!(env::Union{MountainCar, MountainCarContinuous}) = (env.current .= randpos())


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ISCSIPersistentVolumeSource <: SwaggerModel
    chapAuthDiscovery::Any # spec type: Union{ Nothing, Bool } # spec name: chapAuthDiscovery
    chapAuthSession::Any # spec type: Union{ Nothing, Bool } # spec name: chapAuthSession
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    initiatorName::Any # spec type: Union{ Nothing, String } # spec name: initiatorName
    iqn::Any # spec type: Union{ Nothing, String } # spec name: iqn
    iscsiInterface::Any # spec type: Union{ Nothing, String } # spec name: iscsiInterface
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    portals::Any # spec type: Union{ Nothing, Vector{String} } # spec name: portals
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: secretRef
    targetPortal::Any # spec type: Union{ Nothing, String } # spec name: targetPortal

    function IoK8sApiCoreV1ISCSIPersistentVolumeSource(;chapAuthDiscovery=nothing, chapAuthSession=nothing, fsType=nothing, initiatorName=nothing, iqn=nothing, iscsiInterface=nothing, lun=nothing, portals=nothing, readOnly=nothing, secretRef=nothing, targetPortal=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("chapAuthDiscovery"), chapAuthDiscovery)
        setfield!(o, Symbol("chapAuthDiscovery"), chapAuthDiscovery)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("chapAuthSession"), chapAuthSession)
        setfield!(o, Symbol("chapAuthSession"), chapAuthSession)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("initiatorName"), initiatorName)
        setfield!(o, Symbol("initiatorName"), initiatorName)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("iqn"), iqn)
        setfield!(o, Symbol("iqn"), iqn)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("iscsiInterface"), iscsiInterface)
        setfield!(o, Symbol("iscsiInterface"), iscsiInterface)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("portals"), portals)
        setfield!(o, Symbol("portals"), portals)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("targetPortal"), targetPortal)
        setfield!(o, Symbol("targetPortal"), targetPortal)
        o
    end
end # type IoK8sApiCoreV1ISCSIPersistentVolumeSource

const _property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("chapAuthDiscovery")=>Symbol("chapAuthDiscovery"), Symbol("chapAuthSession")=>Symbol("chapAuthSession"), Symbol("fsType")=>Symbol("fsType"), Symbol("initiatorName")=>Symbol("initiatorName"), Symbol("iqn")=>Symbol("iqn"), Symbol("iscsiInterface")=>Symbol("iscsiInterface"), Symbol("lun")=>Symbol("lun"), Symbol("portals")=>Symbol("portals"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("targetPortal")=>Symbol("targetPortal"))
const _property_types_IoK8sApiCoreV1ISCSIPersistentVolumeSource = Dict{Symbol,String}(Symbol("chapAuthDiscovery")=>"Bool", Symbol("chapAuthSession")=>"Bool", Symbol("fsType")=>"String", Symbol("initiatorName")=>"String", Symbol("iqn")=>"String", Symbol("iscsiInterface")=>"String", Symbol("lun")=>"Int32", Symbol("portals")=>"Vector{String}", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("targetPortal")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ISCSIPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ISCSIPersistentVolumeSource)
    (getproperty(o, Symbol("iqn")) === nothing) && (return false)
    (getproperty(o, Symbol("lun")) === nothing) && (return false)
    (getproperty(o, Symbol("targetPortal")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, name::Symbol, val)
end


text:
@info "loading packages"
using Distributed
addprocs(30; exeflags="--project")

@everywhere begin

using TicketToRide
using Flux
using Flux.Optimise
using TicketToRide: train!
using LuxurySparse
using LinearAlgebra

function task(n, nlayers; verbose=false, nprune=10, nepochs=10, niteration=100, relative_error=true, least_prune=10)
    opt = Optimise.Descent(1e-2)
    circuit = variational_circuit(n, nlayers);
    ham = heisenberg1D(n);
    dispatch!(circuit, :random);
    circuit, history = prune_train(opt, circuit, ham;
        nprune=nprune, nepochs=nepochs,
        niteration=niteration,
        relative_error=relative_error,
        least_prune=least_prune);

    opt = Optimise.Descent(1e-3)
    circuit = variational_circuit(n, nlayers);
    ham = heisenberg1D(n);
    dispatch!(circuit, :random);
    circuit, history = prune_train(opt, circuit, ham;
        nprune=nprune, nepochs=nepochs,
        niteration=niteration,
        relative_error=relative_error,
        least_prune=least_prune);

    return circuit, history
end

function run_task(total, n, nlayers; verbose=false, nprune=10, nepochs=10, niteration=100, relative_error=true, least_prune=10)
    results = pmap(1:total) do k
        task(n, nlayers;
        nprune=nprune, nepochs=nepochs,
        niteration=niteration,
        relative_error=relative_error,
        least_prune=least_prune)
    end

    min_energy = 0
    circuit, history = nothing, nothing
    for each in results
        if each[2][end] < min_energy
            circuit = each[1]
            history = each[2]
        end
    end
    return circuit, history
end

end

using JLD2
# @time c, h = run_task(2, 2, 2; nepochs=2, nprune=5, niteration=2, least_prune=5)

@time c, h = run_task(nprocs(), 10, 2; nepochs=2, nprune=5, niteration=2, least_prune=5)
c, h = run_task(nprocs(), 10, 100; nepochs=20, nprune=50, niteration=2000, least_prune=5)

jldopen("reinitialize-data-10-10.jld", "w+") do f
    f["circuit"] = c
    f["history"] = h
end

# opt = Optimise.Descent(1e-3)
# circuit = variational_circuit(10, 100);
# nparameters(circuit)
# dispatch!(circuit, :random);
# history = train!(opt, 2000, heisenberg1D(10), circuit; verbose=true)
# prune_train(opt, circuit, heisenberg1D(10);
#     nprune=50, nepochs=20,
#     niteration=500,
#     relative_error=true,
#     least_prune=5);

# nparameters(circuit)


text:
# Preamble

const _DEFAULT_PREAMBLE = Ref{Vector{String}}()
_CUSTOM_PREAMBLE = String[]

"""
    preamble()

Returns the current preamble.
"""
function preamble()
    pre = String[]
    append!(pre, _DEFAULT_PREAMBLE[])
    append!(pre, _CUSTOM_PREAMBLE)
    return pre
end

"""
    reset_preamble!()

Resets the preamble to the default.
"""
function reset_preamble!()
    empty!(_CUSTOM_PREAMBLE)
    return preamble()
end

"""
    add_preamble!(pre::AbstractString...; reset=false)

Adds preamble statements on top of the default preamble. `pre` is any number of `AbstractString`s.

If you already have a number of preamble statements and you would like to get rid of them first before adding new statements, pass `reset=true`.
"""
function add_preamble!(pre::AbstractString...; reset::Bool=false)
    reset && reset_preamble!()
    append!(_CUSTOM_PREAMBLE, String.(pre))
    return preamble()
end

# Document

function _assemble_document(latex_content::AbstractString; standalone::Bool=false)
    if standalone
        return String(latex_content)
    else
        doc = """
        \\documentclass[12pt]{article}
        $(join(preamble(), "\n"))
        \\pagestyle{empty}
        \\begin{document}
        $(String(latex_content))
        \\end{document}
        """
        return doc
    end
end


text:
function reflectioncoeff(n₁::T, n₂::T, θᵢ::Real) where T <: Number
    θₜ = snellius(n₁, n₂, θᵢ)
    rₛ = (n₁*cos(θᵢ) - n₂*cos(θₜ)) ./ (n₁*cos(θᵢ) + n₂*cos(θₜ))
    rₚ = (n₂*cos(θᵢ) - n₁*cos(θₜ)) ./ (n₁*cos(θₜ) + n₂*cos(θᵢ))
    rₛ, rₚ
end

function reflection(n₁::Number, n₂::Number, θᵢ::Real)
    rₛ, rₚ = reflectioncoeff(n₁, n₂, θᵢ)
    Rₛ = abs(rₛ).^2
    Rₚ = abs(rₚ).^2
    Rₛ, Rₚ
end

function reflection(m::RefractiveIndex.Material, λ::Real, θ::Real)
    k₂ = m.kfun(λ)
    typeof(k₂) == Void ? (n₂ = m.nfun(λ)) : (n₂ = m.nfun(λ) + k₂*im)
    reflection(one(typeof(n₂)), n₂, θ)
end

function phaseangle(n₁::Number, n₂::Number, θᵢ::Real)
    rₛ, rₚ = reflectioncoeff(n₁, n₂, θᵢ)
    Φₛ = angle(rₛ)
    Φₚ = angle(rₚ)
    Φₛ, Φₚ
end

function phaseangle(m::RefractiveIndex.Material, λ::Real, θ::Real)
    k₂ = m.kfun(λ)
    typeof(k₂) == Void ? (n₂ = m.nfun(λ)) : (n₂ = m.nfun(λ) + k₂*im)
    phaseangle(one(typeof(n₂)), n₂, θ)
end

function get_complex_ior(m::RefractiveIndex.Material, λ::Real)
    n = m.nfun(λ)
    k = m.kfun(λ)
    n̲ = n + k*im
end

function snellius(n₁::Number, n₂::Number, θᵢ::Real)
    θₜ = asin(n₁ ./ n₂ .* sin(θᵢ))
end


text:
# module Bukdu.Octo

export default, Changeset, change, cast, validate_length

import Base: ==
import ..ApplicationController

default(T::Type, ::Type{String}) = ""
default(T::Type, ::Type{Int}) = 0
default(T::Type, ::Type{Float64}) = 0.0
default(T::Type, ::Type{Float32}) = 0.0

function default(T::Type)::T
    # broadcast #
    # fields = fieldtype.(T,fieldnames(T))
    # T(default.(T, fields)...)
    fields = map(x->fieldtype(T, x), fieldnames(T))
    T(map(x-> default(T, x), fields)...)
end

type Changeset
    model
    changes::Assoc
    function Changeset(model, changes::Assoc)
        T = typeof(model)
        lhs = Assoc(map(x->(x,getfield(model, x)), fieldnames(T)))
        rhs = typed_assoc(T, changes)
        new(model, setdiff(rhs, lhs))
    end
end

function ==(lhs::Changeset, rhs::Changeset)
    T = typeof(lhs.model)
    !isa(rhs.model, T) && return false
    all(x -> ==(getfield(lhs.model, x), getfield(rhs.model, x)), fieldnames(T)) && ==(lhs.changes, rhs.changes)
end

function |>(changeset::Changeset, func::Function)
    func(changeset)
end

function |>(model, func::Function)
    func(change(model))
end

function typed_assoc(T::Type, changes::Assoc)::Assoc
    typ_fieldnames = fieldnames(T)
    Assoc(
        map(filter(kv -> first(kv) in typ_fieldnames, changes)) do kv
            (name,v) = kv
            if name in typ_fieldnames
                fieldT = fieldtype(T, name)
                if isa(v, fieldT)
                    (name, v)
                else
                    (name, parse(fieldT, v))
                end
            end
        end)
end

function cutout_brackets(typ::Type, param::Tuple{Symbol,Any})::Tuple{Symbol,Any}
    t = lowercase(string(typ))
    (key,value) = param
    k = string(key)
    if startswith(k, "$t[") && endswith(k, "]")
        key = Symbol(last(split(k,"["))[1:end-1])
    elseif startswith(k, "$(t)_")
        key = Symbol(k[length("$(t)_")+1:end])
    end
    (key, value)
end

function change{T<:Any,AC<:ApplicationController}(c::AC, model::T)::Changeset
    Changeset(model, Assoc(map(param->cutout_brackets(T,param), c[:query_params])))
end

function change{AC<:ApplicationController}(c::AC, T::Type)::Changeset
    change(c, default(T))
end

function change(model; kw...)::Changeset
    Changeset(model, Assoc(kw))
end

function change(typ::Type; kw...)::Changeset
    change(default(typ); kw...)
end

function cast(changeset::Changeset, params, required_fields)::Changeset
    changeset
end

function validate_length(changeset::Changeset, field::Symbol; kw...)::Changeset
    changeset
end

function cast(params, required_fields)::Function
    (changeset) -> cast(changeset, params, required_fields)
end

function validate_length(field::Symbol; kw...)::Function
    (changeset) -> validate_length(changeset, field; kw...)
end


text:
tree = ParseNewick("(A,B,(C,(D,E)F)G)H;")

@testset "Pruning leaf" begin
    node_list = ["C"]
    pt = prune_tree(tree, node_list)
    @test newick(pt) == newick(ParseNewick("(A,B,((D,E)F)G)H;"))

    node_list = ["C", "C"]
    @test_throws ArgumentError prune_tree(tree, node_list)

    node_list = ["K", "L"]
    @test_throws ArgumentError prune_tree(tree, node_list)

    node_list = ["L", "C"]
    pt = prune_tree(tree, node_list)
    @test newick(pt) == newick(ParseNewick("(A,B,((D,E)F)G)H;"))
end

@testset "Pruning root" begin
    node_list = ["H"]
    @test_throws ArgumentError prune_tree(tree, node_list)
end

@testset "Pruning misc" begin
    node_list1 = ["B", "C"]
    node_list2 = ["F"]
    node_list3 = ["A", "C", "G"]

    prune_tree1 = prune_tree(tree, node_list1)
    prune_tree2 = prune_tree(tree, node_list2)
    prune_tree3 = prune_tree(tree, node_list3)

    @test newick(prune_tree1) == newick(ParseNewick("(A,((D,E)F)G)H;"))
    @test newick(prune_tree2) == newick(ParseNewick("(A,B,(C)G)H;"))
    @test newick(prune_tree3) == newick(ParseNewick("(B)H;"))
end

@testset "Pruning inplace" begin
    node_list1 = ["B", "C"]
    prune_tree!(tree, node_list1)
    @test newick(tree) ==  newick(ParseNewick("(A,((D,E)F)G)H;"))
end


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.
#*
#*--------------------------------------------------------------------
#* Models to simulate a power plant.
#*--------------------------------------------------------------------
#* Author: Argimiro R. Secchi
#* $Id: power_plant.mso 195 2007-03-07 20:30:12Z arge $
#*-------------------------------------------------------------------
# Declaracao de tipos
type SumidouroQ
	SumidouroQ()=begin
		new(
			Potencia(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>0,
				:PosY=>0.5
			)),
			[:Qin,]
		)
	end
	Qin::Potencia
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export SumidouroQ
function atributes(in::SumidouroQ,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/sumidouroQ"
	fields[:Brief]="Sumidouro de calor"
	drive!(fields,_)
	return fields
end
SumidouroQ(_::Dict{Symbol,Any})=begin
	newModel=SumidouroQ()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import Statistics: mean

function getstart(s::TSeries{T}) where T <: Monthly
    if period(s.firstdate) in [1, 4, 7, 10]
        return s.firstdate
    elseif 1 < period(s.firstdate) < 4
        y = year(s.firstdate)
        p = 4
    elseif 4 < period(s.firstdate) < 7
        y = year(s.firstdate)
        p = 7
    elseif 7 < period(s.firstdate) < 10
        y = year(s.firstdate)
        p = 10
    else 7 < period(s.firstdate) < 10
        y = year(s.firstdate) + 1
        p = 1
    end

    return MIT{T}(y*ppy(T) + p - 1)
        # return qq(year(s.firstdate) + 1, 1)
end

function getend(s::TSeries{T}) where T <: Monthly
    if period(lastdate(s)) in [12, 3, 6, 9]
        return lastdate(s)
    elseif 1 <= period(lastdate(s)) < 3
        y = year(lastdate(s)) - 1
        p = 12
    elseif 4 <= period(lastdate(s)) < 6
        y = year(lastdate(s))
        p = 3
    elseif 7 <= period(lastdate(s)) < 9
        y = year(lastdate(s))
        p = 6
    else 10 <= period(lastdate(s)) < 12
        y = year(lastdate(s))
        p = 9
    end

    return MIT{T}(y*ppy(T) + p - 1)
        # return qq(year(s.firstdate) + 1, 1)
end

function getstart(s::TSeries{T}) where T <: Quarterly
    if period(s.firstdate) == 1
        return s.firstdate
    else
        y = year(s.firstdate) + 1
        p = 1

        return MIT{T}(y*ppy(T) + p - 1)
        # return qq(year(s.firstdate) + 1, 1)
    end
end

function getend(s::TSeries{T}) where T <: Quarterly
    if TimeSeriesEcon.period(lastdate(s)) == ppy(s)
        return lastdate(s)
    else
        y = year(lastdate(s)) - 1
        p = ppy(s)
        return MIT{T}(y*ppy(T) + p - 1)
        # return qq(year(lastdate(s)) - 1, ppy(s))
    end
end

function crop(s::TSeries{T}) where T <: Union{Quarterly, Monthly}
    getstart(s) <= getend(s) || error("TSeries can't be cropped. start:$(getstart(s)) > end:$(getend(s))")

    return s[getstart(s):getend(s)]
end

function Base.convert(::Type{TSeries{Yearly}}, s::TSeries{Quarterly})
    values = Vector{Float64}()

    for i in 1:4:length(crop(s))
        push!(values, mean(crop(s)[i:i+3]))
    end

    firstdate = year(crop(s).firstdate) |> yy

    return TSeries(firstdate, values)
end

function Base.convert(::Type{TSeries{Quarterly}}, s::TSeries{Monthly})
    values = Vector{Float64}()

    for i in 1:3:length(crop(s))
        push!(values, mean(crop(s)[i:i+2]))
    end

    y = year(crop(s).firstdate) |> yy
    p = period(crop(s).firstdate) |>
        x -> div(x - 1, 3) + 1

    return TSeries(qq(y, p), values)
end

function Base.convert(::Type{TSeries{Yearly}}, s::TSeries{Monthly})
    convert(TSeries{Quarterly}, s) |>
        x -> convert(TSeries{Yearly}, x)
end

# Arbitrary function

function Base.convert(::Type{TSeries{Yearly}}, s::TSeries{Quarterly}, f::Symbol)
    values = Vector{Float64}()
    for i in 1:4:length(crop(s))
        push!(values, eval(f)(crop(s)[i:i+3]))
    end
    firstdate = year(crop(s).firstdate) |> yy
    return TSeries(firstdate, values)
end

# convert(TSeries{Yearly}, s)
# convert(TSeries{Yearly}, p)

# function mode(a)
#     isempty(a) && throw(ArgumentError("mode is not defined for empty collections"))
#     cnts = Dict{eltype(a),Int}()
#     # first element
#     mc = 1
#     mv, st = iterate(a)
#     cnts[mv] = 1
#     # find the mode along with table construction
#     y = iterate(a, st)
#     while y !== nothing
#         x, st = y
#         if haskey(cnts, x)
#             c = (cnts[x] += 1)
#             if c > mc
#                 mc = c
#                 mv = x
#             end
#         else
#             cnts[x] = 1
#             # in this case: c = 1, and thus c > mc won't happen
#         end
#         y = iterate(a, st)
#     end
#     return mv
# end

# mode([1, 2, 33, 33, 33])


text:
using TensorNetworkEvolve, Yao, Graphs
using LinearAlgebra, Test

@testset "random gate" begin
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end
    # simple
    peps = rand_simplepeps(ComplexF64, g, 2; Dmax=10)
    reg = Yao.ArrayReg(vec(peps))
    rydberg_tebd!(peps, g; t=0.3, C=10.0, Δ=0.2, Ω=0.3, nstep=10)
    @test !isapprox(statevec(peps), statevec(reg); atol=1e-3)
    rydberg_tebd!(reg, g; t=0.3, C=10.0, Δ=0.2, Ω=0.3, nstep=10)
    @test isapprox(statevec(peps), statevec(reg); atol=1e-3)
    @show norm(statevec(peps) .- statevec(reg))

    # vidal
    peps = rand_vidalpeps(ComplexF64, g, 2; Dmax=10, ϵ=1e-10)
    reg = Yao.ArrayReg(vec(peps))
    rydberg_tebd!(peps; t=0.3, C=10.0, Δ=0.2, Ω=0.3, nstep=10)
    @test !isapprox(statevec(peps), statevec(reg); atol=1e-3)
    rydberg_tebd!(reg, g; t=0.3, C=10.0, Δ=0.2, Ω=0.3, nstep=10)
    @test isapprox(statevec(peps), statevec(reg); atol=1e-3)
    @show norm(statevec(peps) .- statevec(reg))
end


text:
using CommandLineEquationsOfState
using Test

@testset "CommandLineEquationsOfState.jl" begin
    # Write your own tests here.
    include("Settings.jl")
end


text:
using WiltonInts84
using StaticArrays
using Test

#include("num_quad.jl")

if !@isdefined record
    record = false
end

v1p = SVector(4.12, 1.74, 1.8);
v2p = SVector(2.04, 3.47, -3.5);
v3p = SVector(-4.21, 3.2, 0.08);
lvp = v2p-v1p; # an edge vector
observationPointp = v1p + 2.84*lvp;
I, J = wiltonints(v1p,v2p,v3p,observationPointp,Val{0})

using JLD2
fn = joinpath(dirname(@__FILE__),"issue1.jld2")

if record == true
    K, L = dblquadints1(v1p, v2p, v3p, observationPointp, Val{0})
    K = [k for k in K]
    L = [L[d][i] for i in 1:3, d in eachindex(L)]
    jldopen(fn, "w") do file
        write(file, "K", K)
        write(file, "L", L)
    end
else
    K, L = jldopen(fn,"r") do file
        read(file, "K"),
        read(file, "L")
    end
end

for d in eachindex(K)
    @eval @test nearlyequal(I[$d], K[$d], 10e-8)
    for i in 1:3
        @eval @test nearlyequal(J[$d][$i], L[$i,$d], 10e-8)
    end
end