text:
# generate examples
import Literate
using Base.Filesystem: cp, basename

println("Start generating Literate.jl examples")

EXAMPLEDIR = joinpath(@__DIR__, "..", "examples")
GENERATEDDIR = joinpath(@__DIR__, "src", "examples")
mkpath(GENERATEDDIR)
for name in readdir(EXAMPLEDIR)
    path = abspath(joinpath(EXAMPLEDIR, name))
    if endswith(path, ".jl")
        example = path
        script = Literate.script(example, GENERATEDDIR)
        code = strip(read(script, String))
        mdpost(str) = replace(str, "@__CODE__" => code)
        Literate.markdown(example, GENERATEDDIR, postprocess=mdpost)
        Literate.notebook(example, GENERATEDDIR, execute=true)
    elseif any(endswith.(path, [".png", ".jpg", ".gif"]))
        cp(path, joinpath(GENERATEDDIR, name); force=true)
    elseif isdir(path)
        folder = path
        target = joinpath(GENERATEDDIR, basename(folder))
        cp(folder, target, force=true)
    else
        @warn "ignoring $name"
    end
end

# remove any .vtu files in the generated dir (should not be deployed)
cd(GENERATEDDIR) do
    foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
    foreach(file -> endswith(file, ".pvd") && rm(file), readdir())
end

println("Finished generating Literate.jl examples")


text:
module Reexport

macro reexport(ex::Expr)
    esc(reexport(__module__, ex))
end

reexport(m::Module, l::LineNumberNode) = l

function reexport(m::Module, ex::Expr)
    # unpack any macros
    ex = macroexpand(m, ex)
    # recursively unpack any blocks
    if ex.head === :block
        return Expr(:block, map(e -> reexport(m, e), ex.args)...)
    end

    ex.head in (:module, :using, :import) ||
        ex.head === :toplevel && all(e -> isa(e, Expr) && e.head === :using, ex.args) ||
        error("@reexport: syntax error")

    if ex.head === :module
        # @reexport {using, import} module Foo ... end
        modules = Any[ex.args[2]]
        ex = Expr(:toplevel, ex, :(using .$(ex.args[2])))
    elseif ex.head in (:using, :import) && ex.args[1].head == :(:)
        # @reexport {using, import} Foo: bar, baz
        symbols = [e.args[end] for e in ex.args[1].args[2:end]]
        return Expr(:toplevel, ex, :(eval(Expr(:export, $symbols...))))
    elseif ex.head === :import && all(e -> e.head === :., ex.args)
        # @reexport import Foo.bar, Baz.qux
        symbols = Any[e.args[end] for e in ex.args]
        return Expr(:toplevel, ex, :(eval(Expr(:export, $symbols...))))
    else
        # @reexport using Foo, Bar, Baz
        modules = Any[e.args[end] for e in ex.args]
    end

    names = GlobalRef(@__MODULE__, :exported_names)
    out = Expr(:toplevel, ex)
    for mod in modules
        push!(out.args, :(Core.eval($m, Expr(:export, $names($mod)...))))
    end
    return out
end

exported_names(m::Module) = filter!(x -> Base.isexported(m, x), names(m; all=true, imported=true))

export @reexport

end # module

ext: jl


text:
# test basic plotting

using Test
using KernelDensityEstimate
using KernelDensityEstimatePlotting
using Colors

@testset "Test basic plotting functions..." begin

# check for errors on plotting code
p = kde!(rand(100));
q = kde!(rand(100).+1.0);

plotKDE([p],c=["red"]);
plot(p);

plot([p;q], c=["red";"green"])

p2 = kde!(rand(3,100));
q2 = kde!(rand(3,100).+1.0);

plot(p2);
plot([p2;q2]);
plot([p2;q2];c=["red";"green"],levels=3);

end


text:
# Contains the physical specification of the model

# A physical specification of a model.
# Contains the geometry information, but no discretization parameters.
# The exact model used is defined by the list of terms.
struct Model{T <: Real}
    # Lattice and reciprocal lattice vectors in columns
    lattice::Mat3{T}
    recip_lattice::Mat3{T}
    unit_cell_volume::T
    recip_cell_volume::T
    n_dim::Int  # Dimension of the system; 3 unless `lattice` has zero columns

    # Electrons, occupation and smearing function
    n_electrons::Int  # usually consistent with `atoms` field, but doesn't have to

    # spin_polarization values:
    #     :none       No spin polarization, αα and ββ density identical,
    #                 αβ and βα blocks zero, 1 spin component treated explicitly
    #     :collinear  Spin is polarized, but everywhere in the same direction.
    #                 αα ̸= ββ, αβ = βα = 0, 2 spin components treated
    #     :full       Generic magnetization, non-uniform direction.
    #                 αβ, βα, αα, ββ all nonzero, different (not supported)
    #     :spinless   No spin at all ("spinless fermions", "mathematicians' electrons").
    #                 The difference with :none is that the occupations are 1 instead of 2
    spin_polarization::Symbol
    n_spin_components::Int  # 2 if :collinear, 1 otherwise

    # If temperature==0, no fractional occupations are used.
    # If temperature is nonzero, the occupations are
    # `fn = max_occ*smearing((εn-εF) / temperature)`
    temperature::T
    smearing::Smearing.SmearingFunction # see Smearing.jl for choices

    # Vector of pairs Element => vector of vec3 (positions, fractional coordinates)
    # Possibly empty. `atoms` just contain the information on the atoms, it's up to
    # the `terms` to make use of it (or not)
    atoms::Vector{Pair{Any, Vector{Vec3{T}}}}

    # each element t must implement t(basis), which instantiates a
    # term in a given basis and gives back a term (<: Term)
    # see terms.jl for some default terms
    term_types::Vector

    # list of symmetries of the model
    symmetries::Vector{SymOp}
end

"""
    Model(lattice; n_electrons, atoms, magnetic_moments, terms, temperature,
                   smearing, spin_polarization, symmetry)

Creates the physical specification of a model (without any
discretization information).

`n_electrons` is taken from `atoms` if not specified

`spin_polarization` is :none by default (paired electrons)
unless any of the elements has a non-zero initial magnetic moment.
In this case the spin_polarization will be :collinear.

`magnetic_moments` is only used to determine the symmetry and the
`spin_polarization`; it is not stored inside the datastructure.

`smearing` is Fermi-Dirac if `temperature` is non-zero, none otherwise

The `symmetries` kwarg allows (a) to pass `true` / `false` to enable / disable
the automatic determination of lattice symmetries or (b) to pass an explicit list
of symmetry operations to use for lowering the computational effort.
The default behaviour is equal to `true`, namely that the code checks the
specified model in form of the Hamiltonian `terms`, `lattice`, `atoms` and `magnetic_moments`
parameters and from these automatically determines a set of symmetries it can safely use.
If you want to pass custom symmetry operations (e.g. a reduced or extended set) use the
`symmetry_operations` function. Notice that this may lead to wrong results if e.g. the
external potential breaks some of the passed symmetries. Use `false` to turn off
symmetries completely.
"""
function Model(lattice::AbstractMatrix{T};
               n_electrons=nothing,
               atoms=[],
               magnetic_moments=[],
               terms=[],
               temperature=T(0.0),
               smearing=nothing,
               spin_polarization=default_spin_polarization(magnetic_moments),
               symmetries=default_symmetries(lattice, atoms, magnetic_moments, terms, spin_polarization),
               ) where {T <: Real}
    lattice = Mat3{T}(lattice)
    temperature = T(austrip(temperature))

    if n_electrons === nothing
        # get it from the atom list
        isempty(atoms) && error("Either n_electrons or a non-empty atoms list should be provided")
        n_electrons = sum(length(pos) * n_elec_valence(el) for (el, pos) in atoms)
    else
        @assert n_electrons isa Int
    end

    # Special handling of 1D and 2D systems, and sanity checks
    n_dim = count(!iszero, eachcol(lattice))
    n_dim > 0 || error("Check your lattice; we do not do 0D systems")
    for i = n_dim+1:3
        norm(lattice[:, i]) == norm(lattice[i, :]) == 0 || error(
            "For 1D and 2D systems, the non-empty dimensions must come first")
    end
    _is_well_conditioned(lattice[1:n_dim, 1:n_dim]) || @warn (
        "Your lattice is badly conditioned, the computation is likely to fail.")

    # Note: In the 1D or 2D case, the volume is the length/surface
    recip_lattice = compute_recip_lattice(lattice)
    unit_cell_volume  = compute_unit_cell_volume(lattice)
    recip_cell_volume = compute_unit_cell_volume(recip_lattice)

    spin_polarization in (:none, :collinear, :full, :spinless) ||
        error("Only :none, :collinear, :full and :spinless allowed for spin_polarization")
    spin_polarization == :full && error("Full spin polarization not yet supported")
    !isempty(magnetic_moments) && !(spin_polarization in (:collinear, :full)) && @warn(
        "Non-empty magnetic_moments on a Model without spin polarization detected."
    )
    n_spin = length(spin_components(spin_polarization))

    if smearing === nothing
        @assert temperature >= 0
        # Default to Fermi-Dirac smearing when finite temperature
        smearing = temperature > 0.0 ? Smearing.FermiDirac() : Smearing.None()
    end

    if !allunique(string.(nameof.(typeof.(terms))))
        error("Having several terms of the same name is not supported.")
    end

    # Determine symmetry operations to use
    symmetries == true  && (symmetries = default_symmetries(lattice, atoms, magnetic_moments,
                                                            terms, spin_polarization))
    symmetries == false && (symmetries = [identity_symop()])
    @assert !isempty(symmetries)  # Identity has to be always present.

    Model{T}(lattice, recip_lattice, unit_cell_volume, recip_cell_volume, n_dim, n_electrons,
             spin_polarization, n_spin, T(temperature), smearing, atoms, terms, symmetries)
end
Model(lattice::AbstractMatrix{T}; kwargs...) where {T <: Integer} = Model(Float64.(lattice); kwargs...)
Model(lattice::AbstractMatrix{Q}; kwargs...) where {Q <: Quantity} = Model(austrip.(lattice); kwargs...)

normalize_magnetic_moment(::Nothing)  = Vec3{Float64}(zeros(3))
normalize_magnetic_moment(mm::Number) = Vec3{Float64}(0, 0, mm)
normalize_magnetic_moment(mm::AbstractVector) = Vec3{Float64}(mm)

"""
:none if no element has a magnetic moment, else :collinear or :full
"""
function default_spin_polarization(magnetic_moments)
    isempty(magnetic_moments) && return :none
    all_magmoms = (normalize_magnetic_moment(magmom)
                   for (_, magmoms) in magnetic_moments
                   for magmom in magmoms)

    all(iszero, all_magmoms) && return :none
    all(iszero(magmom[1:2]) for magmom in all_magmoms) && return :collinear

    :full
end

"""
Default logic to determine the symmetry operations to be used in the model.
"""
function default_symmetries(lattice, atoms, magnetic_moments, terms, spin_polarization;
                        tol_symmetry=1e-5)
    dimension = count(!iszero, eachcol(lattice))
    if spin_polarization == :full || dimension != 3
        return [identity_symop()]  # Symmetry not supported in spglib
    elseif spin_polarization == :collinear && isempty(magnetic_moments)
        # Spin-breaking due to initial magnetic moments cannot be determined
        return [identity_symop()]
    elseif any(breaks_symmetries, terms)
        return [identity_symop()]  # Terms break symmetry
    else
        magnetic_moments = [el => normalize_magnetic_moment.(magmoms)
                            for (el, magmoms) in magnetic_moments]
        return symmetry_operations(lattice, atoms, magnetic_moments,
                                   tol_symmetry=tol_symmetry)
    end
end

"""
Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise).
"""
function filled_occupation(model)
    if model.spin_polarization in (:spinless, :collinear)
        return 1
    elseif model.spin_polarization == :none
        return 2
    else
        error("Not implemented $(model.spin_polarization)")
    end
end

"""
Explicit spin components of the KS orbitals and the density
"""
function spin_components(spin_polarization::Symbol)
    spin_polarization == :collinear && return (:up, :down  )
    spin_polarization == :none      && return (:both,      )
    spin_polarization == :spinless  && return (:spinless,  )
    spin_polarization == :full      && return (:undefined, )
end
spin_components(model::Model) = spin_components(model.spin_polarization)

_is_well_conditioned(A; tol=1e5) = (cond(A) <= tol)

"""
Compute the reciprocal lattice. Takes special care of 1D or 2D cases.
We use the convention that the reciprocal lattice is the set of G vectors such
that G ⋅ R ∈ 2π ℤ for all R in the lattice.
"""
function compute_recip_lattice(lattice::AbstractMatrix{T}) where {T}
    # Note: pinv pretty much does the same, but the implied SVD causes trouble
    #       with interval arithmetic and dual numbers, so we go for this version.
    n_dim = count(!iszero, eachcol(lattice))
    @assert 1 ≤ n_dim ≤ 3
    if n_dim == 3
        2T(π) * inv(lattice')
    else
        2T(π) * Mat3{T}([
            inv(lattice[1:n_dim, 1:n_dim]')   zeros(T, n_dim, 3 - n_dim);
            zeros(T, 3 - n_dim, 3)
        ])
    end
end

"""
Compute unit cell volume volume. In case of 1D or 2D case, the volume is the length/surface.
"""
function compute_unit_cell_volume(lattice)
    n_dim = count(!iszero, eachcol(lattice))
    abs(det(lattice[1:n_dim, 1:n_dim]))
end


text:
using TensorToolbox
using Base.Test


text:
# implements a rollout policy
module RolloutPolicies

using DiscreteMDPs

export RolloutPolicy
export action

import DiscreteMDPs.Policy
import DiscreteMDPs.action

type RolloutPolicy <: Policy
    d::Int64               # search depth
    n::Int64               # number of iterations
    c::Float64             # explortaion constant
    Q::Vector{Float64}     # Q values at the root state
    N::Vector{Int64}       # Number of visits to each state
    actions::Vector{Int64} # Vector of allowable actions
    mdp::DiscreteMDP
end

function RolloutPolicy(mdp::DiscreteMDP; d::Int64=20, n::Int64=100, c::Float64=0.1)
    nActions = numActions(mdp)
    Q = zeros(nActions)
    N = zeros(Int64, nActions)
    actions = [1:nActions]
    return RolloutPolicy(d, n, c, Q, N, actions, mdp)
end

function action(policy::RolloutPolicy, sIdx::Int64)
    n = policy.n
    fill!(policy.Q, 0.0)
    fill!(policy.N, 0)
    for i = 1:n
        simulate(policy, sIdx)
    end
    a = policy.actions
    return a[indmax(policy.Q)]
end

function simulate(policy::RolloutPolicy, s::Int64)
    Q = policy.Q
    N = policy.N
    c = policy.c
    d = policy.d
    actions = policy.actions
    mdp = policy.mdp
    discount = mdp.discount
    i = indmax(Q + c .* real(sqrt(complex(log(sum(N)) ./ N))))
    a = actions[i]
    sp = nextState(mdp, s, a)
    q = reward(mdp, s, a) + discount*rollout(policy, sp, d)
    N[i] += 1
    Q[i] += (q - Q[i]) / N[i]
    return
end

function rollout(policy::RolloutPolicy, s::Int64, depth::Int64)
    mdp = policy.mdp
    discount = mdp.discount
    if depth == 0
        return 0.0
    else
        a  = getAction(mdp, s)
        sp = nextState(mdp, s, a)
        return reward(mdp, s, a) + discount*rollout(policy, sp, depth-1)
    end
end

end # module


text:
using Pythia
pythia = newPythia()

icxx"""$pythia->readString("Beams:idA = 11");"""
icxx"""$pythia->readString("Beams:idB = -11");"""
icxx"""$pythia->readString("Beams:eCM = 10.");"""
icxx"""$pythia->readString("WeakSingleBoson:ffbar2ffbar(s:gm) = on");"""
icxx"""$pythia->init();"""

# Check if this is a e⁺e⁻ -> qqbar process
function isqqbar(e)
    # Look through the particles to find those in the hard process
    for p in particles(e)
        if abs(icxx"$p.status();") == 23
            id = icxx"$p.id();"
            # quarks
            return 1 <= abs(id) <= 5
        end
    end
end

count_hadrons(e) = sum([icxx"$p.isHadron();" for p in particles(e)])

energies = linspace(10,120,6)
nhadrons = map(energies) do eCM
   pythia = newPythia()

   icxx"""$pythia->readString("Beams:idA = 11");"""
   icxx"""$pythia->readString("Beams:idB = -11");"""
   readString(pythia,"Beams:eCM = $eCM")
   icxx"""$pythia->readString("WeakSingleBoson:ffbar2ffbar(s:gm) = on");"""
   icxx"""$pythia->init();"""

   mean(filter(x->x != -1, map(e->isqqbar(e) ? count_hadrons(e) : -1, take(events(pythia),1000))))
end
plot(x = energies, y = nhadron, Geom.point, Geom.smooth(method=:loess,smoothing=0.9), Guide.xlabel("√s [GeV]"), Guide.ylabel("⟨Nhadron⟩"))


text:
#=
JoiceTypeConstructorFunctions:
- Julia version: 1.6
- Author: connorforsythe
- Date: 2021-11-05
=#

export ConstructModelInfo

function checkModelInfo(Data::DataFrame, ChoiceColumn::Symbol, Parameters::Vector{Symbol}, ChoiceSetIDColumn::Symbol, RandomParameters::RandColsType=nothing, PanelIDColumn::ColsType=nothing, WeightColumn::ColsType=nothing)

    dataCols = Symbol.(names(Data))
    missingCols = []
    if !(ChoiceColumn in dataCols)
        push!(missingCols, ChoiceColumn)
    end
    if !(ChoiceSetIDColumn in dataCols)
        push!(missingCols, ChoiceSetIDColumn)
    end
    if !(PanelIDColumn in dataCols) && (PanelIDColumn!=nothing)
        push!(missingCols, PanelIDColumn)
    end

    for tempSymbol in Parameters
        if !(tempSymbol in dataCols)
            push!(missingCols, tempSymbol)
        end
    end

    if RandomParameters!=nothing
        for tempSymbol in keys(RandomParameters)
            if !(tempSymbol in dataCols)
                push!(missingCols, tempSymbol)
            end
        end
    end

    if length(missingCols)>0
        error("You have specified the following cols that aren't in the provided dataset: $missingCols")
    end
end

function convertParameterColumns(x::Vector{String}) return Symbol.(x) end
function convertParameterColumns(x::Vector{Symbol}) return x end
function convertParameterColumns(x::Symbol) return [x] end
function convertParameterColumns(x::String) return convertParameterColumns(Symbol(x)) end
function convertParameterColumns(x::Nothing) return x end
function convertParameterColumn(x::String) return Symbol(x) end
function convertParameterColumn(x::Symbol) return x end
function convertParameterColumn(x::Nothing) return x end
function convertRandomParameterColumns(x::Dict{Symbol, String}) return x end
function convertRandomParameterColumns(x::Dict{String, String}) return Dict(zip(convertParameterColumns([keys(x)...]), values(x))) end
function convertRandomParameterColumns(x::Nothing) return x end

function ConstructModelInfo(Data::DataFrame, ChoiceColumn::ColType, Parameters::ColsType, ChoiceSetIDColumn::ColType; RandomParameters::RandColsType=nothing, PanelIDColumn::ColType=nothing,
    Space::String="preference", Name::NameType=nothing, InterDrawCount::Int=2048, WeightColumn::ColType=nothing)

    TimeCreated = Dates.now()

    ChoiceColumn = convertParameterColumn(ChoiceColumn)
    Parameters = convertParameterColumns(Parameters)
    ChoiceSetIDColumn = convertParameterColumn(ChoiceSetIDColumn)
    WeightColumn = convertParameterColumn(WeightColumn)
    PanelIDColumn = convertParameterColumn(PanelIDColumn)
    RandomParameters = convertRandomParameterColumns(RandomParameters)
    RandomDraws = getDraws(RandomParameters, InterDrawCount)
    typedSpace::JoiceModelSpace = Preference()

    colsToKeep = copy(Parameters)

    appendParamNames!(colsToKeep, RandomParameters)

    push!(colsToKeep, ChoiceColumn)
    push!(colsToKeep, ChoiceSetIDColumn)

    if(lowercase(Space)=="wtp")
        typedSpace = WTP()
    end

    typedDataStructure::JoiceDataStructure = Panel()
    ObservationColumn::Symbol = ChoiceSetIDColumn
    if(PanelIDColumn==nothing)
        typedDataStructure = Cross()
    else
        push!(colsToKeep, PanelIDColumn)
        ObservationColumn = PanelIDColumn
    end

    if(WeightColumn!=nothing)
        push!(colsToKeep, WeightColumn)
    end

    checkModelInfo(Data, ChoiceColumn, Parameters, ChoiceSetIDColumn, RandomParameters, PanelIDColumn, WeightColumn)
    PreferenceSearchBounds = [-1, 1]

    println("Model being constructed with the following properties:")
    if(PanelIDColumn!=nothing)
        println("Panel Variable: $PanelIDColumn")
    end
    println("Choice Set Variable: $ChoiceSetIDColumn")
    println("Choice Indicator Variable: $ChoiceColumn")
    println("Fixed Parameters: $Parameters")
    if(RandomParameters!=nothing)
        println("Random Parameters: $RandomParameters")
    end
    print("Data has ")
    print(size(Data)[1])
    println(" rows")

    cleanDataFrame = cleanData(Data, colsToKeep, ChoiceSetIDColumn, PanelIDColumn)

    r = ModelInfo(cleanDataFrame, ChoiceColumn, Parameters, :choiceSetID, RandomParameters, RandomDraws,
                    PanelIDColumn, ObservationColumn, typedDataStructure, typedSpace, Name, InterDrawCount,
                    WeightColumn, TimeCreated, PreferenceSearchBounds)

    return r
end

function cleanData(data::DataFrame, colsToKeep::Vector{Symbol}, choiceSetIDColumn::Symbol, panelIDColumn::Symbol)

    data = data[:, colsToKeep]

    #Construct new id to allow for repeated choice set ids given
    maxChoiceSetID = max(data[:, choiceSetIDColumn]...)
    panelIDMultiplier = 10^(ceil(log10(maxChoiceSetID)))

    data[:, :choiceSetID] .= (data[:, panelIDColumn].*panelIDMultiplier).+data[:, choiceSetIDColumn]

    return data

end

function cleanData(data::DataFrame, colsToKeep::Vector{Symbol}, choiceSetIDColumn::Symbol, panelIDColumn::Nothing)

    data = data[:, colsToKeep]

    data[:, :choiceSetID] .= data[:, choiceSetIDColumn]

    return data

end

function appendParamNames!(currentVect::Vector{Symbol}, randomParameters::Nothing) end

function appendParamNames!(currentVect::Vector{Symbol}, randomParameters::Dict{Symbol, String})
    for tempCol in keys(randomParameters)
        push!(currentVect, tempCol)
    end
end

max_stars_repo_name: crforsythe/Joice-Releases

max_issues_repo_name: crforsythe/Joice-Releases

max_forks_repo_name: crforsythe/Joice-Releases


text: using MartaCT, Test


text:
module CDF

# loads a graph from a IEEE CDF file.
# http://www2.ee.washington.edu/research/pstca/formats/cdf.txt
# http://www2.ee.washington.edu/research/pstca/pf30/ieee30cdf.txt

using LightGraphs
using LightGraphs: AbstractGraphFormat

import LightGraphs: loadgraph, loadgraphs, savegraph

export CDFFormat

struct CDFFormat <: AbstractGraphFormat end

function _loadcdf(io::IO)
    srcs = Vector{Int}()
    dsts = Vector{Int}()
    vertices = Vector{Int}()
    inbusdata = false
    inbranchdata = false
    while !eof(io)
        line = strip(chomp(readline(io)))
        if inbusdata
            if occursin("-999", line)
                inbusdata = false
            else
                v = parse(Int, split(line)[1])
                push!(vertices, v)
            end
        elseif inbranchdata
            if occursin("-999", line)
                inbranchdata = false
            else
                (src_s, dst_s) = split(line)[1:2]
                src = something(findfirst(isequal(parse(Int, src_s)), vertices), 0)
                dst = something(findfirst(isequal(parse(Int, dst_s)), vertices), 0)
                push!(srcs, src)
                push!(dsts, dst)
            end
        else
            inbusdata = startswith(line, "BUS DATA FOLLOWS")
            inbranchdata = startswith(line, "BRANCH DATA FOLLOWS")
        end
    end
    n_v = length(vertices)
    g = LightGraphs.Graph(n_v)
    for p in zip(srcs, dsts)
        add_edge!(g, p)
    end
    return g
end

loadcdf(io::IO, gname::String) = _loadcdf(io)
loadgraph(io::IO, gname::String, ::CDFFormat) = loadcdf(io, gname)
loadgraphs(io::IO, ::CDFFormat) = Dict("graph" => loadcdf(io, "graph"))

end # module


text:
function identical_matrix(n::Integer)
    ret = zeros(n, n)
    for itr1 in 1 : n
        ret[itr1, itr1] = 1
    end
    return ret
end

function inner_product(mat::Matrix, a::Integer, b::Integer)
    # if $a$ or $b$ are out of bound, return $0$
    if a > size(mat, 2) || b > size(mat, 2)
        return 0
    end
    # initialize the result
    ret = 0
    # calculate the inner (dot) product of column $a$ and $b$
    for itr1 in 1 : size(mat, 1)
        ret += mat[itr1, a] * mat[itr1, b]
    end
    return ret
end

function column_elimination!(mat::Matrix, a::Integer, b::Integer, k::Number, row::Integer = 1)
    for itr1 in row : size(mat, 1)
        mat[itr1, b] += mat[itr1, a] * k
    end
    return mat
end

function mgs(mat::Matrix)
    # make a copy of the matrix as $Q$
    Q = float(copy(mat))
    # it's assumed that the matrix has $m$ rows and $n$ columns
    m, n = size(mat)
    # make a identical matrix as $R$
    R = identical_matrix(n)
    for itr1 in 1 : n
        # calculate the norm of the column vector
        norm = sqrt(inner_product(Q, itr1, itr1))
        # if the column is $0$, just skip it
        norm == 0 && continue
        # normalize the column vector
        for itr2 in 1 : m
            Q[itr2, itr1] /= norm
        end
        # multiply the $itr1$th row with the norm
        for itr2 in itr1 : n
            R[itr1, itr2] *= norm
        end
        # substract the projections of the current vector
        # from the remaining vectors
        for itr2 in (itr1 + 1) : n
            scaler = inner_product(Q, itr1, itr2)
            column_elimination!(Q, itr1, itr2, -scaler)
            # add the corresponding vectors to some rows
            for itr3 in itr2 : n
                R[itr1, itr3] += R[itr2, itr3] * scaler
            end
        end
    end
    return Q, R
end

max_stars_repo_name: Shimushushushu/Etude-Rhapsody

max_issues_repo_name: Shimushushushu/Etude-Rhapsody

max_forks_repo_name: Shimushushushu/Etude-Rhapsody


text:
using RandomStrainDistributions
using Documenter

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

makedocs(;
    modules=[RandomStrainDistributions],
    authors="W. Joe Meese <meese022@umn.edu> and contributors",
    repo="https://github.com/meese-wj/RandomStrainDistributions.jl/blob/{commit}{path}#{line}",
    sitename="RandomStrainDistributions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
## Quantile Plot

"""
    qqbuild()
"""
function qqbuild(d::UnivariateDistribution, x::AbstractVector)
    n = length(x)
    ind = collect(1:n)
    grid = ind ./ (n+1)

    qd = quantile.(Ref(d), grid)
    qx = sort(x, rev=false)

    return Distributions.QQPair(qd, qx)   # QQpair(x=Model, y=Empirical)
end

"""
    qqplot()

Quantile-Quantile Plot
"""
function qqplot(qq::QQPair, elements::Gadfly.ElementOrFunction...)
    return plot(x=qq.qx,  # Model
                y=qq.qy,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                #Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                #Guide.title("Quantile Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

function qqplot(fm::UnivariateDistribution, data::AbstractVector, elements::Gadfly.ElementOrFunction...)
    qq = qqbuild(fm, data)
    return plot(x=qq.qx,  # Model
                y=qq.qy,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                Guide.title("Quantile Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

## Non-stationnary Quantile Plot (À RETRAVAILLER)

struct QQPairNS{U<:AbstractVector, V<:AbstractVector}
    qx::U
    qy::V
end

"""
    qqbuild() (standardized)
"""
function qqbuild(fm::Array{<:UnivariateDistribution,1}, data::Array{<:Real,1})
    ztilde = -log.(1 .- cdf.(fm, data));  # Transformation to a standard exponential distribution
    ztildei = sort(ztilde)

    m = length(ztilde)
    x = -log.(1 .- ((1:m) / (m+1)));

    return QQPairNS(x, ztildei)
end

function qqplot(qq::QQPairNS, elements::Gadfly.ElementOrFunction...)
    return plot(x=qq.qx,  # Model
                y=qq.qy,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                #Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                #Guide.title("Residual Quantile Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

function qqplot(fm::Array{<:UnivariateDistribution,1}, data::Array{<:Real,1}, elements::Gadfly.ElementOrFunction...)
    qq = qqbuild(fm, data)
    return plot(x=qq.qx,  # Model
                y=qq.qy,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                Guide.title("Residual Quantile Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

## Probability Plot

struct PPPair{U<:AbstractVector, V<:AbstractVector}
    px::U
    py::V
end

"""
    ppbuild()
"""
function ppbuild(d::UnivariateDistribution, x::AbstractVector)

    n = length(x)
    ind = collect(1:n)
    px = ind ./ (n+1)
    pd = cdf.(Ref(d), sort(x, rev=false))

    return PPPair(pd, px)
end

"""
    probplot()

Probability plot
"""
function probplot(pp::PPPair, elements::Gadfly.ElementOrFunction...)
    return plot(x=pp.px,  # Model
                y=pp.py,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                #Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                #Guide.title("Probability Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

function probplot(fm::UnivariateDistribution, data::AbstractVector, elements::Gadfly.ElementOrFunction...)
    pp = ppbuild(fm, data)
    return plot(x=pp.px,  # Model
                y=pp.py,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                Guide.title("Probability Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

## Non-stationnary Probability Plot (À RETRAVAILLER)

struct PPPairNS{U<:AbstractVector, V<:AbstractVector}
    px::U
    py::V
end

"""
    ppbuild() (standardized)
"""
function ppbuild(fm::Array{<:UnivariateDistribution,1}, data::Array{<:Real,1})
    # TO-DO : - vérifier l'ordre des variables (x vs y)

    ztilde = -log.(1 .- cdf.(fm, data));  # Transformation to a standard exponential distribution
    ztildei = sort(ztilde)

    m = length(ztilde)
    x = (1:m) ./ (m+1);
    y = 1 .- exp.(-ztildei)

    return PPPairNS(y, x)
end

function probplot(pp::PPPairNS, elements::Gadfly.ElementOrFunction...)
    return plot(x=pp.px,  # Model
                y=pp.py,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                #Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                #Guide.title("Residual Probability Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

function probplot(fm::Array{<:UnivariateDistribution,1}, data::Array{<:Real,1}, elements::Gadfly.ElementOrFunction...)
    pp = ppbuild(fm, data)
    return plot(x=pp.px,  # Model
                y=pp.py,  # Empirical
                Geom.point,
                Geom.abline(color="red", style=:dash),
                Guide.xlabel("Model"), Guide.ylabel("Empirical"),
                Guide.title("Residual Probability Plot"),
                Theme(discrete_highlight_color=c->nothing),
                elements...)
end

## Histogram Plot

"""
    histplot(fm::UnivariateDistribution, data::AbstractVector)

Histogram plot
"""
function histplot(fm::UnivariateDistribution, data::AbstractVector)

    #TODO: arranger max/min pour la grille

    nbin = floor(sqrt(length(data)))
    d = layer(x -> pdf(fm, x), 0, maximum(data), Geom.line, Theme(default_color="red"))
    h = layer(x = data, Geom.histogram(density=true, bincount=nbin))
    return plot(d,
                h,
                Coord.cartesian(xmin = minimum(data), xmax = maximum(data)),
                Guide.xlabel("Data"),
                Guide.ylabel("Density"),
                Guide.title("Density Plot")
                )
end

## Return Level Plot

function ecdf(y::Vector{<:Real})::Tuple{Vector{<:Real}, Vector{<:Real}}
    ys = sort(y)
    n = length(ys)
    p = collect(1:n)/(n+1)

    return ys, p
end

function returnlevelplot_data(fm::UnivariateDistribution, data::AbstractVector)
    y, p = ecdf(data)

    T = 1 ./ (1 .- p)

    n = length(y)

    q = Vector{Float64}(undef, n)

    for i in eachindex(p)
       q[i] = quantile(fm, p[i])[1]
    end

    return y, T, q

end

"""
    returnlevelplot(fm::UnivariateDistribution, data::AbstractVector)

Return level plot
"""
function returnlevelplot(fm::UnivariateDistribution, data::AbstractVector)
    y, T, q = returnlevelplot_data(fm, data)
    l1 = layer(x=T, y=q,Geom.line, Theme(default_color="red", line_style=[:dash]))
    l2 = layer(x=T, y=y, Geom.point)
    return plot(l1,
                l2,
                Scale.x_log10,
                Guide.xlabel("Return Period"),
                Guide.ylabel("Return Level"),
                Guide.title("Return Level Plot"),
                Theme(discrete_highlight_color=c->nothing))
end

## Diagnostic plots

"""
    diagnosticplots(fm::UnivariateDistribution, data::AbstractVector)

"""
function diagnosticplots(fm::UnivariateDistribution, data::AbstractVector)::Gadfly.Compose.Context
    f1 = probplot(fm, data)
    f2 = qqplot(fm, data)
    f3 = histplot(fm, data)
    f4 = returnlevelplot(fm, data)

    return gridstack([f1 f2; f3 f4])
end


text:
# Load Julia packages (libraries) needed

using StatisticalRethinking, CmdStan
#gr(size=(600,600));

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path("..", "scripts", "02")
cd(ProjDir)

# Define the Stan language model

binomialstanmodel = "
// Inferring a Rate
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

stanmodel = Stanmodel(name="binomial", monitors = ["theta"], model=binomialstanmodel,
  output_format=:mcmcchains);

# Use 16 observations

N2 = 15
d = Binomial(9, 0.66)
n2 = Int.(9 * ones(Int, N2));

# Show first 5 (generated) observations

k2 = rand(d, N2);
k2[1:min(5, N2)]

# Input data for cmdstan

binomialdata = Dict("N" => length(n2), "n" => n2, "k" => k2);

# Sample using cmdstan

rc, chn, cnames = stan(stanmodel, binomialdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME);

# Describe the draws

MCMCChains.describe(chn)

# Allocate array of Normal fits

fits = Vector{Normal{Float64}}(undef, 4)
for i in 1:4
  fits[i] = fit_mle(Normal, convert.(Float64, chn.value[:, 1, i]))
  println(fits[i])
end

# Plot the 4 chains

mu_avg = sum([fits[i].μ for i in 1:4]) / 4.0;
sigma_avg = sum([fits[i].σ for i in 1:4]) / 4.0;

if rc == 0
  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
  x = 0:0.001:1
  for i in 1:4
    vals = convert.(Float64, chn.value[:, 1, i])
    μ = round(fits[i].μ, digits=2)
    σ = round(fits[i].σ, digits=2)
    p[i] = density(vals, lab="Chain $i density",
       xlim=(0.45, 1.0), title="$(N2) data points")
    plot!(p[i], x, pdf.(Normal(fits[i].μ, fits[i].σ), x), lab="Fitted Normal($μ, $σ)")
  end
  plot(p..., layout=(4, 1))
  #savefig("m2.1s.pdf")
end

# Show the hpd region

hpd(chn, alpha=0.055)

# Compute the hpd bounds for plotting

d, p, c = size(chn);
theta = convert(Vector{Float64}, reshape(chn.value, (d*p*c)));
bnds = quantile(theta, [0.045, 0.945])

# Show hpd region

println("hpd bounds = $bnds\n")

# quadratic approximation

# Compute MAP, compare with CmndStan & MLE

tmp = convert(Array{Float64,3}, chn.value)
draws = reshape(tmp, (size(tmp, 1)*size(tmp, 3)),)

# Compute MAP

using Optim

x0 = [0.5]
lower = [0.2]
upper = [1.0]

inner_optimizer = GradientDescent()

function loglik(x)
  ll = 0.0
  ll += log.(pdf.(Beta(1, 1), x[1]))
  ll += sum(log.(pdf.(Binomial(9, x[1]), k2)))
  -ll
end

res = optimize(loglik, lower, upper, x0, Fminbox(inner_optimizer))

# Summarize mean and sd estimates

# CmdStan mean and sd:

[mean(chn.value), std(chn.value)]

# MAP estimate and associated sd:

[Optim.minimizer(res)[1], std(draws, mean=mean(chn.value))]

# MLE of mean and sd:

[mu_avg, sigma_avg]

# Turing Chain &  89% hpd region boundaries

plot( x, pdf.(Normal( mu_avg , sigma_avg  ) , x ),
xlim=(0.0, 1.2), lab="Normal approximation using MLE")
plot!( x, pdf.(Normal( Optim.minimizer(res)[1] , std(draws, mean=mean(chn.value))) , x),
lab="Normal approximation using MAP")
density!(draws, lab="CmdStan chain")
vline!([bnds[1]], line=:dash, lab="hpd lower bound")
vline!([bnds[2]], line=:dash, lab="hpd upper bound")

# End of `02/clip_08s.jl`


text:
using Test
using CliOptions

testfiles = [ent for ent in readdir(dirname(@__FILE__))
             if endswith(ent, ".jl") && ent != basename(@__FILE__)]

targets = String[]
let i = 1
    while i ≤ length(ARGS)
        if ARGS[i] == "-h" || ARGS[i] == "--help"
            println("Usage: julia \"$PROGRAM_FILE\" [TEST_NAME...]")
            println()
            println("Available Tests:")
            for testfile in testfiles
                println("    " * testfile[1:end-3])
            end
            exit()
        else
            for t in testfiles
                if startswith(t, ARGS[i])
                    push!(targets, t)
                end
            end
        end

        i += 1
    end
end
if length(targets) == 0
    append!(targets, testfiles)
end

@testset "CliOptions" begin
    foreach(t -> include(t), targets)
end


text:
using LinearAlgebra#dot to get the answer
#Julia core using BigFloat
# begin at 1

fname = "results.txt"
function printtofile(fout,result)#file output
    println(fout,result)
end

function commandMode()
    fout=open(fname,"w")
        println("您正在使用的是 命令行 模式，请输入您的向量的维度")
        dim = parse(Int,readline())
        vector = zeros(BigFloat,dim)
        println("请输入您的向量个数")
        n = parse(Int,readline())
        vectors = Vector{Any}(undef,n)
        println("请输入的时候以行为单位，每个数字请回车键入，不要用空格分割两个数字")
        for i in 1:n
            for j in 1:dim
                vector[j] = parse(BigFloat,readline())
            end
            vectors[i]=copy(vector)
        end
        k = Int(n/2)
        t1=time_ns()
        for i = 1:k
            result = dot(vectors[2*i],vectors[2*i-1])
            printtofile(fout,result)
        end
        t2=time_ns()
        println(string((t2-t1)/1e9,"seconds to finish this job"))
        close(fout)
end
function fileMode()
    fout=open(fname,"w")
    println("您正在使用的是 文件 模式，请输入您的向量的维度")
    dim = parse(Int,readline())
    vector = zeros(BigFloat,dim)
    println("请输入您的向量个数")
    n = parse(Int,readline())
    vectors = Vector{Any}(undef,n)
    println("请输入您想打开的文件路径，使用绝对路径")
    pathfile = readline()
    fin = open(pathfile,"r")
    println("正在读取文件请稍后")
    for i in 1:n
        for j in 1:dim
            vector[j] = BigFloat(readline(fin))
        end
        vectors[i]=copy(vector)
    end
    println("正在进行计算")
    k = Int(n/2)
    t1=time_ns()
    for i = 1:k
        result = dot(vectors[2*i],vectors[2*i-1])
        printtofile(fout,result)
    end
    t2=time_ns()
    println(string((t2-t1)/1e9-0.03*k,"seconds to finish this job"))
    close(fout)
end
println("Julia核心已经初始化完成，正在进入主程序")
println("欢迎继续使用本程序，输入-q退出程序，输入-f使用文件读取模式，输入-c进入命令行读取模式")
global cmd = readline()
while cmd!="-q"
    if cmd[1]!='-'
        print("wrong command, please try again")
    elseif cmd[2]=='f'
    fileMode()
    elseif cmd[2]=='c'
    commandMode()
    else
        print("wrong command, please try again")
    end
    println("请输入您的下一个指令，-q退出程序，输入-f使用文件读取模式，输入-c进入命令行读取模式")
    global cmd = readline()
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ b6b1b198-7bd8-11eb-10f5-a98dd0a59fcf
begin
	cd("..")
	using Pkg
	Pkg.activate(".")
end

# ╔═╡ 594227e6-7bd9-11eb-3ffb-ab8e6e7b4d95
begin
	using ColabINAOE2021
	using PyPlot
	using SparseArrays
	using LinearAlgebra
end

# ╔═╡ 2141105e-7bda-11eb-3519-1d5cff10d945
ion()

# ╔═╡ 6310aeb4-7bd9-11eb-0bd4-81a074081964
# Global Constants
begin
	k = 1
	N = 32
end;

# ╔═╡ 81e6e088-7bd9-11eb-17a6-03ef56dc026b
# Initial amplitude vector constructors

pint(q) = sparsevec(Dict(q + 1 => 1), N + 1); # Single amplitude point at 0 ≤ q ≤ N

# ╔═╡ e3ceb6d0-7bda-11eb-115f-4584602f29df
# Initial amplitude vector

q0 = normalize(pint(8) + pint(24));

# ╔═╡ 203134d6-7bdb-11eb-2df1-374c190cc45b
t_list = range(0, stop = 61, length = 200);

# ╔═╡ b2181d84-7fcf-11eb-2400-23cdf5b8fff7
md"### 1. Circular array"

# ╔═╡ b49a93f8-806b-11eb-2dc1-3f8e020d0f8c
evolc = ColabINAOE2021.Qc(1, t_list, Array(q0));

# ╔═╡ e30513a6-7bdb-11eb-0e31-cd6f16283698
begin
	fig_1 = figure(figsize = (10,5))
	ax_1 = gca()
	ax_1.imshow(evolc)
	ax_1.set_xlabel("Evolution time: t")
	ax_1.set_ylabel("Element index: n")
	ax_1.set_aspect("2")
	ax_1.set_xticks([0,33,66,99,132,165,200])
	ax_1.set_xticklabels([0,10,20,30,40,50,60])

	pcm1 = ax_1.get_children()[10]
	cb1 = colorbar(pcm1,orientation="vertical",shrink=0.7,aspect=35,fraction=0.015)

	tight_layout()
	#savefig("notebooks/circ.png",dpi=150,transparent=true)
end

# ╔═╡ 07fccfdc-7bdc-11eb-1a00-a345a06f8517
sum(abs2.(evolc))

# ╔═╡ ac233d02-7fcf-11eb-1640-c172d67a0cfc
md"### 2. Lineal array"

# ╔═╡ cbc567c2-8095-11eb-0996-05753de896af
evoll = ColabINAOE2021.Ql(1, t_list, Array(q0));

# ╔═╡ fb886f34-8095-11eb-0b4f-e58f8dcb8964
begin
	fig_2 = figure(figsize = (10,5))
	ax_2 = gca()
	ax_2.imshow(evoll)
	ax_2.set_xlabel("Evolution time: t")
	ax_2.set_ylabel("Element index: n")
	ax_2.set_aspect("2")
	ax_2.set_xticks([0,33,66,99,132,165,200])
	ax_2.set_xticklabels([0,10,20,30,40,50,60])

	pcm2 = ax_2.get_children()[10]
	cb2 = colorbar(pcm2,orientation="vertical",shrink=0.7,aspect=35,fraction=0.015)

	tight_layout()
	#savefig("notebooks/lin.png",dpi=150,transparent=true)
end

# ╔═╡ 64040048-8096-11eb-2c6f-ed400b9dbf42
md"### 3. Kravchuk array "

# ╔═╡ 6ec1d92e-8096-11eb-3980-7107878eecb7
evolk = ColabINAOE2021.Qk(t_list, Array(q0));

# ╔═╡ 7f1a6296-8096-11eb-12b1-a14e3a7eb885
begin
	fig_3 = figure(figsize = (10,5))
	ax_3 = gca()
	ax_3.imshow(evolk)
	ax_3.set_xlabel("Evolution time: t")
	ax_3.set_ylabel("Element index: n")
	ax_3.set_aspect("2")
	ax_3.set_xticks([0,33,66,99,132,165,200])
	ax_3.set_xticklabels([0,10,20,30,40,50,60])

	pcm3 = ax_2.get_children()[10]
	cb3 = colorbar(pcm3,orientation="vertical",shrink=0.7,aspect=35,fraction=0.015)

	tight_layout()
	#savefig("notebooks/krav.png",dpi=150,transparent=true)
end

# ╔═╡ Cell order:
# ╠═b6b1b198-7bd8-11eb-10f5-a98dd0a59fcf
# ╠═594227e6-7bd9-11eb-3ffb-ab8e6e7b4d95
# ╠═2141105e-7bda-11eb-3519-1d5cff10d945
# ╠═6310aeb4-7bd9-11eb-0bd4-81a074081964
# ╠═81e6e088-7bd9-11eb-17a6-03ef56dc026b
# ╠═e3ceb6d0-7bda-11eb-115f-4584602f29df
# ╠═203134d6-7bdb-11eb-2df1-374c190cc45b
# ╟─b2181d84-7fcf-11eb-2400-23cdf5b8fff7
# ╠═b49a93f8-806b-11eb-2dc1-3f8e020d0f8c
# ╠═e30513a6-7bdb-11eb-0e31-cd6f16283698
# ╠═07fccfdc-7bdc-11eb-1a00-a345a06f8517
# ╟─ac233d02-7fcf-11eb-1640-c172d67a0cfc
# ╠═cbc567c2-8095-11eb-0996-05753de896af
# ╠═fb886f34-8095-11eb-0b4f-e58f8dcb8964
# ╟─64040048-8096-11eb-2c6f-ed400b9dbf42
# ╠═6ec1d92e-8096-11eb-3980-7107878eecb7
# ╠═7f1a6296-8096-11eb-12b1-a14e3a7eb885


text:
function MapFeature(X1,X2,d)

    m=size(X1,1)
    n=((d+d^2)/2)+d

    p=zeros(m,round(Int,n))
    count=0
    for i =1:d
        for j=0:i
            count=count+1
            p[:,count]=(X1.^(i-j)).*(X2.^j)
        end
    end

    return p
end

function PolyFeature(X,d)

    m=size(X,1)
    n=size(X,2)
    p=zeros(m,d*n)
c=0
for j=1:n
    for i=1:d
        c=c+1
        p[:,c]=X[:,j].^i
    end
end
    return p

end


text:
#! format: off

using SIIPExamples;
using PowerSystems;

print_struct(Bus)

print_tree(PowerSystems.IS.InfrastructureSystemsType)

print_tree(TimeSeriesData)

print_struct(Deterministic)

print_struct(System)

BASE_DIR = abspath(joinpath(dirname(Base.find_package("PowerSystems")), ".."))
include(joinpath(BASE_DIR, "test", "data_5bus_pu.jl")) #.jl file containing 5-bus system data
nodes_5 = nodes5() # function to create 5-bus buses

sys = System(
    100.0,
    nodes_5,
    vcat(thermal_generators5(nodes_5), renewable_generators5(nodes_5)),
    loads5(nodes_5),
    branches5(nodes_5),
)

@show get_component(Bus, sys, "nodeA")
@show get_component(Line, sys, "1")

get_components(Bus, sys) |> collect

get_components(Branch, sys) |> collect

bus1 = get_component(Bus, sys, "nodeA")
@show get_name(bus1);
@show get_magnitude(bus1);

loads = collect(get_components(PowerLoad, sys))
for (l, ts) in zip(loads, load_timeseries_DA[2])
    add_time_series!(
        sys,
        l,
        Deterministic(
            "activepower",
            Dict(TimeSeries.timestamp(load_timeseries_DA[2][1])[1] => ts),
        ),
    )
end

get_forecast_initial_times(sys)

ts_names = get_time_series_names(Deterministic, loads[1])

ta = get_time_series_array(Deterministic, loads[1], ts_names[1])

ts = get_time_series_values(Deterministic, loads[1], ts_names[1])

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# using Combinatorics
# Core data structure for the state of an SRA analysis
#   See also: sra_model_plotting.jl for plotting methods
#   See also: sra_model_function.jl for the analysis steps
# include("../utils/sindy_statistics_utils.jl")

#####
##### Helper structure, for the analysis parameters
#####
mutable struct sra_parameters
    debug_mode::Bool
    debug_folder::String
    # Initial, naive model
    initial_subsampling::Bool
    # Calculating residuals of the previous model
    use_turing::Bool
    noise_factor::Number
    # Parameters shared between all iterations of SINDy model
    # variable_names
    # sindy_library::Dict
    model_template::DynamicalSystemModel
    sindy_terms_list
    sindyc_ensemble_parameters
    # Subsampling to get the next iteration
    num_pts::Number
    start_ind::Number
end

# Initializer with defaults
function get_sra_defaults(use_control=true)
    debug_mode = true;
    debug_folder = "";
    # Initial, naive model
    initial_subsampling = false;
    # Calculating residuals of the previous model
    use_turing = false;
    noise_factor = 1.0;
    # Parameters shared between all iterations of SINDy model
    variable_names = ["x", "y", "z"]
    sindy_library = Dict("cross_terms"=>[2],"constant"=>nothing);
    optimizer = slstQuantile(0.1);
    if use_control
        sindyc_ensemble_parameters = get_sindyc_ensemble_parameters();
        model_template = sindycModel(
            convert_string2function(sindy_library), variable_names, optimizer);
    else
        sindyc_ensemble_parameters = Dict();
        model_template = sindyModel(
            convert_string2function(sindy_library), variable_names, optimizer);
    end
    sindy_terms_list = Iterators.product(1:3,1:3,1:3)
    # Subsampling to get the next iteration
    num_pts = 1000;
    start_ind = 1;

    return sra_parameters(debug_mode, debug_folder,
        initial_subsampling,
        use_turing,
        noise_factor,
        model_template,
        sindy_terms_list,
        sindyc_ensemble_parameters,
        num_pts,
        start_ind)
end

#####
##### Helper structure, for the true values
#####
struct sra_truth_object
    dat_raw::Array # Uncontrolled data
    true_grad::Array
    U_true
    core_dyn_true
end

#####
##### Main structure, for the current state
#####
mutable struct sra_stateful_object
    is_saved::Bool
    parameters::sra_parameters
    i::Number
    # Initial data and gradients
    ts
    tspan
    dat::Array
    u0::Vector # Note: initial conditions may not be same as training data
    numerical_grad::Array
    # Analysis history
    all_sindy_models
    all_ensemble_sindy_criteria
    all_best_sindy_criteria
    all_sindy_dat
    all_residuals
    all_subsample_ind
    all_noise_guesses
    # Current analysis state
    sindy_model
    ensemble_sindy_criteria
    best_sindy_criteria
    sindy_dat
    residual
    subsample_ind
    noise_guess
end

# Initializer with defaults
function sra_stateful_object(ts, tspan, dat, u0, numerical_grad)
    return sra_stateful_object(
        false,
        get_sra_defaults(),
        0, # i.e. uninitialized
        ts, tspan, dat, u0, numerical_grad,
        [], [], [], [], [], [], [],
        nothing, nothing, nothing, nothing, nothing, nothing, nothing
        )
end

export sra_stateful_object, sra_truth_object,
    sra_parameters, get_sra_defaults