text:
@testset "PolyProducts" begin

@info("--------- Testing Products of Polynomials ----------")

##
using ACE, Test, ForwardDiff, JuLIP, JuLIP.Testing

using LinearAlgebra: norm, cond
using ACE.OrthPolys: OrthPolyBasis
using JuLIP: evaluate, evaluate_d

##

maxn = 30
tdf = rand(200)
ww = 1.0 .+ rand(200)
Jd = OrthPolyBasis(maxn, 2, 1.0, 1, -1.0, tdf, ww)

##

@info("Testing the Expansion of products of radial polynomials")
coeffs = ACE.OrthPolys.OrthPolyProdCoeffs(Jd)

for n1 = 1:5, n2 = 1:5
   # product
   f1(x) = evaluate(Jd, x)[n1] * evaluate(Jd, x)[n2]
   # expansion
   f2(x) = ( P = coeffs(n1, n2);
             sum( evaluate(Jd, x)[ν] * P[ν]  for ν = 1:length(P) ) )
   for ntest = 1:10
      x = 2 * rand() - 1
      print_tf(@test f1(x) ≈ f2(x))
   end
end

##

@info("Testing Expansion of Products of Spherical Harmonics")
using ACE.Orth: SHProdCoeffs
using ACE.SphericalHarmonics: SHBasis, index_y

coeffs = SHProdCoeffs()
for l1 = 0:4, l2=l1:4, m1 = -l1:l1, m2 = -l2:l2
   f1 = let SH = SHBasis(10)
      x -> ( Y = evaluate(SH, x);
             Y[index_y(l1, m1)] * Y[index_y(l2, m2)] )
   end
   f2 = let P = coeffs(l1, m1, l2, m2), SH = SHBasis(10)
      x -> sum( p * evaluate(SH, x)[index_y(L, M)]
                for (L, M, p) in P )
   end
   nfail = 0
   for ntest = 1:10
      x = randn(JVec); x /= norm(x)
      nfail += (abs(f1(x) - f2(x)) > 1e-12)
   end
   print_tf(@test nfail == 0)
end

##

end


text:
push!(LOAD_PATH, "./femmodule/")
push!(LOAD_PATH, "./gridmodule/")

using BaseFEM
using BaseGrid

# m=grid1D("./mesh1d/mesh3")
# m=grid2D("./mesh2d/test")
m=grid3D("./mesh3d/ball")

top  = m.topology
ndof = ndof_P1(top)
geo  = restrain_geometry(m.geometry,ndof)

A,M=assembly_P1(geo,top)

# solve with homogeneous dirichlet boundary conditions
reduce = topology_id(geo.topology,1).==0

# elimination of boundary conditions
reduce = reduce[1:ndof]
u      = zeros(ndof)
rhs    = M*ones(ndof)

# solve reduced system
A_reduced    = A[reduce,reduce]
rhs_reduced  = rhs[reduce]
u_reduced    = A_reduced\rhs_reduced

# extend to full system again
u[reduce] = u_reduced

# output to ParaView vtk file
outvtk_scalar("output/example_poisson.vtk",geo,top,"sol",u)

ext: jl

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
# mapping a distance to a bin id
function dist2bin(ℓ, r, diam, n_bins)
    floor(Int, min(1 - log(r, ℓ / diam), n_bins))
end

# A struct that contains the assignment of the different points to the different bins of length scales
struct Rack
    bins::Vector{Vector{Int}}
    r::Float64
    diam::Float64
    ℓ::Vector{Float64}
    addresses::Vector{Tuple{Int,Int}}
end

# Constructor for Rack
function Rack(n_bins, r, diam, ℓ)
    bins = Vector{Vector{Int}}(undef, n_bins)
    for i = 1 : length(bins)
        bins[i] = Vector{Int}(undef,0)
    end
    addresses = Vector{Tuple{Int,Int}}(undef, length(ℓ))
    for (i, l) in enumerate(ℓ)
        bin = dist2bin(l, r, diam, n_bins)
        push!(bins[bin], i)
        addresses[i] = (bin, length(bins[bin]))
    end
    return Rack(bins, r, diam, ℓ, addresses)
end

# function that modifies distance value. Will only reduce the distance value and otherwise not to anything
function reduce_distance!(rack::Rack, id, l)
    if l < rack.ℓ[id]
        # updating the distance
        rack.ℓ[id] = min(l, rack.ℓ[id])
        # finding the new bin
        new_bin = dist2bin(l, rack.r, rack.diam, length(rack.bins))
        # finding the old address (and therefore the old bin)
        old_address = rack.addresses[id]
        # only act, if the bin changes
        if new_bin > old_address[1]
            # replacing the old position of id with the last element of the bin
            rack.bins[old_address[1]][old_address[2]] = rack.bins[old_address[1]][end]
            # updating the address of the last element of the bin
            rack.addresses[rack.bins[old_address[1]][end]] = old_address
            # dropping the last element of the original bin
            pop!(rack.bins[old_address[1]])
            # add id to new bin
            push!(rack.bins[new_bin], id)
            # overwrite address of id
            rack.addresses[id] = (new_bin, length(rack.bins[new_bin]))
        end
    end
end

# looks for the next dof to include in the ordering, starting at the bin first_bin
# returns the new first_bin, where it found the dof, as well as the id of the dof.
# If all remaining dofs are in the last bin, it signals it by returning zero
function next_dof(first_bin::Int, rack::Rack)
    while first_bin < length(rack.bins)
        if !isempty(rack.bins[first_bin])
            return first_bin, first(rack.bins[first_bin])
        end
        first_bin += 1
    end
    return 0, 0
end


text:
using BHAtp, Statistics

ProjDir = @__DIR__
!isdir(joinpath(ProjDir, "plots")) && mkdir(joinpath(ProjDir, "plots"))

ProjName = split(ProjDir, "/")[end]
bhaj = BHAtp.BHAJ(ProjName, ProjDir)
bhaj.ratio = 0.5

segs = [
	# Element type,  Material,    Length,     OD,   ID,  fc
	(:bit,           :steel,        0.00,   2.75,  9.0,  0.0),
	(:collar,        :steel,       45.00,   2.75,  7.0,  0.0),
	(:stabilizer,    :steel,        0.00,   2.75,  9.0,  0.0),
	(:pipe,          :steel,       30.00,   2.75,  7.0,  0.0),
	(:stabilizer,    :steel,        0.00,   2.75,  9.0,  0.0),
	(:pipe,          :steel,       90.00,   2.75,  7.0,  0.0),
	(:stabilizer,    :steel,        0.00,   2.75,  9.0,  0.0),
	(:pipe,          :steel,      100.00,   2.75,  7.0,  0.0)
]

traj = [
	#Heading,  Diameter
	( 60.0,      9.0)
]

wobs = 20:10:40
incls = 20:10:50  # Or e.g. incls = [5 10 20 30 40 45 50]

@time bhaj(segs, traj, wobs, incls)

display("Fetch tp=false, wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false, tp=false)
df[:,[2; 5:6; 9:12]] |> display
println()

println("Fetch tp=true, wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false);
df[:,[2; 5:6; 9:12]] |> display
println()

tpf = create_final_tp_df(ProjDir, wobs, incls; ofu=false)
tpf |> display
println()

# Collect computed side forces on NBS @ 45ft

sfs = []
for wob in wobs
  for incl in incls
    df, df_tp = show_solution(ProjDir, wob, incl, show=false);
    sf = @linq df |> where(:dtb .== 45)
    append!(sfs, sf.fx)
  end
end

# Simulated actual tps

noise = 0.1rand(2*length(incls))
noise = noise .- mean(noise)
actual = DataFrame(
  incls = LinRange(15.0, 50.0, 2*length(incls)),
  atp = LinRange(-2.7, -4.5, 2*length(incls)) .+ noise
)

# Function to simulated formation correction factors

intp_fcf  = LinearInterpolation(LinRange(0.0, 3000, 10), LinRange(0.63, 0.72, 10));
ctp = tpf.tp .- intp_fcf.(sfs)

result_df = DataFrame(
  wob = sort(repeat(wobs, length(incls))),
  incl = repeat(incls, length(wobs)),
  tp = tpf.tp,
  side_force = sfs,
  corr_factor = intp_fcf.(sfs),
  corr_tp = ctp
);
result_df |> display
println()

name = "case_01a"

case_01a = @pgf GroupPlot({
  group_style = {
    group_size = "1 by 1",
    vertical_sep = "0pt",
    xticklabels_at = "edge bottom" }
  },
  Axis({ height = "6cm", width = "8cm", xlabel = "Inclination [°]",
      ylabel = "TP [°/100ft]", title = "Actual vs. Corrected TP" },
      Plot({ color="blue", mark = "*" }, Table(incls, result_df.tp[9:12])),
      LegendEntry("Isotropic TP"),
      Plot({ color="red", mark = "*" }, Table(incls, result_df.corr_tp[9:12])),
      LegendEntry("Predicted"),
      Plot({ thick, color="black", mark = "+" }, Table(actual.incls, actual.atp)),
      LegendEntry("Actual")))

savefig(joinpath(ProjDir, "plots", name), case_01a, ProjDir)

name = "case_01b"

case_01b = @pgf GroupPlot({
  group_style = {
    group_size = "1 by 3",
    vertical_sep = "0pt",
    xticklabels_at = "edge bottom" }
  },
  Axis({
        height = "6cm", width = "8cm",
        ylabel = "TP [°/100ft]", title = "TP and Formatioin Correction (ΔTP)"
      },
      Plot({ color="red", mark = "*" }, Table(incls, tpf.tp[1:4])),
      LegendEntry("WOB 20"),
      Plot({ color="blue", mark = "+" }, Table(incls, tpf.tp[5:8])),
      LegendEntry("WOB 30"),
      Plot({ color="green", mark = "o" }, Table(incls, tpf.tp[9:12])),
      LegendEntry("WOB 40")),
  Axis({
        height = "6cm", width = "8cm",
        ylabel = "Sideforce NBS"
      },
      Plot({ color="red", mark = "*", title ="TP" }, Table(incls, sfs[1:4])),
      LegendEntry("WOB 20"),
      Plot({ color="blue", mark = "+" }, Table(incls, sfs[5:8])),
      LegendEntry("WOB 30"),
      Plot({ color="green", mark = "o" }, Table(incls, sfs[9:12])),
      LegendEntry("WOB 40")),
   Axis({ height = "6cm", width = "8cm", xlabel = "Inclination [°]",
      ylabel = "ΔTP [°/100ft]" },
      Plot({ color="red", mark = "*" }, Table(incls, result_df.corr_factor[1:4])),
      LegendEntry("WOB 20"),
      Plot({ color="blue", mark = "+" }, Table(incls, result_df.corr_factor[5:8])),
      LegendEntry("WOB 30"),
      Plot({ color="green", mark = "o" }, Table(incls, result_df.corr_factor[9:12])),
      LegendEntry("WOB 40")))

savefig(joinpath(ProjDir, "plots", name), case_01b, ProjDir)


text:
################################################################################
#
#   HELPER FUNCTIONS FOR CONSTRUCTION OF LATTICES / NEIGHBORS
#
#   FILE CONTAINS
#   - bond concanation
#
################################################################################

################################################################################
#
#   CONCANATION OF BONDS
#
################################################################################

# N=0
function concanateBonds(
            b1 :: B,
            b2 :: B,
            labelfunction :: Function = (bo1,bo2) -> getDefaultLabel(L)
        ) :: B where {L,B<:AbstractBond{L,0}}

    # assert that the bonds can indeed by concanated
    @assert to(b1) == from(b2) "Bonds cannot be concanated as "*string(from(b1))*"-"*string(to(b1))*" (-) "*string(from(b2))*"-"*string(to(b2))*" not possible"

    # create a new bond
    return newBond(
        B,
        from(b1) :: Int64,
        to(b2) :: Int64,
        labelfunction(b1,b2) :: L,
        NTuple{0,Int64}()
    )
end
# N=1
function concanateBonds(
            b1 :: B,
            b2 :: B,
            labelfunction :: Function = (bo1,bo2) -> getDefaultLabel(L)
        ) :: B where {L,B<:AbstractBond{L,1}}

    # assert that the bonds can indeed by concanated
    @assert to(b1) == from(b2) "Bonds cannot be concanated as "*string(from(b1))*"-"*string(to(b1))*" (-) "*string(from(b2))*"-"*string(to(b2))*" not possible"

    # create a new bond
    return newBond(
        B,
        from(b1) :: Int64,
        to(b2) :: Int64,
        labelfunction(b1,b2) :: L,
        (
            wrap(b1)[1] + wrap(b2)[1],
        )
    )
end
# N=2
function concanateBonds(
            b1 :: B,
            b2 :: B,
            labelfunction :: Function = (bo1,bo2) -> getDefaultLabel(L)
        ) :: B where {L,B<:AbstractBond{L,2}}

    # assert that the bonds can indeed by concanated
    @assert to(b1) == from(b2) "Bonds cannot be concanated as "*string(from(b1))*"-"*string(to(b1))*" (-) "*string(from(b2))*"-"*string(to(b2))*" not possible"

    # create a new bond
    return newBond(
        B,
        from(b1) :: Int64,
        to(b2) :: Int64,
        labelfunction(b1,b2) :: L,
        (
            wrap(b1)[1] + wrap(b2)[1],
            wrap(b1)[2] + wrap(b2)[2]
        )
    )
end
# N=3
function concanateBonds(
            b1 :: B,
            b2 :: B,
            labelfunction :: Function = (bo1,bo2) -> getDefaultLabel(L)
        ) :: B where {L,B<:AbstractBond{L,3}}

    # assert that the bonds can indeed by concanated
    @assert to(b1) == from(b2) "Bonds cannot be concanated as "*string(from(b1))*"-"*string(to(b1))*" (-) "*string(from(b2))*"-"*string(to(b2))*" not possible"

    # create a new bond
    return newBond(
        B,
        from(b1) :: Int64,
        to(b2) :: Int64,
        labelfunction(b1,b2) :: L,
        (
            wrap(b1)[1] + wrap(b2)[1],
            wrap(b1)[2] + wrap(b2)[2],
            wrap(b1)[3] + wrap(b2)[3]
        )
    )
end

# export the concanation function
export concanateBonds


text:
module TestCommon

using StructTypes
using SimplePNGs
using JSON3
using ColorTypes
using FixedPointNumbers
using CodecZlib

mutable struct JsonIMG
    T::String
    width::Int
    height::Int
    data::Vector{UInt16}

    JsonIMG() = new()
    JsonIMG(T, w, h, data) = new(T, w, h, data)
end
StructTypes.StructType(::Type{JsonIMG}) = StructTypes.Mutable()

function JsonIMG(x::AbstractArray{T0}) where {T0 <: AbstractGray}
    T = eltype(x) |> Symbol |> String
    T = replace(T, "ColorTypes." => "")
    T = replace(T, "FixedPointNumbers." => "")
    width, height = size(x)
    data = vec(reinterpret.(gray.(x)))

    return JsonIMG(T, width, height, data)
end

function JsonIMG(x::AbstractArray{T0}) where {T0 <: GrayA}
    T = eltype(x) |> Symbol |> String
    T = replace(T, "ColorTypes." => "")
    T = replace(T, "FixedPointNumbers." => "")
    width, height = size(x)
    data = vcat(vec(reinterpret.(gray.(x))),
                vec(reinterpret.(alpha.(x))))

    return JsonIMG(T, width, height, data)
end

function JsonIMG(x::AbstractArray{T0}) where {T0 <: AbstractRGB}
    T = eltype(x) |> Symbol |> String
    T = replace(T, "ColorTypes." => "")
    T = replace(T, "FixedPointNumbers." => "")
    width, height = size(x)
    data = vcat(vec(reinterpret.(red.(x))),
                vec(reinterpret.(green.(x))),
                vec(reinterpret.(blue.(x))))

    return JsonIMG(T, width, height, data)
end

function JsonIMG(x::AbstractArray{T0}) where {T0 <: RGBA}
    T = eltype(x) |> Symbol |> String
    T = replace(T, "ColorTypes." => "")
    T = replace(T, "FixedPointNumbers." => "")
    width, height = size(x)
    data = vcat(vec(reinterpret.(red.(x))),
                vec(reinterpret.(green.(x))),
                vec(reinterpret.(blue.(x))),
                vec(reinterpret.(alpha.(x))))

    return JsonIMG(T, width, height, data)
end

function reconstruct(x::JsonIMG)
    T0 = replace(x.T, "ColorTypes." => "")
    T0 = replace(T0, "FixedPointNumbers." => "")
    res = if T0 == "Gray{Normed{UInt8,8}}"
        T = Gray{N0f8}
        data = UInt8.(x.data)
        SimplePNGs.pixel.(T, data)
    elseif T0 == "RGB{Normed{UInt8,8}}"
        T = RGB{N0f8}
        lng = x.width * x.height
        data = UInt8.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng], data[2*lng+1:3*lng])
    elseif T0 == "Gray{Normed{UInt16,16}}"
        T = Gray{N0f16}
        data = UInt16.(x.data)
        SimplePNGs.pixel.(T, data)
    elseif T0 == "RGB{Normed{UInt16,16}}"
        T = RGB{N0f16}
        lng = x.width * x.height
        data = UInt16.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng], data[2*lng+1:3*lng])
    elseif T0 == "GrayA{Normed{UInt8,8}}"
        T = GrayA{N0f8}
        lng = x.width * x.height
        data = UInt8.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng])
    elseif T0 == "GrayA{Normed{UInt16,16}}"
        T = GrayA{N0f16}
        lng = x.width * x.height
        data = UInt16.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng])
    elseif T0 == "RGBA{Normed{UInt8,8}}"
        T = RGBA{N0f8}
        lng = x.width * x.height
        data = UInt8.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng], data[2*lng+1:3*lng], data[3*lng+1:4*lng])
    elseif T0 == "RGBA{Normed{UInt16,16}}"
        T = RGBA{N0f16}
        lng = x.width * x.height
        data = UInt16.(x.data)
        SimplePNGs.pixel.(T, data[1:lng], data[lng+1:2*lng], data[2*lng+1:3*lng], data[3*lng+1:4*lng])
    end

    reshape(res, x.width, x.height)
end

function save_json(name)
    namein = joinpath("..", "test", "PngSuite", name*".png")
    nameout = joinpath("..", "test", "jsons", name*"_json.gz")

    img = SimplePNGs.load(namein)
    out = JSON3.write(JsonIMG(img))
    compressed = transcode(GzipCompressor, out)
    write(nameout, compressed)
    nothing
end

function load_json(name, test = true)
    if test
        namein = joinpath("jsons", name*"_json.gz")
    else
        namein = joinpath("..", "test", "jsons", name*"_json.gz")
    end
    json = transcode(GzipDecompressor, read(namein))
    json_img = JSON3.read(json, JsonIMG)
    return reconstruct(json_img)
end

end # module


text:
using RecipesBase

@testset "Plot recipes" begin
    data = randn(10, 10)
    img = AstroImage(data)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img)
    @test getfield(rec[1], 1) == Dict{Symbol, Any}(:seriestype   => :heatmap,
                                                   :aspect_ratio => :equal,
                                                   :color        => :grays)
    @test rec[1].args == (img.data,)
end


text:
using Statistics

abstract type AbstractFiniteDifferencePolicySearch <: AbstractPolicyAgent end

# finite differences policy gradient search
mutable struct BasicRandomSearch{T,TP} <: AbstractFiniteDifferencePolicySearch where{T<:Real,TP<:AbstractPolicy}
    π::TP
    α::T
    γ::T
    ν::T
    num_samples::Int
    max_episodes
    θ::Array{T}
    δs::Array{T}
    Gs::Array{T}
    ep_number::Int
    dnum::Int
    pnum::Int
    G::T
    cur_gamma::T

    function BasicRandomSearch(π::AbstractPolicy, α::T, γ::T, ν::T, num_samples::Int, rng::AbstractRNG, max_episodes=Inf) where {T<:Real}
        θ = deepcopy(get_params(π))
        δs = ν .* randn(rng, (num_samples, size(θ)...))
        Gs = zeros((num_samples, 2))
        dnum = 1
        ep_number = 1
        pnum = 1
        new{T,typeof(π)}(clone(π), α, γ, ν, num_samples, max_episodes, θ, δs, Gs, ep_number, dnum, pnum, 0., 1.)
    end

    function BasicRandomSearch(π::AbstractPolicy, α::T, γ::T, ν::T, num_samples::Int, max_episodes=Inf) where {T<:Real}
        θ = deepcopy(get_params(π))
        δs = zeros(T, (num_samples, size(θ)...))
        Gs = zeros(T, (num_samples, 2))
        dnum = 1
        ep_number = 1
        pnum = 1
        new{T,typeof(π)}(clone(π), α, γ, ν, num_samples, max_episodes, θ, δs, Gs, ep_number, dnum, pnum, 0., 1.)
    end
end

function act!(agent::AbstractFiniteDifferencePolicySearch, env::AbstractEnvironment, rng::AbstractRNG)
    get_action!(agent.π, env.state, rng)

    reward = step!(env, agent.π.action, rng)

    agent.G += agent.cur_gamma * reward
    agent.cur_gamma *= agent.γ

    if is_terminal(env)
        agent.Gs[agent.dnum, agent.pnum] = agent.G
        agent.G = 0
        agent.ep_number += 1
        if agent.pnum == 1
            agent.pnum = 2
        else
            agent.pnum = 1
            agent.dnum += 1
        end
    end
end

function new_episode!(agent::AbstractFiniteDifferencePolicySearch, rng::AbstractRNG)
    agent.G = 0.
    agent.cur_gamma = 1.0

    set_params!(agent.π, agent.θ)

    if agent.ep_number ≤ agent.max_episodes
        if agent.dnum > agent.num_samples
            update_theta!(agent)
            agent.δs .= agent.ν .* randn(rng, (agent.num_samples, size(agent.θ)...))
            agent.dnum = 1
            agent.pnum = 1
        end

        Δ = view(reshape(agent.δs, (agent.num_samples, :)), agent.dnum, :)
        if agent.pnum == 1
            add_to_params!(agent.π, -Δ)
        else
            add_to_params!(agent.π, Δ)
        end
    end

end

function update_theta!(agent::BasicRandomSearch)
    θ = agent.θ
    Gs = agent.Gs
    δs = agent.δs

    g = (view(Gs, :, 2) .- view(Gs, :, 1)) .* reshape(δs, (agent.num_samples, :))
    ḡ =  reshape(mean(g, dims=1), size(θ))
    @. agent.θ += agent.α * ḡ
end

function clone(agent::BasicRandomSearch)::BasicRandomSearch
    a = BasicRandomSearch(agent.π, agent.α, agent.γ, agent.ν, agent.num_samples, agent.max_episodes)
    a.θ .= agent.θ
    a.δs .= agent.δs
    a.Gs .= agent.Gs
    a.ep_number = agent.ep_number
    a.dnum = agent.dnum
    a.pnum = agent.pnum
    a.G = agent.G
    a.cur_gamma = agent.cur_gamma

    return a
end

mutable struct AugmentedRandomSearch{T,TP} <: AbstractFiniteDifferencePolicySearch where{T<:Real,TP<:AbstractPolicy}
    π::TP
    α::T
    γ::T
    ν::T
    num_samples::Int
    num_elite::Int
    max_episodes
    θ::Array{T}
    δs::Array{T}
    Gs::Array{T}
    ep_number::Int
    dnum::Int
    pnum::Int
    G::T
    cur_gamma::T

    function AugmentedRandomSearch(π::AbstractPolicy, α::T, γ::T, ν::T, num_samples::Int, num_elite::Int, rng::AbstractRNG, max_episodes=Inf) where {T<:Real}
        θ = deepcopy(get_params(π))
        δs = ν .* randn(rng, (num_samples, size(θ)...))
        Gs = zeros((num_samples, 2))
        dnum = 1
        ep_number = 1
        pnum = 1
        new{T,typeof(π)}(clone(π), α, γ, ν, num_samples, num_elite, max_episodes, θ, δs, Gs, ep_number, dnum, pnum, 0., 1.)
    end

    function AugmentedRandomSearch(π::AbstractPolicy, α::T, γ::T, ν::T, num_samples::Int, num_elite::Int, max_episodes=Inf) where {T<:Real}
        θ = deepcopy(get_params(π))
        δs = zeros(T, (num_samples, size(θ)...))
        Gs = zeros(T, (num_samples, 2))
        dnum = 1
        ep_number = 1
        pnum = 1
        new{T,typeof(π)}(clone(π), α, γ, ν, num_samples, num_elite, max_episodes, θ, δs, Gs, ep_number, dnum, pnum, 0., 1.)
    end
end

function update_theta!(agent::AugmentedRandomSearch)
    θ = agent.θ
    Gs = agent.Gs
    δs = agent.δs

    idxs = sortperm(vec(maximum(Gs, dims=2)), rev=true)[1:agent.num_elite]

    g  = (view(Gs, idxs, 2) .- view(Gs, idxs, 1)) .* view(reshape(δs, (agent.num_samples, :)), idxs, :)
    ḡ  = mean(g, dims=1)
    ḡ  = reshape(ḡ, size(θ))
    σᵣ = std(view(Gs, idxs, :))  # standard deviation of the 2*num_elite returns used
    @. agent.θ += agent.α / (σᵣ + 1e-8) * ḡ
end

function clone(agent::AugmentedRandomSearch)::AugmentedRandomSearch
    a = AugmentedRandomSearch(agent.π, agent.α, agent.γ, agent.ν, agent.num_samples, agent.num_elite, agent.max_episodes)
    a.θ .= agent.θ
    a.δs .= agent.δs
    a.Gs .= agent.Gs
    a.ep_number = agent.ep_number
    a.dnum = agent.dnum
    a.pnum = agent.pnum
    a.G = agent.G
    a.cur_gamma = agent.cur_gamma

    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using YaoLocations
using BenchmarkTools
using BQCESubroutine
using StrideArrays
using LoopVectorization
using ThreadingUtilities
using ArrayInterface
#using CheapThreads      # @batch is not working properly. see https://github.com/JuliaSIMD/Polyester.jl/issues/32
using BQCESubroutine: BitSubspace, bsubspace, bcomspace

function subspace_mul4x4!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        for b in 1:size(st, 1)
            Base.Cartesian.@nexprs 4 i -> begin
                y_i = zero(T)
                Base.Cartesian.@nexprs 4 j -> begin
                    y_i += U_i_j * st[b, idx_j]
                end
            end

            Base.Cartesian.@nexprs 4 i -> begin
                st[b, idx_i] = y_i
            end
        end
    end
    return st
end

@inline function subspace_mul4x4_task!(st::AbstractMatrix{T}, indices, U, subspace, range, offset) where T
    D = StrideArrays.static_length(indices)
    y = StrideArray{T}(undef, (D, ))

    @inbounds begin
        Base.Cartesian.@nexprs 4 i -> begin
            indices_i = indices[i]
        end

        Base.Cartesian.@nexprs 4 i -> begin
            Base.Cartesian.@nexprs 4 j -> begin
                U_i_j = U[i][j]
            end
        end
    end

    @inbounds for (s, b) in range
        k = subspace[s]
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        Base.Cartesian.@nexprs 4 i -> begin
            y[i] = zero(T)
            Base.Cartesian.@nexprs 4 j -> begin
                y[i] += U_i_j * st[b, idx_j]
            end
        end

        Base.Cartesian.@nexprs 4 i -> begin
            st[b, idx_i] = y[i]
        end
    end
    return
end

function threaded_subspace_mul4x4!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    space = BQCESubroutine.CartesianSpace(1:length(subspace), 1:size(st, 1))
    nthreads = Threads.nthreads()
    U = (U[1, 1], U[1, 2], U[1, 3], U[1, 4]),
        (U[2, 1], U[2, 2], U[2, 3], U[2, 4]),
        (U[3, 1], U[3, 2], U[3, 3], U[3, 4]),
        (U[4, 1], U[4, 2], U[4, 3], U[4, 4])

    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1
    indices = (indices_1, indices_2, indices_3, indices_4)

    #@batch for tid in 1:nthreads
    for tid in 1:nthreads
        range = BQCESubroutine.schedule_task(space, tid, nthreads)
        subspace_mul4x4_task!(st, indices, U, subspace, range, offset)
    end
    return st
end

T = Float64
n = 10
S = rand(ComplexF64, 100, 1<<n);
U = rand(ComplexF64, 1<<2, 1<<2);
locs = Locations((1, 3))

subspace = bsubspace(n, locs)
comspace = bcomspace(n, locs)
inds = StrideArray{Int}(undef, (length(comspace), ))
@simd ivdep for i in eachindex(inds)
    inds[i] = comspace[i] + 1
end

S1 = BQCESubroutine.subspace_mul_generic!(copy(S), inds, U, subspace);
S2 = BQCESubroutine.threaded_subspace_mul_generic!(copy(S), inds, U, subspace)
S3 = subspace_mul4x4!(copy(S), comspace, U, subspace)
S4 = threaded_subspace_mul4x4!(copy(S), comspace, U, subspace)

S1 ≈ S2
S1 ≈ S3
S1 ≈ S4

@benchmark BQCESubroutine.subspace_mul_generic!($(copy(S)), inds, $U, subspace)
@benchmark subspace_mul4x4!($(copy(S)), comspace, $U, subspace)
@benchmark BQCESubroutine.threaded_subspace_mul_generic!($(copy(S)), inds, $U, subspace)
@benchmark threaded_subspace_mul4x4!($(copy(S)), comspace, $U, subspace)

n = 15
S = rand(ComplexF64, 1<<n);
U = rand(ComplexF64, 1<<3, 1<<3);
locs = Locations((1, 3, 5))
subspace = bsubspace(n, locs)
comspace = bcomspace(n, locs)
indices = StrideArray{Int}(undef, (StaticInt{length(comspace)}(), ))
@simd ivdep for i in eachindex(indices)
    indices[i] = comspace[i] + 1
end

S2 = BQCESubroutine.subspace_mul_generic!(copy(S), indices, U, subspace)
S3 = BQCESubroutine.threaded_subspace_mul_generic!(copy(S), indices, U, subspace)
S2 ≈ S3

using YaoArrayRegister

@benchmark instruct!($(copy(S)), $U, (1, 3, 5))

@benchmark BQCESubroutine.subspace_mul_generic!($(copy(S)), indices, U, subspace)
@benchmark BQCESubroutine.threaded_subspace_mul_generic!($(copy(S)), indices, U, subspace)

@benchmark BQCESubroutine.subspace_mul_generic!(S, $indices, $U, $subspace) setup=(S=$(copy(S)))
@benchmark BQCESubroutine.threaded_subspace_mul_generic!(S, $indices, $U, $subspace) setup=(S=$(copy(S))) samples=1000 evals=1
t = @benchmark threaded_subspace_mm!(S, $indices, $U, $subspace) setup=(S=$(copy(S))) samples=1000 evals=1

@time threaded_subspace_mm!(S, indices, U, subspace);
BQCESubroutine.subspace_mul_generic!(S, indices, U, subspace);

D = ArrayInterface.static_length(indices)
U_re = StrideArray{T}(undef, (D, D))
U_im = StrideArray{T}(undef, (D, D))
@inbounds @simd ivdep for i in 1:length(U)
    U_re[i] = real(U[i])
    U_im[i] = imag(U[i])
end

subspace_ref = Ref(subspace)
Bmax = size(S, 1)
total = length(subspace) * (((Bmax-1)>>>3)+1)
nthreads = 6
len, rem = divrem(total, nthreads)
S4 = copy(S)
Sr = reinterpret(reshape, T, S4)
offset = 0
range1 = BQCESubroutine.div_thread(1, len, rem)
range2 = BQCESubroutine.div_thread(2, len, rem)
range3 = BQCESubroutine.div_thread(3, len, rem)
range4 = BQCESubroutine.div_thread(4, len, rem)
range5 = BQCESubroutine.div_thread(5, len, rem)
range6 = BQCESubroutine.div_thread(6, len, rem)

BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range1, offset)
BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range2, offset)
BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range3, offset)
BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range4, offset)
BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range5, offset)
BQCESubroutine.subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range6, offset)

p = ThreadingUtilities.taskpointer(1)
BQCESubroutine.setup_subspace_mul_generic(p, S, indices, U_re, U_im, subspace_ref, range1, 0)

f = BQCESubroutine.SubspaceMulGeneric{T, D}()
f(p)

BQCESubroutine.launch_subspace_mul_generic(1, S, indices, U_re, U_im, subspace_ref, range2, offset)
ThreadingUtilities.wait(1)

GC.@preserve S U_re U_im subspace_ref begin
    for tid in 1:nthreads-1
        range = div_thread(tid, len, rem)
        launch_subspace_mul_generic(tid, S, indices, U_re, U_im, subspace_ref, range, offset)
    end

    range = div_thread(nthreads, len, rem)
    subspace_mul_generic_task!(Sr, indices, U_re, U_im, subspace, range, offset)

    for tid in 1:nthreads-1
        ThreadingUtilities.wait(tid)
    end
end


text:
using Documenter

makedocs(
    sitename = "Wandb",
    authors = "Avik Pal",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages = [
        "Home"              => "index.md",
        "QuickStart"        => "quickstart.md",
        "Miscellaneous"     => "misc.md",
        "Examples" => [
            "Getting Started"             => "examples/demo.md",
            "Flux.jl Intergration"        => "examples/flux.md",
            "FluxTraining.jl Integration" => "examples/fluxtraining.md",
            "HyperParameter Sweeps"       => "examples/hparams.md",
            "Artifacts API"               => "examples/artifacts.md",
            "MPI.jl Integration"          => "examples/mpi.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/avik-pal/Wandb.jl.git",
    push_preview = true,
    devbranch = "main",
)


text:
function DenseOutputModel()
    return DenseOutputModel(())
end

function append(obj::DenseOutputModel, arg0::DenseOutputModel)
    return jcall(obj, "append", void, (DenseOutputModel,), arg0)
end

function get_final_time(obj::DenseOutputModel)
    return jcall(obj, "getFinalTime", jdouble, ())
end

function get_initial_time(obj::DenseOutputModel)
    return jcall(obj, "getInitialTime", jdouble, ())
end

function get_interpolated_state(obj::DenseOutputModel, arg0::jdouble)
    return jcall(obj, "getInterpolatedState", ODEStateAndDerivative, (jdouble,), arg0)
end

function handle_step(obj::DenseOutputModel, arg0::ODEStateInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (ODEStateInterpolator, jboolean), arg0, arg1)
end

function init(obj::DenseOutputModel, arg0::ODEStateAndDerivative, arg1::jdouble)
    return jcall(obj, "init", void, (ODEStateAndDerivative, jdouble), arg0, arg1)
end


text:
# policy file better sanity fewer contacts
pinft = 0.45.*pinf_o

pshopt = [[1,1,1],[1,1,1]] - alpha_open*[[0.1,0.33,0],[0.1,0.33,0]]

shopmeett = [8,28,10]

tauadj = tauadj_o

socialmaxyyt = socialmaxyy_o
#socialmaxoyt = 2
#socialmaxoot = 1
#socialmaxyot = 1


text:
module SparseIterateTest

using Test
using ProximalBase
using Random
using LinearAlgebra
using SparseArrays

@testset "sparseIterate" begin

    @testset "issparse" begin
        @test issparse(SparseIterate(sparse(fill(1,5,5))))
        @test issparse(SparseIterate(fill(1,5,5)))
        @test issparse(SparseIterate(fill(1,5)))
    end

    @testset "dropzeros and nnz" begin
        x = SparseIterate(3)
        x[3] = 2.
        x[1] = 1.
        @test nnz(x) == 2
        x[1] = 0.
        @test nnz(x) == 2
        dropzeros!(x)
        @test nnz(x) == 1
        @test Vector(x) == [0., 0., 2.]

        x = SparseIterate(3)
        x[3] = 2.
        x[1] = 1.
        x[3] = 0.
        x[1] = 0.
        @test nnz(x) == 2
        dropzeros!(x)
        @test nnz(x) == 0
        @test Vector(x) == [0., 0., 0.]
    end

    @testset "show" begin
        io = IOBuffer()
        x = SparseIterate(3)

        show(io, MIME"text/plain"(), x)
        @test String(take!(io)) == "3-element ProximalBase.SparseIterate{Float64,1} with 0 stored entries"

        x[1] = 1.
        show(io, MIME"text/plain"(), x)
        @test String(take!(io)) == "3-element ProximalBase.SparseIterate{Float64,1} with 1 stored entry:\n  [1]  =  1.0\n"

        x[2] = 2.
        show(io, MIME"text/plain"(), x)
        @test String(take!(io)) == "3-element ProximalBase.SparseIterate{Float64,1} with 2 stored entries:\n  [1]  =  1.0\n  [2]  =  2.0\n"

        show(io, x)
        @test String(take!(io)) == "  [1]  =  1.0\n  [2]  =  2.0\n"

        y = SparseIterate(2,3)

        show(io, MIME"text/plain"(), y)
        @test String(take!(io)) == "2×3 ProximalBase.SparseIterate{Float64,2} with 0 stored entries"
        y[1,3] = 1.
        show(io, MIME"text/plain"(), y)
        @test String(take!(io)) == "2×3 ProximalBase.SparseIterate{Float64,2} with 1 stored entry:\n  [1, 3]  =  1.0"

        y[2,2] = 2.
        show(io, MIME"text/plain"(), y)
        @test String(take!(io)) == "2×3 ProximalBase.SparseIterate{Float64,2} with 2 stored entries:\n  [2, 2]  =  2.0\n  [1, 3]  =  1.0"
    end

    @testset "copy" begin
        v = sprand(20, 0.4)
        x = SparseIterate(v)
        y = copy(x)
        @test Vector(y) == Vector(v)
        @test Vector(y) == Vector(x)
        @test !(x === y)

        M = sprand(10, 20, 0.2)
        x = SparseIterate(M)
        y = copy(x)
        @test Matrix(y) == Matrix(M)
        @test Matrix(y) == Matrix(x)
        @test !(x === y)
    end

    @testset "other" begin
        x = SparseIterate(3)
        x[3] = 2.
        @test numCoordinates(x) == 3
        @test nnz(x) == 1
        @test convert(Vector, x) == [0., 0., 2.]
        @test convert(Array, x) == [0., 0., 2.]
        @test Vector(x) == [0., 0., 2.]

        x = SparseIterate(10)
        xv = randn(5)
        for i=1:5
            x[i] = xv[i]
        end
        y = randn(10)
        @test dot(xv, y[1:5]) ≈ dot(x, y)
        @test dot(y[1:5], xv) ≈ dot(x, y)

        # convert
        n, m = 10, 100
        M = sprandn(n, m, 0.2)
        x = SparseIterate(M)
        @test numCoordinates(x) == n*m
        @test Matrix(x) == Matrix(M)
        @test Array(x) == Matrix(M)

        # constructor
        @test typeof(SparseIterate(2, 3)) == SparseIterate{Float64, 2}

        # multiplication
        v = sprandn(50, 0.2)
        M = randn(30, 50)
        x = SparseIterate(v)

        @test mul!(zeros(30), M, x) == M*v

    end

end

@testset "sparseSymmetricIterate" begin

    @testset "dropzeros and nnz" begin
        x = SymmetricSparseIterate(2)
        x[3] = 2.
        x[1] = 1.
        @test nnz(x) == 2
        x[1] = 0.
        @test nnz(x) == 2
        dropzeros!(x)
        @test nnz(x) == 1
        @test Matrix(x) == [0.  0.; 0. 2.]

        x = SymmetricSparseIterate(2)
        x[3] = 2.
        x[1] = 1.
        x[3] = 0.
        x[1] = 0.
        @test nnz(x) == 2
        dropzeros!(x)
        @test nnz(x) == 0
        @test Matrix(x) == [0. 0.; 0. 0.]
    end

    @testset "other" begin
      # convert
      n = 10
      M = sprandn(n, n, 0.1)
      M = (M + M') / 2
      x = SymmetricSparseIterate(M)
      @test numCoordinates(x) == n * (n+1) / 2
      @test Matrix(x) == Matrix(M)
      @test Array(x) == Matrix(M)

      # constructor
      @test typeof(SymmetricSparseIterate(2)) == SymmetricSparseIterate{Float64}
    end

    @testset "mul" begin
      n = 10
      M = sprandn(n, n, 0.1)
      M = (M + M') / 2
      x = SymmetricSparseIterate(M)

      B = randn(10, 20)
      C = M * B
      CC = zeros(10, 20)
      @test mul!(CC, x, B) == C
    end

end

@testset "AtomIterate" begin
    x = AtomIterate((2,3), 2, true)

    @test length(x) == 6
    @test size(x) == (2,3)
    @test IndexStyle(typeof(x)) == IndexLinear()

    x[:] = collect(1:6)
    y = [1. 3 5; 2 4 6]

    @test x.storage == y

    g = ProxL2(1., [10., 10.])
    out = prox(g, x)
    @test out.storage == zeros(2,3)
    out = prox(g, x, 1.)
    @test out.storage == zeros(2,3)

    g = ProxL2(1., [1., 10.])
    out = prox(g, x)
    @test out.storage[1,:] ≈ (1. - 1. / norm(y[1,:]))*y[1,:]
    @test out.storage[2,:] == zeros(3)

    g = ProxL2(1., [1., 1.])
    out = prox(g, x)
    @test out.storage[1,:] ≈ (1. - 1. / norm(y[1,:]))*y[1,:]
    @test out.storage[2,:] ≈ (1. - 1. / norm(y[2,:]))*y[2,:]

    @test_throws ArgumentError AtomIterate((2,3), 2, false)

    x = AtomIterate((2,3), 3, false)
    x[:] = collect(1:6)

    g = ProxL2(1., [1., 1.])
    @test_throws DimensionMismatch prox(g, x)

    g = ProxL2(1., [0., 1., 0.])
    out = prox(g, x)
    @test out.storage[:, 1] == y[:,1]
    @test out.storage[:, 2] ≈ (1. - 1. / norm(y[:,2]))*y[:,2]
    @test out.storage[:, 3] == y[:,3]
end

end


text:
module SVD

const documents = begin
"""
the svd algorithm
suppose we have a datasets with M length:

.. math ::

   \{ (x_i, y_i) | i = 1, 2, 3, \cdots, M \}

where :math:`x_i` has :math:`N` features.

original data :

.. math ::

    & A \triangleq [x_1^T, x_2^T, \cdots, x_M^T ]^T \\
    & A = U_{M \times M} \Sigma_{M \times N} V_{N \times N}

where both :math:`U, V` are orthogonal and :math:`\Sigma_{M \times N}` is diagonal.

So, :math:`A A^T` = :math:`\Sigma \Sigma^T` leads to

.. math ::

    & A A^T U = \Sigma \Sigma^T U \\
    & \forall i = 1, 2, \cdots, M, \;\; A^T A u_i = \sigma_i^2 u_i

where :math:`v_i` is the :math:`i` -th one of :math:`V` 's column partition.

As the fomula above indicates the way to get eigenvalues(:math:`\sigma_i^2`)
and characteristic vectors( :math:`u_i` ) of :math:`A A^T` , we can calculate out
:math:`abs \Sigma` and :math:`U` .
"""
end
end

max_stars_repo_name: thautwarm/FlowerQ

max_issues_repo_name: thautwarm/FlowerQ

max_forks_repo_name: thautwarm/FlowerQ


text:
include("../src/twoscales_pure_ab.jl")
include("../src/henon_heiles.jl")
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i = 1:size(tab, 1)
        for j = (i+1):size(tab, 1)
            nm = norm(tab[i] - tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, nmmin
end

function fctmain(n_tau, prec)

    setprecision(prec)
    Random.seed!(5612)
    u0 = rand(BigFloat, 4)

    t_max = big"1.0"
    epsilon = big"1e-3"
    println("epsilon=$epsilon")
    nbmaxtest = 10
    ordmax = 12
    debord = 3
    pasord = 1
    y = ones(Float64, nbmaxtest, div(ordmax - debord, pasord) + 1)
    fct = u -> [u[2]^2, 1 / (1 + u[3]^2), 0.5 - u[4], 0.8u[1]]
    res_gen = Array{Array{BigFloat,1},2}(undef, nbmaxtest, div(ordmax - debord, pasord) + 1)
    x = zeros(Float64, nbmaxtest)
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    parphi = PreparePhi(epsilon, n_tau, A, fct)
    nm = NaN
    ordmax += 1
    ordprep = undef
    nb = 20 * 2^(nbmaxtest)
    solref = undef
    while isnan(nm)
        ordmax -= 1
        ordprep = ordmax + 2
        @time par_u0 = PrepareU0(parphi, ordprep, u0)
        @time pargen = PrepareTwoScalePureAB(nb, t_max, ordmax, par_u0)
        @time solref = twoscales_pure_ab(pargen, only_end = true, diff_fft = false)
        nm = norm(solref, Inf)
    end

    tabsol = Array{Array{BigFloat,1},1}(undef, 1)

    tabsol[1] = solref

    println("ordmax=$ordmax")

    println("solref=$solref")
    solref_gen = solref
    indref = 1

    ind = 1
    for order = debord:pasord:ordmax
        ordprep = order + 2
        println("eps=$epsilon solRef=$solref order=$order")
        nb = 10
        indc = 1
        labels = Array{String,2}(undef, 1, order - debord + 1)
        resnorm = 0
        resnormprec = 1
        sol = undef
        println("preparation ordre $order + 2")
        @time par_u0 = PrepareU0(parphi, ordprep, u0)
        while indc <= nbmaxtest
            @time pargen = PrepareTwoScalePureAB(nb, t_max, order, par_u0)
            @time sol = twoscales_pure_ab(pargen, only_end = true, diff_fft = false)
            push!(tabsol, sol)
            res_gen[indc, ind] = sol
            (a, b), nm = getmindif(tabsol)
            if a != indref
                println("New solref !!!! a=$a, b=$b nm=$nm")
                indref = a
                solref = tabsol[a]
                for i = 1:ind
                    borne = (i < ind) ? size(res_gen, 1) : indc
                    for j = 1:borne
                        nm2 = min(norm(res_gen[j, i] - solref, Inf), 1.1)
                        y[j, i] = nm2 == 0 ? nm : nm2
                    end
                end
            else
                diff = solref - sol
                y[indc, ind] = min(norm(diff, Inf), 1.1)
            end
            x[indc] = 1.0 / nb
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i = debord:pasord:order
            labels[1, (i-debord)÷pasord+1] = " eps,order=$(convert(Float32,epsilon)),$i "
        end
        p = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "delta t",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :topleft,
            label = labels,
            marker = 2,
        )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32, epsilon)
        Plots.savefig(p, "out/r2_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_fct3.pdf")
        if resnorm > resnormprec
            break
        end
        resnormprec = resnorm
        ind += 1
    end
end

# testODESolver()

# for i=3:9
#     fctMain(2^i)
# end
fctmain(256, 512)


text:
export place, cm_types, cm_type

################################################################################
#
#  CM theory
#
################################################################################

# A primitive embedding type

# Alternatively, it could just be an Int + Bool
# With the Int referring to complex_places(K)[i]
mutable struct NumFieldEmbedding
  field::AnticNumberField
  plc::InfPlc
  isconj::Bool
end

Base.:(==)(f::NumFieldEmbedding, g::NumFieldEmbedding) = (f.plc == g.plc) &&
                                                         (f.isconj == g.isconj)

number_field(f::NumFieldEmbedding) = f.field

place(f::NumFieldEmbedding) = f.plc

function _print_acb_neatly(io, x::acb)
  if !_isimag(x)
    a = Float64(real(x))
    print(io, @sprintf "%.2f" a)
  end
  if !isreal(x)
    b = Float64(imag(x))
    print(io, _isimag(x) ? "" : " + ", @sprintf("%.2f", b), " * i")
  end
end

function _isimag(x::acb)
  z = arb()
  ccall((:acb_get_imag, libarb), Cvoid, (Ref{arb}, Ref{acb}), z, x)
  return iszero(z)
end

function Base.show(io::IO, ::MIME"text/plain", f::NumFieldEmbedding)
  print(io, "Embedding of\n", number_field(f), "\n")
  print(io, "corresponding to ≈ ")
  _print_acb_neatly(io, f.isconj ? conj(place(f).r) : place(f).r)
end

function Base.show(io::IO, f::NumFieldEmbedding)
  print(io, "Embedding corresponding to ≈ ")
  _print_acb_neatly(io, f.isconj ? conj(place(f).r) : place(f).r)
end

function (f::NumFieldEmbedding)(x::nf_elem)
  @req number_field(f) === parent(x) "Element not in the domain of the embedding"
  xx = evaluate(x, place(f))
  return f.isconj ? conj(xx) : xx
end

conj(E::NumFieldEmbedding) = NumFieldEmbedding(E.field, E.plc, !E.isconj)

function complex_embeddings(K::AnticNumberField)
  cpls = get_attribute(K, :complex_embeddings)
  if cpls !== nothing
    return cpls::Vector{NumFieldEmbedding}
  end
  res = NumFieldEmbedding[]
  for p in infinite_places(K)
    push!(res, NumFieldEmbedding(K, p, false))
    push!(res, NumFieldEmbedding(K, p, true))
  end
  set_attribute!(K, :complex_embeddings => res)
  return res
end

################################################################################
#
#  Restriction
#
################################################################################

function restrict(E::NumFieldEmbedding, f::NfToNfMor)
  k = domain(f)
  @req iscm_field(k)[1] "Subfield must be a CM field"
  kem = complex_embeddings(k)
  a = gen(k)
  b = E(f(a))
  # This is not optimal, but guarded against precision issues
  cn = [overlaps(b, p(a)) for p in kem]
  @assert count(cn) == 1

  i = findfirst(cn)
  return kem[i]
end

################################################################################
#
#  CM types
#
################################################################################

# A good reference is
# [Streng2010]: Marco Streng, Complex multiplication of abelian surfaces, 2010

mutable struct CMType
  field::AnticNumberField
  embeddings::Vector{NumFieldEmbedding}

  function CMType(K::AnticNumberField, embeddings::Vector{NumFieldEmbedding})
    z = new(K, embeddings)
    return z
  end
end

function cm_type(K::AnticNumberField, embeddings::Vector{NumFieldEmbedding})
  @req iscm_field(K)[1] "Field must a CM field"
  @req 2 * length(embeddings) == degree(K) "Wrong number of embeddings"
  @req length(unique(map(x -> x.plc, embeddings))) == div(degree(K), 2) "Embeddings must be pairwise non-conjugated"
  return CMType(K, embeddings)
end

Base.:(==)(f::CMType, g::CMType) = (f.field === g.field) &&
                                      f.embeddings == g.embeddings

number_field(C::CMType) = C.field

embeddings(C::CMType) = C.embeddings

################################################################################
#
#  Induction
#
################################################################################

function induce(C::CMType, f::NfToNfMor)
  @assert C.field == domain(f)
  K = codomain(f)
  res = NumFieldEmbedding[]
  for E in complex_embeddings(K)
    e = restrict(E, f)
    if e in C.embeddings
      push!(res, E)
    end
  end
  return CMType(K, res)
end

function isinduced(C::CMType, f::NfToNfMor)
  k = domain(f)
  fl, _ = Hecke.iscm_field(k)
  for D in cm_types(k)
    if induce(D, f) == C
      return true, D
    end
  end
  return false, C
end

function isprimitive(C::CMType)
  for (k, ktoK) in subfields(C.field)
    if degree(k) == degree(C.field) || !iscm_field(k)[1]
      continue
    end
    if isinduced(C, ktoK)[1]
      return false
    end
  end
  return true
end

function cm_types(K::AnticNumberField)
  fl, _ = iscm_field(K)
  @assert fl
  g = div(degree(K), 2)
  cpl = complex_places(K)
  it = cartesian_product_iterator([true, false], g)
  res = CMType[]
  for x in it
    push!(res, CMType(K, [NumFieldEmbedding(K, cpl[i], x[i]) for i in 1:g]))
  end
  return res
end

function Base.:(*)(f::NfToNfMor, E::NumFieldEmbedding)
  k = domain(f)
  @assert k === codomain(f)
  @assert k === E.field
  a = gen(k)
  b = E(f(a))
  kplc = complex_places(k)
  for p in kplc
    if overlaps(b, evaluate(a, p))
      return NumFieldEmbedding(k, p, false)
    elseif overlaps(b, conj(evaluate(a, p)))
      return NumFieldEmbedding(k, p, true)
    end
  end
end

function Base.:(*)(f::NfToNfMor, C::CMType)
  return CMType(domain(f), [f * E for E in C.embeddings])
end

function reflex(c::CMType)
  K = c.field
  N, KtoN = normal_closure(K)
  cind = induce(c, KtoN)
  A = automorphisms(N)
  a = gen(N)
  cp = complex_places(N)
  P = cp[1] # lets one place of N to identify N with a subset of C
  res = NumFieldEmbedding[]

  for E in cind.embeddings
    b = E(a)
    _i = [ overlaps(evaluate(f(a), P), b) for f in A ]
    @assert count(_i) == 1
    i = findfirst(_i)
    f = A[i]
    finv = inv(f)
    c = evaluate(finv(a), P)
    _i = vcat([ overlaps(c, evaluate(a, P)) for P in cp ], [overlaps(c, conj(evaluate(a, P))) for P in cp])
    @assert count(_i) == 1
    for P in cp
      if overlaps(c, evaluate(a, P))
        push!(res, NumFieldEmbedding(N, P, false))
      elseif overlaps(c, conj(evaluate(a, P)))
        push!(res, NumFieldEmbedding(N, P, true))
      end
    end
  end
  cinv = CMType(N, res)

  for (k, ktoK) in subfields_normal(N)
    if degree(k) == degree(N) || !iscm_field(k)[1]
      continue
    end
    fl, D = isinduced(cinv, ktoK)
    if fl
      return D
    end
  end
end


text:
function assemble(
    problem::StiffnessTopOptProblem{dim,T},
    elementinfo::ElementFEAInfo{dim,T},
    vars=ones(T, getncells(getdh(problem).grid)),
    penalty=PowerPenalty(T(1)),
    xmin=T(0.001),
) where {dim,T}
    globalinfo = GlobalFEAInfo(problem)
    assemble!(globalinfo, problem, elementinfo, vars, penalty, xmin)
    return globalinfo
end

function assemble!(
    globalinfo::GlobalFEAInfo{T},
    problem::StiffnessTopOptProblem{dim,T},
    elementinfo::ElementFEAInfo{dim,T,TK},
    vars=ones(T, getncells(getdh(problem).grid)),
    penalty=PowerPenalty(T(1)),
    xmin=T(0.001);
    assemble_f=true,
) where {dim,T,TK}
    ch = problem.ch
    dh = ch.dh
    K, f = globalinfo.K, globalinfo.f
    if assemble_f
        f .= elementinfo.fixedload
    end
    Kes, fes = elementinfo.Kes, elementinfo.fes
    black = problem.black
    white = problem.white
    varind = problem.varind

    _K = K isa Symmetric ? K.data : K
    _K.nzval .= 0
    assembler = Ferrite.AssemblerSparsityPattern(_K, f, Int[], Int[])

    global_dofs = zeros(Int, ndofs_per_cell(dh))
    fe = zeros(typeof(fes[1]))
    Ke = zeros(T, size(rawmatrix(Kes[1])))

    celliterator = CellIterator(dh)
    for (i, cell) in enumerate(celliterator)
        # get global_dofs for cell#i
        celldofs!(global_dofs, dh, i)
        fe = fes[i]
        _Ke = rawmatrix(Kes[i])
        Ke = _Ke isa Symmetric ? _Ke.data : _Ke
        if black[i]
            if assemble_f
                Ferrite.assemble!(assembler, global_dofs, Ke, fe)
            else
                Ferrite.assemble!(assembler, global_dofs, Ke)
            end
        elseif white[i]
            if PENALTY_BEFORE_INTERPOLATION
                px = xmin
            else
                px = penalty(xmin)
            end
            Ke = px * Ke
            if assemble_f
                fe = px * fe
                Ferrite.assemble!(assembler, global_dofs, Ke, fe)
            else
                Ferrite.assemble!(assembler, global_dofs, Ke)
            end
        else
            if PENALTY_BEFORE_INTERPOLATION
                px = density(penalty(vars[varind[i]]), xmin)
            else
                px = penalty(density(vars[varind[i]], xmin))
            end
            Ke = px * Ke
            if assemble_f
                fe = px * fe
                Ferrite.assemble!(assembler, global_dofs, Ke, fe)
            else
                Ferrite.assemble!(assembler, global_dofs, Ke)
            end
        end
    end

    #* apply boundary condition
    _K = TK <: Symmetric ? K.data : K
    apply!(_K, f, ch)

    return nothing
end

function assemble_f(
    problem::StiffnessTopOptProblem{dim,T},
    elementinfo::ElementFEAInfo{dim,T},
    vars::AbstractVector{T},
    penalty,
    xmin=T(1) / 1000,
) where {dim,T}
    f = get_f(problem, vars)
    assemble_f!(f, problem, elementinfo, vars, penalty, xmin)
    return f
end
get_f(problem, vars::Array) = zeros(T, ndofs(problem.ch.dh))

function assemble_f!(
    f::AbstractVector,
    problem::StiffnessTopOptProblem,
    elementinfo::ElementFEAInfo,
    vars::AbstractVector,
    penalty,
    xmin,
)
    black = elementinfo.black
    white = elementinfo.white
    varind = elementinfo.varind
    fes = elementinfo.fes

    dof_cells = elementinfo.metadata.dof_cells

    update_f!(
        f, fes, elementinfo.fixedload, dof_cells, black, white, penalty, vars, varind, xmin
    )
    return f
end

function update_f!(
    f::Vector, fes, fixedload, dof_cells, black, white, penalty, vars, varind, xmin
)
    @inbounds for dofidx in 1:length(f)
        f[dofidx] = fixedload[dofidx]
        r = dof_cells.offsets[dofidx]:(dof_cells.offsets[dofidx + 1] - 1)
        for i in r
            cellidx, localidx = dof_cells.values[i]
            if black[cellidx]
                f[dofidx] += fes[cellidx][localidx]
            elseif white[cellidx]
                if PENALTY_BEFORE_INTERPOLATION
                    px = xmin
                else
                    px = penalty(xmin)
                end
                f[dofidx] += px * fes[cellidx][localidx]
            else
                if PENALTY_BEFORE_INTERPOLATION
                    px = density(penalty(vars[varind[cellidx]]), xmin)
                else
                    px = penalty(density(vars[varind[cellidx]], xmin))
                end
                f[dofidx] += px * fes[cellidx][localidx]
            end
        end
    end

    return nothing
end

function assemble_f!(f::AbstractVector, problem, dloads)
    metadata = problem.metadata
    dof_cells = metadata.dof_cells
    update_f!(f, dof_cells, dloads)
    return f
end

function update_f!(f::Vector, dof_cells, dloads)
    for dofidx in 1:length(f)
        r = dof_cells.offsets[dofidx]:(dof_cells.offsets[dofidx + 1] - 1)
        for i in r
            cellidx, localidx = dof_cells.values[i]
            f[dofidx] += dloads[cellidx][localidx]
        end
    end
    return nothing
end

#=
function update_f!(f::CuVector, dof_cells, dloads)
    args = (f, dof_cells.offsets, dof_cells.values, dloads)
    callkernel(dev, assemble_kernel2, args)
    CUDAdrv.synchronize(ctx)

    return
end

function assemble_kernel2(f, dof_cells_offsets, dof_cells_values, dloads)
    i = @thread_global_index()
    offset = @total_threads()
    @inbounds while i <= length(f)
        r = dof_cells_offsets[i] : dof_cells_offsets[i+1]-1
        for i in r
            cellidx, localidx = dof_cells_values[i]
            f[i] += dloads[cellidx][localidx]
        end
        i += offset
    end
    return
end
=#


text:
"""
```
IntegrationParameters(;
    method::AbstractIntegrator,
    tmin = 0.0,
    tmax = 1e13,
    dtmin = 1e-3,
    dtmax = Inf,
    abstol = 1e-6,
    reltol = 1e-6,
    maxchange = 1.2,
    exponent = 20.0,
    maxiter = 10,
)
```
Creates [`IntegrationParameters`](@ref).
"""
function IntegrationParameters(;
    method::AbstractIntegrator,
    tmin = 0.0,
    tmax = 1e13,
    dtmin = 1e-3,
    dtmax = Inf,
    abstol = 1e-6,
    reltol = 1e-6,
    maxchange = 1.2,
    exponent = 20.0,
    maxiter = 10,
)
    return IntegrationParameters(
        method,
        tmin,
        tmax,
        dtmin,
        dtmax,
        abstol,
        reltol,
        maxchange,
        exponent,
        maxiter,
    )
end

"""
```
IntegrationTime(; dt = 0.0, time = 0.0, step = 0)
```
Creates [`IntegrationTime`](@ref).
"""
function IntegrationTime(; dt = 0.0, time = 0.0, step = 0)
    return IntegrationTime(dt, time, step)
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
macro expect(ch...)
    N = length(ch)
    uint8s = map(UInt8, ch)
    return esc(quote
        Base.@nexprs $N i->begin
            eof(io) && throw(ArgumentError("eof encountered before $T was able to be fully parsed"))
            b = JSON2.readbyte(io)
            b == $(uint8s)[i] || throw(JSON2.invalid(T, b))
        end
    end)
end

macro expectoneof(ch...)
    uint8s = map(UInt8, ch)
    return esc(quote
        if !eof(io)
            b = JSON2.readbyte(io)
            b in $uint8s || throw(ArgumentError("invalid JSON, encountered '$(Char(b))' expected one of: '$(join($ch, "' '"))'"))
        end
    end)
end

@inline iswh(b) = b == UInt8('\t') || b == UInt8(' ') || b == UInt8('\n') || b == UInt8('\r')

@inline function wh!(io)
    b = peekbyte(io)
    while iswh(b)
        readbyte(io)
        b = peekbyte(io)
    end
    return
end

invalid(T, b) = ArgumentError("invalid JSON detected parsing type '$T': encountered '$(Char(b))'")

read(str::AbstractString, T=Any; kwargs...) = read(IOBuffer(str), T; kwargs...)

# read generic JSON: detect null, Bool, Int, Float, String, Array, Dict/Object into a NamedTuple
function read(io::IO, ::Type{Any}=Any; kwargs...)
    eof(io) && return NamedTuple()
    wh!(io)
    b = peekbyte(io)
    if b == UInt('{')
        # object
        return read(io, NamedTuple; kwargs...)
    elseif b == UInt('[')
        # array
        # TODO: we could try and be fancy here and parse a more specific vector
        return read(io, Vector{Any}; kwargs...)
    elseif (NEG_ONE < b < TEN) || (b == MINUS || b == PLUS)
        # int or float
        fl = read(io, Float64; kwargs...)
        int = unsafe_trunc(Int, fl)
        return ifelse(int == fl, int, fl)
    elseif b == UInt8('"')
        # string)
        return read(io, String; kwargs...)
    elseif b == UInt8('n')
        # null
        return read(io, Nothing; kwargs...)
    elseif b == UInt8('t')
        # true
        return read(io, Bool; kwargs...)
    elseif b == UInt8('f')
        # false or function literal
        pos = position(io)
        readbyte(io)
        eof(io) && throw(ArgumentError("early EOF"))
        func = peekbyte(io) == UInt8('u')
        seek(io, pos)
        return func ? read(io, Function; kwargs...) : read(io, Bool; kwargs...)
    else
        throw(ArgumentError("error detecting type of JSON object to parse: invalid JSON detected"))
    end
end

nonNothingT(::Type{Union{Nothing, T}}) where {T} = T
function read(io::IO, U::Union; kwargs...)
    if U.a === Nothing || U.b === Nothing
        b = peekbyte(io)
        if b == LITTLE_N
            return read(io, Nothing; kwargs...)
        else
            return read(io, nonNothingT(U); kwargs...)
        end
    else
        pos = position(io)
        try
            return read(io, U.a; kwargs...)
        catch
            seek(io, pos)
            return read(io, U.b; kwargs...)
        end
    end
end

function read(io::IO, T::Type{NamedTuple}; kwargs...)
    @expect '{'
    wh!(io)
    keys = Symbol[]
    vals = Any[]
    peekbyte(io) == CLOSE_CURLY_BRACE && (readbyte(io); @goto done)
    while true
        push!(keys, read(io, Symbol; kwargs...))
        wh!(io)
        @expect ':'
        wh!(io)
        push!(vals, read(io, Any; kwargs...)) # recursively reads value
        wh!(io)
        @expectoneof ',' '}'
        b == CLOSE_CURLY_BRACE && @goto done
        wh!(io)
    end
    @label done
    return NamedTuple{Tuple(keys)}(Tuple(vals))
end

function read(io::IO, ::Type{T}; kwargs...) where {T <: NamedTuple{names, types}} where {names, types}
    @expect '{'
    wh!(io)
    keys = Symbol[]
    vals = Any[]
    peekbyte(io) == CLOSE_CURLY_BRACE && (readbyte(io); @goto done)
    typemap = Dict(k=>v for (k, v) in zip(names, types.parameters))
    while true
        key = read(io, Symbol; kwargs...)
        push!(keys, key)
        wh!(io)
        @expect ':'
        wh!(io)
        push!(vals, read(io, typemap[key]; kwargs...)) # recursively reads value
        wh!(io)
        @expectoneof ',' '}'
        b == CLOSE_CURLY_BRACE && @goto done
        wh!(io)
    end
    @label done
    return NamedTuple{names}(NamedTuple{Tuple(keys)}(Tuple(vals)))
end

read(io::IO, ::Type{T}; kwargs...) where {T <: AbstractDict} = read(io, T(); kwargs...)
function read(io::IO, dict::Dict{K,V}; kwargs...) where {K, V}
    T = typeof(dict)
    @expect '{'
    wh!(io)
    peekbyte(io) == CLOSE_CURLY_BRACE && (readbyte(io); @goto done)
    while true
        key = read(io, K; kwargs...)
        wh!(io)
        @expect ':'
        wh!(io)
        dict[key] = read(io, V; kwargs...) # recursively reads one element
        wh!(io)
        @expectoneof ',' '}'
        b == CLOSE_CURLY_BRACE && @goto done
        wh!(io)
    end
    @label done
    return dict
end

read(io::IO, ::Type{T}; kwargs...) where {T <: Tuple} = Tuple(read(io, Array; kwargs...))
read(io::IO, ::Type{T}; kwargs...) where {T <: AbstractSet} = T(read(io, Array; kwargs...))
read(io::IO, ::Type{T}; kwargs...) where {T <: AbstractArray} = read(io, T([]); kwargs...)

function read(io::IO, A::AbstractArray{eT}; kwargs...) where {eT}
    T = typeof(A)
    @expect '['
    wh!(io)
    peekbyte(io) == CLOSE_SQUARE_BRACE && (readbyte(io); @goto done)
    while true
        push!(A, read(io, eT; kwargs...)) # recursively reads one element
        wh!(io)
        @expectoneof ',' ']'
        b == CLOSE_SQUARE_BRACE && @goto done
        wh!(io)
    end
    @label done
    return A
end

function read(io::IO, ::Type{Function}; kwargs...)
    buf = IOBuffer()
    wh!(io)
    while !eof(io)
        b = readbyte(io)
        Base.write(buf, b)
        b == UInt8('{') && break
    end
    bracketcount = 1
    while !eof(io) && bracketcount > 0
        b = readbyte(io)
        Base.write(buf, b)
        bracketcount += b == UInt8('{') ? 1 : b == UInt8('}') ? -1 : 0
    end
    return Function(String(take!(buf)))
end

# read Number, String, Bool
read(io::IO, ::Type{T}; kwargs...) where {T <: Integer} = Parsers.parse(T, io)
read(io::IO, ::Type{T}; kwargs...) where {T <: AbstractFloat} = Parsers.parse(T, io)

function read(io::IO, T::Type{Char}; kwargs...)
    @expect '"'
    c = Char(readbyte(io))
    @expect '"'
    return c
end
read(io::IO, ::Type{Date}; dateformat=Dates.ISODateFormat, kwargs...) = Date(read(io, String; kwargs...), dateformat)
function read(io::IO, ::Type{DateTime}; read_datetimeformats=[], dateformat=Dates.ISODateTimeFormat, kwargs...)
    dfs = [read_datetimeformats; dateformat]
    s = read(io, String; kwargs...)
    for df in dfs
        try
            return DateTime(s, df)
        catch ex
        end
    end
    error("Unable to parse datetime: $s")
end
read(io::IO, ::Type{T}; kwargs...) where {T <: Enum} = Core.eval(parentmodule(T), read(io, Symbol; kwargs...))

function read(io::IO, T::Type{Nothing}; kwargs...)
    @expect 'n' 'u' 'l' 'l'
    return nothing
end

function read(io::IO, T::Type{Missing}; kwargs...)
    @expect 'n' 'u' 'l' 'l'
    return missing
end

function read(io::IO, T::Type{Bool}; kwargs...)
    b = peekbyte(io)
    if b == UInt8('t')
        @expect 't' 'r' 'u' 'e'
        return true
    else
        @expect 'f' 'a' 'l' 's' 'e'
        return false
    end
end

function generate_default_read_body(N, types, jsontypes, isnamedtuple)
    inner = Expr(:block, :(kwargs = JSON2.mergedefaultkwargs(T; kwargs...)))
    keys = ((((Symbol("key_$j") for j = 1:i)...,) for i = 1:N)...,)
    vals = ((((Symbol("val_$j") for j = 1:i)...,) for i = 1:N)...,)
    foreach(1:N) do i
        if isnamedtuple
            ret = :(T(($((vals[i])...),)))
        else
            ret = :(T($((vals[i])...)))
        end
        push!(inner.args, quote
            $(keys[i][i]) = JSON2.read(io, String; kwargs...)
            JSON2.wh!(io)
            JSON2.@expect ':'
            JSON2.wh!(io)
            $(vals[i][i]) = $(JSON2.getconvert(types[i]))(JSON2.read(io, $(jsontypes[i]); kwargs...))
            JSON2.wh!(io)
            JSON2.@expectoneof ',' '}'
            b == JSON2.CLOSE_CURLY_BRACE && return $ret
            JSON2.wh!(io)
        end)
    end
    if isnamedtuple
        ret = :(T(($((vals[N])...),)))
    else
        ret = :(T($((vals[N])...)))
    end
    body = quote
        JSON2.@expect '{'
        JSON2.wh!(io)
        JSON2.peekbyte(io) == JSON2.CLOSE_CURLY_BRACE && (JSON2.readbyte(io); return T())
        $inner
        # in case there are extra fields, just ignore
        curlies = 1
        while !eof(io)
            b = JSON2.readbyte(io)
            if b == JSON2.OPEN_CURLY_BRACE
                curlies += 1
            elseif b == JSON2.CLOSE_CURLY_BRACE
                curlies -= 1
            end
            curlies == 0 && return $ret
        end
        throw(ArgumentError("failed to parse $T from JSON"))
    end
    # @show body
    return body
end

const EMPTY_SYMBOL = Symbol("")

function read_args(io, b, name, names, T, types, jT, jsontypes, defaults, fulltypes, fulljsontypes, args=(), argnm=EMPTY_SYMBOL, arg=nothing; kwargs...)
    # @show name, names, T, types, jT, jsontypes, defaults, args, argnm, arg
    # @show name == argnm
    # @show names
    # @show isempty(names)
    if argnm != EMPTY_SYMBOL
        if name == argnm
            if isempty(names)
                return ((args..., arg), b)
            else
                return JSON2.read_args(io, b,
                    names[1], Base.tail(names),
                    types[1], Base.tail(types),
                    jsontypes[1], Base.tail(jsontypes),
                    defaults, fulltypes, fulljsontypes,
                    (args..., arg), EMPTY_SYMBOL, nothing;
                    kwargs...)
            end
        else
            if isempty(names)
                return ((args..., defaults[name]), b)
            else
                return JSON2.read_args(io, b,
                    names[1], Base.tail(names),
                    types[1], Base.tail(types),
                    jsontypes[1], Base.tail(jsontypes),
                    defaults, fulltypes, fulljsontypes,
                    (args..., defaults[name]), argnm, arg;
                    kwargs...)
            end
        end
    end
    JSON2.wh!(io)
    if b == JSON2.CLOSE_CURLY_BRACE
        if isempty(names)
            return ((args..., defaults[name]), b)
        else
            return JSON2.read_args(io, b, names[1], Base.tail(names),
                                   types[1], Base.tail(types),
                                   jsontypes[1], Base.tail(jsontypes),
                                   defaults, fulltypes, fulljsontypes,
                                   (args..., defaults[name]), argnm, arg;
                                   kwargs...)
        end
    end
    key = JSON2.read(io, Symbol; kwargs...)
    JSON2.wh!(io)
    JSON2.@expect ':'
    JSON2.wh!(io)
    if key == name
        val = JSON2.getconvert(T)(JSON2.read(io, jT; kwargs...))
    else
        val = defaults[name]
        if haskey(fulltypes, key)
            argnm = key
            arg = JSON2.getconvert(fulltypes[key])(JSON2.read(io, fulljsontypes[key]; kwargs...))
        else
            # ignore unknown field
            JSON2.read(io, Any; kwargs...)
        end
    end
    JSON2.wh!(io)
    JSON2.@expectoneof ',' '}'
    if isempty(names)
        return ((args..., val), b)
    else
        return JSON2.read_args(io, b,
            names[1], Base.tail(names),
            types[1], Base.tail(types),
            jsontypes[1], Base.tail(jsontypes),
            defaults, fulltypes, fulljsontypes,
            (args..., val), argnm, arg;
            kwargs...)
    end
end

function generate_missing_read_body(names, types, jsontypes, defaults)
    fulltypes = NamedTuple{names}(types)
    fulljsontypes = NamedTuple{names}(jsontypes)
    body = quote
        kwargs = JSON2.mergedefaultkwargs(T; kwargs...)
        JSON2.@expect '{'
        JSON2.wh!(io)
        JSON2.peekbyte(io) == JSON2.CLOSE_CURLY_BRACE && (JSON2.readbyte(io); return T())
        args, b = JSON2.read_args(io, b,
            $(QuoteNode(names[1])), $(Base.tail(names)),
            $(types[1]), $(Base.tail(types)),
            $(jsontypes[1]), $(Base.tail(jsontypes)),
            $defaults, $fulltypes, $fulljsontypes; kwargs...)
        b == JSON2.CLOSE_CURLY_BRACE && return T(args...)
        # in case there are extra fields, just ignore
        curlies = 1
        while !eof(io)
            b = JSON2.readbyte(io)
            if b == JSON2.OPEN_CURLY_BRACE
                curlies += 1
            elseif b == JSON2.CLOSE_CURLY_BRACE
                curlies -= 1
            end
            curlies == 0 && return T(args...)
        end
        throw(ArgumentError("failed to parse $T from JSON"))
    end
    # @show body
    return body
end

function generate_read_body_noargs(N, names, jsontypes, defaults)
    body = quote
        kwargs = JSON2.mergedefaultkwargs(T; kwargs...)
        x = T()
        eof(io) && return x
        JSON2.@expect '{' # start of object
        JSON2.wh!(io)
        JSON2.peekbyte(io) == JSON2.CLOSE_CURLY_BRACE && (JSON2.readbyte(io); return x)
        names = $names
        defaults = $defaults
        allnames = Set($(names))
        seen = Set{Symbol}()
        jsontypes = $jsontypes
        while !eof(io)
            key = JSON2.read(io, Symbol; kwargs...)
            key = get(names, key, key)
            JSON2.wh!(io)
            JSON2.@expect ':'
            JSON2.wh!(io)
            val = JSON2.read(io, get(jsontypes, key, Any); kwargs...)
            key in allnames && Core.setfield!(x, key, val)
            push!(seen, key)
            JSON2.wh!(io)
            JSON2.@expectoneof ',' '}'
            b == JSON2.CLOSE_CURLY_BRACE && break
            JSON2.wh!(io)
        end
        for i = 1:$N
            key = fieldname(T, i)
            if !(key in seen) && haskey(defaults, key)
                Core.setfield!(x, key, defaults[key])
            end
        end
        b == JSON2.CLOSE_CURLY_BRACE && return x
        throw(ArgumentError("failed to parse $T from JSON"))
    end
    # @show body
    return body
end

@generated function read(io::IO, ::Type{T}; kwargs...) where T
    N = fieldcount(T)
    types = Tuple(fieldtype(T, i) for i = 1:N)
    return generate_default_read_body(N, types, types, T <: NamedTuple)
end


text:
using Printf
import Flux: σ
using ModelingToolkit
using GalacticOptim
using Optim
using DiffEqFlux
using NeuralPDE
using Quadrature, Cubature, Cuba
using Plots
using IfElse

@parameters t,x,y
@variables c(..)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dx = Differential(x)
Dy = Differential(y)

# Parameters
u = 1.0
v = 1.0
v_vector = (u, v)
R = 0
D = 0 # 0.1 # diffusion
t_min = 0.0
t_max = 2.0
x_min = -1.0
x_max = 1.0
y_min = -1.0
y_max = 1.0
d = y_max-y_min

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

# exp(-(x^2+y^2)/0.1)
# div(v_vector * c) = dx(uc) + dy(vc)
# Equations, initial and boundary conditions
eqs = Dt(c(t,x,y)) ~ D * (Dxx(c(t,x,y)) + Dyy(c(t,x,y))) - (u*Dx(c(t,x,y)) + v*Dy(c(t,x,y))) + R

bcs = [
        c(0, x, y) ~ exp(-(x^2+y^2)/0.1), #cos(π*x) * cos(π*y) + 1.0 #,
        c(t, x, y_min) ~ c(t, x, y_max),
        c(t, x_min, y) ~ c(t, x_max, y)
      ]
#first step
domains = [t ∈ IntervalDomain(0.0,0.01),
           x ∈ IntervalDomain(x_min,x_max),
           y ∈ IntervalDomain(y_min,y_max)
          ]

dim = length(domains)
hidden = 16
chain = FastChain(FastDense(dim, hidden, Flux.σ),
                  FastDense(hidden, hidden, Flux.σ),
                  FastDense(hidden, 1))

quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=CubaCuhre(),
                                                   reltol = 1e-4, abstol = 1e-2,
                                                   maxiters = 20, batch=100)

discretization = NeuralPDE.PhysicsInformedNN(chain,quadrature_strategy)
pde_system = PDESystem(eqs, bcs, domains, [t,x,y], [c])
prob = NeuralPDE.discretize(pde_system,discretization)
sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)

# res = GalacticOptim.solve(prob,ADAM(0.1);cb=cb,maxiters=2500)
# prob = remake(prob,u0=res.minimizer)
res = GalacticOptim.solve(prob,Optim.BFGS();cb=cb,maxiters=2000)

nx = 32
ny = 32
dx = (x_max-x_min) / (nx - 1)
dy = (y_max-y_min) / (ny -1)
dt = 0.01
phi = discretization.phi
ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

analytic_solve(t,x,y) = exp(-((x-t)^2+(y-t)^2)/0.1)

function plots_init()
    c_predict = reshape([ phi([0.0, x, y], res.minimizer)[1] for x in xs for y in ys], length(xs), length(ys))
    extract_f(x, y) = exp(-((x)^2+(y)^2)/0.1)
    c_real = reshape([extract_f(x,y) for x in xs for y in ys], (length(xs),length(ys)))
    p1 = plot(xs, ys, c_predict)
    p2 = plot(xs, ys, c_real)
    plot(p1,p2)
end

plots_init()
savefig("advection_diffusion_2d_init_comparison")

function plot_comparison()
    # Animate
    anim = @animate for (i, t) in enumerate(0:dt:t_max)
        @info "Animating frame $i..."
        title = @sprintf("Advection-diffusion t = %.3f", t)
        c_real = reshape([analytic_solve(t,x,y) for x in xs for y in ys], (length(xs),length(ys)))
        c_predict = reshape([phi([t, x, y], res.minimizer)[1] for x in xs for y in ys], length(xs), length(ys))
        error_c = abs.(c_predict .- c_real)
        title = @sprintf("predict")
        p1 = heatmap(xs, ys, c_predict, label="", title=title , xlims=(-1, 1), ylims=(-1, 1), color=:thermal, clims=(0, 1.))
        title = @sprintf("real")
        p2 = heatmap(xs, ys, c_real, label="", title=title , xlims=(-1, 1), ylims=(-1, 1), color=:thermal, clims=(0, 1.))
        title = @sprintf("error")
        p3 = heatmap(xs, ys, error_c, label="", title=title , xlims=(-1, 1), ylims=(-1, 1), color=:thermal)
        plot(p1,p2,p3)
    end
    gif(anim, "advection_diffusion_2d_pinn_comparison.gif", fps=15)
end

plot_comparison()

quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=CubaCuhre(),
                                                   reltol = 1e-8, abstol = 1e-6,
                                                   maxiters = 30)
for t_ in collect(0.1:0.2:t_max)
    domains = [t ∈ IntervalDomain(0.0,t_),
               x ∈ IntervalDomain(x_min,x_max),
               y ∈ IntervalDomain(y_min,y_max)
              ]

    pde_system = remake(pde_system; domain = domains )
    discretization = remake(discretization; init_params = res.minimizer)
    prob = NeuralPDE.discretize(pde_system,discretization)
    res = GalacticOptim.solve(prob,ADAM(0.1);cb=cb,maxiters=3000)
    # prob = remake(prob,u0=res.minimizer)
    # res = GalacticOptim.solve(prob,Optim.BFGS();cb=cb,maxiters=10)
    plot_comparison()
    println("done")
end

nx = 128
ny = 128
dx = (x_max-x_min) / (nx - 1)
dy = (y_max-y_min) / (ny -1)
xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains[2:3]]

plot_comparison()

# Animate
anim = @animate for (i, t) in enumerate(0:dt:t_max)
    @info "Animating frame $i..."
    c_real = reshape([analytic_solve(t,x,y) for x in xs for y in ys], (length(xs),length(ys)))
    title = @sprintf("Advection-diffusion extract t = %.3f", t)
    heatmap(xs, ys, c_real, label="", title=title , xlims=(-1, 1), ylims=(-1, 1), color=:thermal, clims=(0, 1.))
end
gif(anim, "advection_diffusion_2d_extract.gif", fps=15)

# Animate
anim = @animate for (i, t) in enumerate(0:dt:t_max)
    @info "Animating frame $i..."
    c_predict = reshape([phi([t, x, y], res.minimizer)[1] for x in xs for y in ys], length(xs), length(ys))
    title = @sprintf("Advection-diffusion predict t = %.3f", t)
    heatmap(xs, ys, c_predict, label="", title=title , xlims=(-1, 1), ylims=(-1, 1), color=:thermal, clims=(0, 1.0))
end
gif(anim, "advection_diffusion_pinn.gif", fps=15)

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
export AbstractSymmetryEmbedding

export SymmetryEmbedding

export embed
export element, elements
export symmetry
export generator_elements, generator_indices
export fractional_momentum

abstract type AbstractSymmetryEmbedding end

const SymmetryOrEmbedding = Union{<:AbstractSymmetry, <:AbstractSymmetryEmbedding}

"""
    SymmetryEmbedding{SymmetryType<:AbstractSymmetry}

# Fields
- `lattice::Lattice`
- `symmetry::SymmetryType`
- `elements::Vector{SitePermutation}`
"""
struct SymmetryEmbedding{SymmetryType}<:AbstractSymmetryEmbedding
    lattice::Lattice
    symmetry::SymmetryType
    elements::Vector{SitePermutation}

    @doc """
        SymmetryEmbedding(lattice::Lattice, symmetry::SymmetryType)

    Construct a `SymmetryEmbedding`.
    """
    function SymmetryEmbedding(
        lattice::Lattice,
        symmetry::SymmetryType,
    ) where {SymmetryType<:AbstractSymmetry}
        if !iscompatible(lattice, symmetry)
            throw(ArgumentError("lattice and symmetry are not compatible"))
        end
        elems = [embed(lattice, elem) for elem in elements(symmetry)]
        if !allunique(elems)
            if isa(symmetry, PointSymmetry)
                reduced_elements = SitePermutation[]
                for elem in elems
                    if !in(elem, reduced_elements)
                        push!(reduced_elements, elem)
                    end
                end
                reduced_group = FiniteGroup(generate_multiplication_table(reduced_elements))
                isomorphic_group_names = String[]
                for i in 1:32
                    target_symmetry = PointSymmetryDatabase.get(i)
                    if !isnothing(group_isomorphism(reduced_group, target_symmetry.group))
                        push!(isomorphic_group_names, target_symmetry.hermann_mauguin)
                    end
                end
                @warn "Lattice $(lattice.hypercube.shape_matrix) is too small for $(symmetry_name(symmetry)) (embedding not faithful).\n"*
                      "Reduced point group isomorphic to: $(join(isomorphic_group_names, ", "))" # COV_EXCL_LINE
            else
                @warn "Lattice $(lattice.hypercube.shape_matrix) is too small for $(symmetry_name(symmetry)) (embedding not faithful)." # COV_EXCL_LINE
            end
        end
        new{SymmetryType}(lattice, symmetry, elems)
    end
end

Base.eltype(::SymmetryEmbedding) = SitePermutation
Base.valtype(::SymmetryEmbedding) = SitePermutation
Base.eltype(::Type{<:SymmetryEmbedding}) = SitePermutation
Base.valtype(::Type{<:SymmetryEmbedding}) = SitePermutation

"""
    elements(symbed::SymmetryEmbedding)

Get elements of the symmetry embedding
"""
elements(symbed::SymmetryEmbedding) = symbed.elements

"""
    element(symbed::SymmetryEmbedding, g...)

Get the `g`th element of the symmetry embedding.
"""
element(symbed::SymmetryEmbedding, g...) = symbed.elements[g...]

"""
    symmetry(symbed::SymmetryEmbedding)

Return the underlying symmetry of the symmetry embedding
"""
symmetry(symbed::SymmetryEmbedding) = symbed.symmetry

"""
    generator_elements(symbed::SymmetryEmbedding)

Get the generators of the symmetry embedding, based on the generators of the underlying
symmetry. The embedding is allowed to be one-to-one, in which case the return value may
contain duplicates.
"""
generator_elements(symbed::SymmetryEmbedding) = element(symbed, generator_indices(symbed))

function Base.:(==)(lhs::SymmetryEmbedding{S}, rhs::SymmetryEmbedding{S}) where S
    return lhs.lattice == rhs.lattice && lhs.symmetry == rhs.symmetry && lhs.elements == rhs.elements
end

Base.iterate(sym::SymmetryEmbedding) = iterate(elements(sym))
Base.iterate(sym::SymmetryEmbedding, i) = iterate(elements(sym), i)

Base.length(sym::SymmetryEmbedding) = length(elements(sym))

for f in [
    :group_order,
    :group_multiplication_table,
    :element_names,
    :element_name,
    :generator_indices,
    :character_table,
    :irreps,
    :irrep,
    :irrep_dimension,
    :num_irreps, :numirreps, :irrepcount,
]
    eval(quote
        @doc """
            $($f)(::SymmetryEmbedding, args...)

        Calls `$($f)` for the underlying symmetry.
        See [`$($f)(::FiniteTranslationSymmetry)`](@ref), [`$($f)(::PointSymmetry)`](@ref).
        """
        ($f)(symbed::SymmetryEmbedding, args...) = ($f)(symmetry(symbed), args...)
    end)
end

"""
    fractional_momentum(symbed::SymmetryEmbedding{FiniteTranslationSymmetry}, args...)

Return fractional momentum(s) of the translation symmetry embedding.
"""
function fractional_momentum(symbed::SymmetryEmbedding{FiniteTranslationSymmetry}, args...)
    return fractional_momentum(symmetry(symbed), args...)
end

"""
    iscompatible(tsymbed, psymbed)

Check whether the point symmetry embedding `psymbed` is compatible with
the translation symmetry embedding `tsymbed`, i.e. whether they have
the same "lattice".
"""
function iscompatible(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    psymbed::SymmetryEmbedding{PointSymmetry},
)::Bool
    return tsymbed.lattice == psymbed.lattice
end

"""
    iscompatible(tsymbed, tsym_irrep_index, psymbed)

Check whether the point symmetry embedding `psymbed` is compatible with
the translation symmetry irrep defined by `tsym_irrep_index` and `symmetry(tsymbed)`.
In other words, the little group elements
"""
function iscompatible(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    tsym_irrep_index::Integer,
    psymbed::SymmetryEmbedding{PointSymmetry},
)::Bool
    !iscompatible(tsymbed, psymbed) && return false
    return little_group_elements(tsymbed, tsym_irrep_index, psymbed) == 1:group_order(psymbed)
end

"""
    little_group_elements(tsymbed, psymbed)
"""
function little_group_elements(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    psymbed::SymmetryEmbedding{PointSymmetry},
)
    if !iscompatible(tsymbed, psymbed)
        throw(ArgumentError("translation and point symmetry-embeddings not compatible"))
    end
    return little_group_elements(symmetry(tsymbed), symmetry(psymbed))
end

"""
    little_group_elements(tsymbed, tsym_irrep_index, psymbed)
"""
function little_group_elements(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    tsym_irrep_index::Integer,
    psymbed::SymmetryEmbedding{PointSymmetry},
)
    if !iscompatible(tsymbed, psymbed)
        throw(ArgumentError("translation and point symmetry-embeddings not compatible"))
    end
    return little_group_elements(symmetry(tsymbed), tsym_irrep_index, symmetry(psymbed))
end

"""
    little_symmetry(tsymbed, psymbed)
"""
function little_symmetry(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    psymbed::SymmetryEmbedding{PointSymmetry},
)
    if !iscompatible(tsymbed, psymbed)
        throw(ArgumentError("translation and point symmetry-embeddings not compatible"))
    end
    psym_little = little_symmetry(symmetry(tsymbed), symmetry(psymbed))
    return SymmetryEmbedding(psymbed.lattice, psym_little)
end

# export little_symmetry_strong

# WIP
# """
#     little_symmetry_strong(tsymbed, psymbed)
# """
# function little_symmetry_strong(tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
#                                 psymbed::SymmetryEmbedding{PointSymmetry})
#     if !iscompatible(tsymbed, psymbed)
#         throw(ArgumentError("translation and point symmetry-embeddings not compatible"))
#     end
#     psym_little = little_symmetry(symmetry(tsymbed), symmetry(psymbed))

#     little_element_indices = Int[]
#     little_element_embedding = Set{SitePermutation}()

#     for (i_elem, elem) in enumerate(elements(psym_little))
#         embed_elem = embed(tsymbed.lattice, elem)
#         if i_elem == 1
#             @assert iscompatible(symmetry(tsymbed), elem)
#             push!(little_element_indices, i_elem)
#             push!(little_element_embedding, embed_elem)
#         elseif (iscompatible(symmetry(tsymbed), elem) &&
#                 !isidentity(embed_elem) &&
#                 !in(embed_elem, little_element_embedding))
#             push!(little_element_indices, i_elem)
#             change = true
#             while change
#                 change = false
#                 push!(little_element_embedding, embed_elem)
#                 for elem1 in elements(tsymbed), elem2 in little_element_embedding
#                     elem3 = elem1 * elem2
#                     if !in(elem3, little_element_embedding)
#                         change = true
#                         push!(reduced_group_elements, elem3)
#                         break
#                     end
#                 end
#             end
#         end
#     end
#     @show little_element_indices
#     @show element_name(symmetry(psymbed), little_element_indices)
#     psym_little2 = little_symmetry(symmetry(tsymbed), symmetry(psymbed), little_element_indices)
#     return SymmetryEmbedding(psymbed.lattice, psym_little2)
# end

"""
    little_symmetry(tsymbed, tsym_irrep_index, psymbed)
"""
function little_symmetry(
    tsymbed::SymmetryEmbedding{FiniteTranslationSymmetry},
    tsym_irrep_index::Integer,
    psymbed::SymmetryEmbedding{PointSymmetry},
)
    if !iscompatible(tsymbed, psymbed)
        throw(ArgumentError("translation and point symmetry-embeddings not compatible"))
    end
    psym_little = little_symmetry(symmetry(tsymbed), tsym_irrep_index, symmetry(psymbed))
    return SymmetryEmbedding(psymbed.lattice, psym_little)
end

"""
    symmetry_name(arg::SymmetryEmbedding)

Return the name of the symmetry embedding.
"""
function symmetry_name(arg::SymmetryEmbedding)
    return "Embed[$(symmetry_name(symmetry(arg))) on $(arg.lattice.hypercube.shape_matrix)"*
        " with $(numsite(arg.lattice.unitcell)) sites]"
end

"""
    embed(lattice::Lattice, tsym::FiniteTranslationSymmetry)

Embed translation symmetry `tsym` to lattice.
"""
embed(lattice::Lattice, tsym::FiniteTranslationSymmetry) = SymmetryEmbedding(lattice, tsym)

"""
    embed(lattice::Lattice, psym::PointSymmetry)

Embed point symmetry `psym` to lattice.
"""
embed(lattice::Lattice, psym::PointSymmetry) = SymmetryEmbedding(lattice, psym)

"""
    translation_symmetry_embedding(lattice::Lattice)
"""
function translation_symmetry_embedding(lattice::Lattice)
    tsym = FiniteTranslationSymmetry(lattice)
    return embed(lattice, tsym)
end

# struct SymmorphicSpaceSymmetryEmbedding<:AbstractSymmetryEmbedding
#     lattice::Lattice
#     translation_symmetry::FiniteTranslationSymmetry
#     point_symmetry::PointSymmetry
#     elements::Matrix{SitePermutation}

#     function SymmorphicSpaceSymmetryEmbedding(
#                 lattice::Lattice,
#                 translation::FiniteTranslationSymmetry,
#                 point::PointSymmetry)
#         if !iscompatible(translation, point)
#             throw(ArgumentError("translation symmetry and point symmetry not compatible"))
#         end
#         tels = [embed(lattice, elem) for elem in translation.elements]
#         pels = [embed(lattice, elem) for elem in point.elements]

#         if !allunique(pels)
#             throw(ArgumentError("lattice too small for the point symmetry operation (not faithful)"))
#         end

#         elements = [pop * top for top in tels, pop in pels]
#         @assert size(elements) == (length(tels), length(pels))
#         new(lattice, translation, point, elements)
#     end
# end

# function embed(lattice::Lattice, tsym::FiniteTranslationSymmetry, psym::PointSymmetry)
#     SymmorphicSpaceSymmetryEmbedding(lattice, tsym, psym)
# end

# function symmetry_name(arg::SymmorphicSpaceSymmetryEmbedding)
#     name1 = symmetry_name(symmetry(arg.component1))
#     name2 = symmetry_name(symmetry(arg.component2))
#     return "Embed[$name1 ⋊ $name2 on $(arg.lattice.hypercube.scale_matrix)]"
# end


text:
"""
Construct model for HydroGen with RunOfRiver Dispatch Formulation
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, D},
    ::Type{S},
) where {H <: HydroCascade, D <: HydroDispatchRunOfRiverCascade, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    #Variables
    PSI.add_variables!(optimization_container, PSI.ActivePowerVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.ReactivePowerVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.SpillageVariable, devices, D())

    #Constraints
    #=
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )=#
    flow_balance_cascade_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    #Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

"""
Construct model for HydroGen with RunOfRiver Dispatch Formulation
with only Active Power.
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, D},
    ::Type{S},
) where {
    H <: HydroCascade,
    D <: HydroDispatchRunOfRiverCascade,
    S <: PM.AbstractActivePowerModel,
}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    #Variables
    PSI.add_variables!(optimization_container, PSI.ActivePowerVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.SpillageVariable, devices, D())

    #Constraints
    #=
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )=#
    flow_balance_cascade_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    #Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end


text:
export global_positions!, global_position

"""
    global_positions!(g::BVHGraph)

Calculate the global positions for every vertex and frame in `g` and store them.

See also: [`global_position`](@ref)
"""
function global_positions!(g::BVHGraph)
    for v in vertices(g)
        positions!(g, v, zeros(Float64, frames(g), 3))
    end

    function calculate_positions(v::Integer, f::Integer, P::Matrix{Float64})
        o = offset(g, v)

        if outneighbors(g, v) != []
            R = rotation(g, v, f)
            M = P * [   R[1, 1] R[1, 2] R[1, 3] o[1];
                        R[2, 1] R[2, 2] R[2, 3] o[2];
                        R[3, 1] R[3, 2] R[3, 3] o[3];
                        0.0 0.0 0.0 1.0]

            positions(g, v)[f, :] = M[1:3, 4]

            for n in outneighbors(g, v)
                calculate_positions(n, f, M)
            end
        else
            vec = P * [o[1:3]..., 1.0]
            positions(g, v)[f, :] = vec[1:3]
        end
    end

    for f in 1:frames(g)
        o = offset(g)
        p = positions(g)[f, :]
        R = rotation(g, 1, f)
        M = [   R[1, 1] R[1, 2] R[1, 3] p[1] + o[1];
                R[2, 1] R[2, 2] R[2, 3] p[2] + o[2];
                R[3, 1] R[3, 2] R[3, 3] p[3] + o[3];
                0.0 0.0 0.0 1.0]

        positions(g, 1)[f, :] = M[1:3, 4]

        for n in outneighbors(g, 1)
            calculate_positions(n, f, M)
        end
    end

    return g
end

global_positions!() = g -> global_positions!(g)

"""
    global_position(g::BVHGraph, v::Integer, f::Integer, N::Matrix{Float64} = Matrix(1.0I, 4, 4))

Return the global position of vertex `v` for frame `f`.

This function does not store the result in `g`. `N` should not be changed.

See also: [`global_positions!`](@ref)
"""
function global_position(g::BVHGraph, v::Integer, f::Integer, N::Matrix{Float64} = Matrix(1.0I, 4, 4))
    if outneighbors(g, v) != []
        R = rotation(g, v, f)
    else
        R = one(RotMatrix{3, Float64})
    end

    if v != 1
        v₋₁ = inneighbors(g, v)[1]
        o = offset(g, v₋₁, v)
        A = [   R[1, 1] R[1, 2] R[1, 3] o[1];
                R[2, 1] R[2, 2] R[2, 3] o[2];
                R[3, 1] R[3, 2] R[3, 3] o[3];
                0.0 0.0 0.0 1.0] * N
        return global_position(g, v₋₁, f, A)
    else
        o = offset(g)
        pos = positions(g)[f, :]
        A = [   R[1, 1] R[1, 2] R[1, 3] pos[1] + o[1];
                R[2, 1] R[2, 2] R[2, 3] pos[2] + o[2];
                R[3, 1] R[3, 2] R[3, 3] pos[3] + o[3];
                0.0 0.0 0.0 1.0] * N
        return A[1:3, 4]
    end
end


text:
# Use $() syntax for string interpolation
# Put $() in a literal string
# If you have "zxywvut" and want to interpolate "Help!" in place of w, you could use "zxy$("Help!")wvut"
# The argument inside $() must be a string value or something that can be formatted as a string

# This works very similar to Python f-strings to mix variables and strings
my_name = "Viki"
ex1 = "My name is $(my_name)"
println(ex1)

# Multiple assignment with string interpolation
name1, name2 = "Joey", "Elliot"
ex2 = "$(name1) and $(name2) are friends"
println(ex2)

# To use a literal $ you need \
ex3 = "I have \$3"
println(ex3)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING


text:
# This file is part of FlagSets.jl project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

module FlagSets

#import Core.Intrinsics.bitcast
export FlagSet, flags, @flagset, @symbol_flagset

"""
    FlagSet{FlagType,BaseType<:Integer} <: AbstractSet{FlagType}

Supertype of sets which are subsets of a finite universe of objects (flags),
where each set is encoded as an integer of type `BaseType`.
Flag sets are implemented similarly to `BitSet`s.

If `T <: FlagSet{B,F}`, then `typemax(T)` is a finite set of flags of type `F`.
If `set isa T`, then `set` is a subset of `typemax(T)`.

See [`@flagset`](@ref) and [`@flagset_type`](@ref) for more information.
"""
abstract type FlagSet{FlagType,BaseType<:Integer} <: AbstractSet{FlagType} end

include("interface.jl")
include("macros.jl")

end # module


text:
function templaterun(
    fol :: AbstractString,
    src :: AbstractString = "";
    dst :: AbstractString = "."
)

    if src == ""; src = joinpath(DEPOT_PATH[1],"files","SAMTools","run_template") end
    if !isdir(src)
        tem = joinpath(@__DIR__,"..","templates","run_template")
        cp(tem,src)
        mkdir(joinpath(src,"LOGS"))
    end
    cp(src,joinpath(dst,fol),force=true)

    return

end

function templateexp(
    dst :: AbstractString = "."
)

    src = joinpath(DEPOT_PATH[1],"files","SAMTools","exp")
    if !isdir(src)
        tem = joinpath(@__DIR__,"..","templates","exp")
        cp(tem,src)
    end
    cp(src,joinpath(dst,"exp"),force=true)

    return

end

function templatehpc()

    src = joinpath(DEPOT_PATH[1],"files","SAMTools","hpc.txt")
    tem = joinpath(@__DIR__,"..","templates","hpc.txt")
    cp(tem,src,force=true)

    return

end

function overwritetemplate(
    overwriteexp :: Bool = false,
    overwriterun :: Bool = false,
    overwritehpc :: Bool = false,
)

    src = joinpath(DEPOT_PATH[1],"files","SAMTools","exp")
    if overwriteexp || !isdir(src)
        tem = joinpath(@__DIR__,"..","templates","exp")
        cp(tem,src,force=true)
    end

    src = joinpath(DEPOT_PATH[1],"files","SAMTools","run_template")
    if overwriterun || !isdir(src)
        tem = joinpath(@__DIR__,"..","templates","run_template")
        cp(tem,src,force=true)
        mkdir(joinpath(src,"LOGS"))
    end

    src = joinpath(DEPOT_PATH[1],"files","SAMTools","hpc.txt")
    if overwritehpc || !isdir(src)
        tem = joinpath(@__DIR__,"..","templates","hpc.txt")
        cp(tem,src,force=true)
    end

    return

end


text:
using MultivariateStats
using Base.Test

srand(34568)

## FA with zero mean

X = randn(5, 10)
Y = randn(3, 10)

W = qr(randn(5, 5))[1][:, 1:3]
Ψ = fill(0.1, 5)
M = FactorAnalysis(Float64[], W, Ψ)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == zeros(5)
@test loadings(M) == W
@test var(M) == Ψ

T = inv(I+W'*diagm(1./var(M))*W)*W'*diagm(1./var(M))
@test transform(M, X[:,1]) ≈ T * X[:,1]
@test transform(M, X) ≈ T * X

R = cov(M)*W*inv(W'W)
@test reconstruct(M, Y[:,1]) ≈ R * Y[:,1]
@test reconstruct(M, Y) ≈ R * Y

## PCA with non-zero mean

mv = rand(5)
M = FactorAnalysis(mv, W, Ψ)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == mv
@test loadings(M) == W
@test var(M) == Ψ

@test transform(M, X[:,1]) ≈ T * (X[:,1] .- mv)
@test transform(M, X) ≈ T * (X .- mv)

@test reconstruct(M, Y[:,1]) ≈ R * Y[:,1] .+ mv
@test reconstruct(M, Y) ≈ R * Y .+ mv

## prepare training data

d = 5
n = 1000

R = qr(randn(d, d))[1]
@test R'R ≈ eye(5)
scale!(R, sqrt.([0.5, 0.3, 0.1, 0.05, 0.05]))

X = R'randn(5, n) .+ randn(5)
mv = vec(mean(X, 2))
Z = X .- mv

# facm (default) & faem
fa_methods = [:cm, :em]

fas = FactorAnalysis[]
for method in fa_methods
    M = fit(FactorAnalysis, X, method=method, tot=5000)
    P = projection(M)
    W = loadings(M)
    push!(fas,M)

    @test indim(M) == 5
    @test outdim(M) == 4
    @test mean(M) == mv
    @test P'P ≈ eye(4)
    @test all(isapprox.(cov(M), cov(X, 2), atol=1e-3))

    M = fit(FactorAnalysis, X; mean=mv, method=method)
    @test loadings(M) ≈ W

    M = fit(FactorAnalysis, Z; mean=0, method=method)
    @test loadings(M) ≈ W

    M = fit(FactorAnalysis, X; maxoutdim=3, method=method)
    P = projection(M)

    @test indim(M) == 5
    @test outdim(M) == 3
    @test P'P ≈ eye(3)
end

# compare two algorithms
M1, M2 = fas
@test all(isapprox.(cov(M1), cov(M2), atol=1e-3)) # noise
LL(m, x) = (-size(x,2)/2)*(size(x,1)*log(2π) + log(det(cov(m))) + trace(inv(cov(m))*cov(x,2)))
@test LL(M1, X) ≈ LL(M2, X) # log likelihood

# test that fit works with Float32 values
X2 = convert(Array{Float32,2}, X)
# Float32 input
M = fit(FactorAnalysis, X2; method=:cm, maxoutdim=3)
M = fit(FactorAnalysis, X2; method=:em, maxoutdim=3)