text:
using CSV
using DataFrames
using Dates
using Statistics

function st_raw_data()
    # no alaska or hawaii, but includes district DC
    contstates = ("AL", "AZ" ,"AR" ,"CA" ,"CO" ,"CT" ,"DE" ,"DC" ,"FL" ,"GA", "ID" ,"IL" ,"IN" ,"IA" ,"KS" ,"KY" ,"LA" ,"ME" ,"MD" ,"MA" ,"MI" ,"MN" ,"MS" ,"MO" ,"MT" ,"NE" ,"NV" ,"NH" ,"NJ" ,"NM" ,"NY" ,"NC" ,"ND" ,"OH" ,"OK" ,"OR" ,"PA" ,"RI" ,"SC" ,"SD" ,"TN" ,"TX" ,"UT" ,"VT" ,"VA" ,"WA" ,"WV" ,"WI" ,"WY")

    # read metadata file for population of states and cities
    city_metadata = CSV.File("data/city_state_metadata.csv") |> DataFrame!
    filter!(row -> row[:abbr] in contstates, city_metadata)
    # convert string numbers to int numbers
    city_metadata[!, :statepop] = parse.(Int, (replace.(city_metadata.statepop, "," => "")))
    city_metadata[!, :citypop] = parse.(Int, (replace.(city_metadata.citypop, "," => "")))

    # read the demography of all states
    _st_demo = CSV.File("data/sc-est2019.csv") |> DataFrame!
    grps = (0:4, 5:16, 17:49, 50:65, 66:99)
    # select only the relevant columns select!(_st_demo, Not(8:17))
    # add age group column
    transform!(_st_demo, :AGE => ByRow(x -> findfirst(grp -> x in grp, grps)) => :AGEGRP)
    # filter to include both male/female, state
    filter!(r -> r[:SEX] == 0 && r[:AGE] <= 99, _st_demo)
    st_demo = combine(groupby(_st_demo, [:AGEGRP, :NAME]), :POPEST2019_CIV => sum)

    # append the demography to city_metadata
    transform!(city_metadata, :state => ByRow(x -> st_demo[st_demo.NAME .== x, :].POPEST2019_CIV_sum) => :age_dist)
    transform!(city_metadata, [:age_dist, :statepop] => ByRow((ag, stpop) -> ag ./ stpop .* 100000 ) => :age_dist_percapita) # convert to per 100000

    # read total vaccination coverage in the USA, need the total population vaccinated (to move to compartment V2) and current vaccination rate
    # this is the preexisting vaccination coverage, i.e. individuals in
    # file downloaded april 22, data as of apr 22 also
    total_vax = DataFrame!(CSV.File("data/cdc_totalvaccinations.csv", normalizenames=true, header=3))
    # total people vaccinated by percent, converted to 100,000 distributed over the age groups
    select!(total_vax, ["State_Territory_Federal_Entity", "People_Fully_Vaccinated_by_State_of_Residence", "Percent_of_Total_Pop_Fully_Vaccinated_by_State_of_Residence"])
    rename!(total_vax, [:state,:total_vaccinated, :percent_vaccinated])
    # fix cdc naming scheme for join
    replace!(total_vax.state, "New York State" => "New York")

    # join dataframes together
    final = leftjoin(city_metadata, total_vax, on = :state)
    # type stability of the joined columns
    final[!, :total_vaccinated] = convert.(Int64, final.total_vaccinated)
    final[!, :percent_vaccinated] = parse.(Float64,  final.percent_vaccinated) ./ 100
     # convert to per 100000
    final[!, :total_vaccination_percapita] = final.percent_vaccinated * 100000

    # use our world in data to get the current state vaccination rollout - i.e. doses per day
    # this gives us the dosage per day
    # file downloaded april 25, data as of april 25
    owid_dod = DataFrame!(CSV.File("data/owid_us_state_vaccinations.csv"))
    filter!(x -> x[:date] == maximum(owid_dod.date), owid_dod)
    select!(owid_dod, [:date, :location, :daily_vaccinations_per_million])
    owid_dod = coalesce.(owid_dod, 0.0) # replace missing with 0.0, note the broadcast
    owid_dod[!, :daily_vaccinations_per_million] = owid_dod.daily_vaccinations_per_million ./ 10 # to convert to per 100000
    final = leftjoin(final, owid_dod, on = (:state => :location))

    # get rvalues from epiforecasts_data downloaded april 25 https://epiforecasts.io/covid/posts/national/united-states/
    # file downlaoded april 25, but data is as of april 20
    rvals = DataFrame!(CSV.File("data/epiforecasts_data.csv", normalizenames=true))
    function _spltx(x)
        xplt = split(x);
        dd = [xplt[1], xplt[2][2:end], xplt[4][1:end-1]]
        parse.(Float64, dd)
    end
    transform!(rvals, :Effective_reproduction_no_ => ByRow(_spltx)  => [:rmean, :rlow, :rhi])
    select!(rvals, [:State, :rmean, :rlow, :rhi])
    #select!(rvals, Not(2:6)) # get rid of columns not need
    final = leftjoin(final, rvals, on = (:state => :State))

    # get the community spread transmission value for initial conditions
    # from the delphi covidcast https://delphi.cmu.edu/covidcast/   filename covidcast-fb.csv
    # data values are per 100 individuals
    # file downloaded april 25, but data is as of april 24 (see date column)
    com_spr = DataFrame!(CSV.File("data/covidcast_db.csv", header=2))
    com_spr.time_value = Date.(com_spr.time_value, "mm/dd/yyyy")  # convert to date column
    com_spr.value = com_spr.value .* 1000  # to get per 100000
    com_spr.geo_value = uppercase.(com_spr.geo_value) # for joining purposes
    filter!(x -> x[:time_value] == maximum(com_spr.time_value), com_spr) # get latest date data

    final = leftjoin(final, com_spr, on=(:abbr => :geo_value))

    # get the total number of positive (confirmed + probable cases) as initial herd immunity
    # file downloaded april 26, data as of march 7... covidtrackingproject stopped on march 7
    inithi = DataFrame!(CSV.File("data/covidtrackingproj_stdata.csv"))
    filter!(x -> x[:date] == maximum(inithi.date), inithi) # get latest date data
    select!(inithi, [:state, :positive]) # only really need the probable+confirmed cases for herd immunity
    inithi[!, :positive] = convert.(Int64, inithi.positive) # type stable
    final = leftjoin(final, inithi, on=(:abbr => :state))

    return final
end
const st_data = st_raw_data()

struct state_information8
    nme::Symbol             # name of state
    pop::Vector{Float64}    # population distribution per 100,000
    inf::Float64   # initial infected population per 100,100
    pre::Float64            # level of presymptomatic (in percent)
    dpd::Float64            # vaccine dpd per 100,000
    dis::NTuple{5, Float64} # distribution of dpd over age groups
    ref::Float64            # r Effective
    #new(:none, zeros(Float64, 5), 0.0, 0.0, 0.0, (0.0, 0.0, 0.3138, 0.2408, 0.4478), 0.0)
    #cov = params.kffcov .* params.pop
    # params.stinit = true
end

function state_information8(st)
    frow = filter(x -> x[:state] == st, st_data)
    nrow(frow) != 1 && error(" $(nrow(frow)) states found")
    stinfo = state_information8(Symbol(st), frow.age_dist_percapita[1], frow.value[1], Float64(frow.positive[1]), Float64(frow.daily_vaccinations_per_million[1]), (0.0, 0.0, 0.3138, 0.2408, 0.4478), frow.rmean[1])
    frow.age_dist_percapita[1]
    return stinfo
end


text:
# Implementation of GeoInterface for Base Types

GeoInterface.isgeometry(::Type{<:AbstractVector{<:Real}}) = true
GeoInterface.geomtrait(::AbstractVector{<:Real}) = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::AbstractVector{<:Real}) = Base.length(geom)
GeoInterface.getcoord(::PointTrait, geom::AbstractVector{<:Real}, i) = getindex(geom, i)

GeoInterface.isgeometry(::Type{<:NTuple{N,<:Real}}) where {N} = true
GeoInterface.geomtrait(::NTuple{N,<:Real}) where {N} = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::NTuple{N,<:Real}) where {N} = N
GeoInterface.getcoord(::PointTrait, geom::NTuple{N,<:Real}, i) where {N} = getindex(geom, i)

for i in 2:4
    sig = NamedTuple{default_coord_names[1:i],NTuple{i,T}} where {T<:Real}
    GeoInterface.isgeometry(::Type{<:sig}) = true
    GeoInterface.geomtrait(::sig) = PointTrait()
    GeoInterface.ncoord(::PointTrait, geom::sig) = i
    GeoInterface.getcoord(::PointTrait, geom::sig, i) = getindex(geom, i)
end

# Custom coordinate order/names NamedTuple
GeoInterface.isgeometry(::Type{<:NamedTuple{Keys,NTuple{N,T}}}) where {Keys,N,T<:Real} = all(in(default_coord_names), Keys)
GeoInterface.geomtrait(::NamedTuple{Keys,NTuple{N,T}}) where {Keys,N,T<:Real} = PointTrait()
GeoInterface.ncoord(::PointTrait, geom::NamedTuple{Keys,NTuple{N,T}}) where {Keys,N,T<:Real} = Base.length(geom)
GeoInterface.getcoord(::PointTrait, geom::NamedTuple{Keys,NTuple{N,T}}, i) where {Keys,N,T<:Real} = getindex(geom, i)
GeoInterface.coordnames(::PointTrait, geom::NamedTuple{Keys,NTuple{N,T}}) where {Keys,N,T<:Real} = Keys

# Default features using NamedTuple and AbstractArray

const NamedTupleFeature = NamedTuple{(:geometry,:properties)}

GeoInterface.isfeature(::Type{<:NamedTupleFeature}) = true
GeoInterface.trait(::NamedTupleFeature) = FeatureTrait()
GeoInterface.geometry(f::NamedTupleFeature) = f.geometry
GeoInterface.properties(f::NamedTupleFeature) = f.properties

const ArrayFeatureCollection = AbstractArray{<:NamedTupleFeature}

GeoInterface.isfeaturecollection(::Type{<:ArrayFeatureCollection}) = true
GeoInterface.trait(::ArrayFeatureCollection) = FeatureCollectionTrait()
GeoInterface.nfeature(::FeatureCollectionTrait, fc::ArrayFeatureCollection) = Base.length(fc)
GeoInterface.getfeature(::FeatureCollectionTrait, fc::ArrayFeatureCollection, i::Integer) = fc[i]

ext: jl


text:
"""
```
function my_speye(n::Integer)
```
Sparse identity matrix - since deprecated in 0.7
"""
@inline function my_speye(n::Integer)
    return SparseMatrixCSC{Float64}(LinearAlgebra.I, n, n)
end

"""
```
function my_speye(T::Type, n::Integer)
```
Sparse identity matrix - since deprecated in 0.7
"""
@inline function my_speye(T::Type, n::Integer)
    return SparseMatrixCSC{T}(LinearAlgebra.I, n, n)
end

"""
```
function stat_dist(λ::AbstractArray{T}) where {T<:Float64}
```
Computes stationary distribution.
"""
function stat_dist(λ::AbstractArray{T}) where {T<:Float64}
    F     = eigen(copy(λ))
    eigvl = F.values
    eigvc = F.vectors
    pos   = argmin(abs.(eigvl))
    statdist = eigvc[:,pos]
    return statdist ./ sum(statdist)
end

"""
```
@inline function changeBasis(basis,inv_basis,g1,psi,pi,c,g0)
```
Update dynamics with new basis.
(by SeHyoun Ahn, June 2016)

"""
@inline function changeBasis(basis,inv_basis,g1,psi,pi,c,g0)
    g1 = basis*g1*inv_basis
    pi = basis*pi
    psi = basis*psi
    c = basis*c
    if (nargin==7)
        g0 = basis*g0*inv_basis
    else
        g0 = my_speye(size(g1, 1))
    end
    return g1, psi, pi, c, g0
end

"""
```
function stateSpaceReduction(g0, g1, n_v, n_g, n_p, n_Z, reduceDist_hor)
```
Reduce dimensionality of state space.

# Arguments
- `g0`: LHS matrix, only used to check it satisfies require form
- `g1`: Dynamics matrix
- `n_v::Int64`: number of jump variables
- `n_g::Int64`: number of state variables
- `n_p::Int64`: number of static constraints
- `n_Z::Int64`: number of exogenous variables

# Output
- `state_red`: tranformation to get from full grid to reduced states
- `inv_state_red`: inverse transform
- `n_g`: number of state variables after reduction
"""
function stateSpaceReduction(g0, g1, n_v, n_g, n_p, n_Z, reduceDist_hor)
    # Check to make sure that LHS satisfies the required form
    loc = [1:n_v+n_g, n_v+n_g+n_p+1:n_v+n_g+n_p+n_Z]
    #if (maximum(abs.(g0[loc,loc] .- my_speye(n_v+n_g+n_Z)))≈0)
    #    error("Make sure that g0 is normalized.")
    #end

    n_total = n_v + n_g

    # Slice Dynamics Equation into Different Parts
    B_pv = g1[n_v+n_g+1:n_v+n_g+n_p,n_v+n_g+1:n_v+n_g+n_p] \ g1[n_v+n_g+1:n_v+n_g+n_p,1:n_v]
    B_pZ = g1[n_v+n_g+1:n_v+n_g+n_p,n_v+n_g+1:n_v+n_g+n_p] \ g1[n_v+n_g+1:n_v+n_g+n_p,n_v+n_g+n_p+1:n_v+n_g+n_p+n_Z]
    B_pg = g1[n_v+n_g+1:n_v+n_g+n_p,n_v+n_g+1:n_v+n_g+n_p] \ g1[n_v+n_g+1:n_v+n_g+n_p,n_v+1:n_v+n_g]
    B_gg = g1[n_v+1:n_v+n_g,n_v+1:n_v+n_g]
    B_gp = g1[n_v+1:n_v+n_g,n_v+n_g+1:n_v+n_g+n_p]

    # Orthonormalize B_pg as the direction to keep
    # Drop redundant Directions
    ~, d0, V_g = svd(B_pg)#,'econ')
    aux = diagm(0 => d0)
    n_Bpg = Int64(sum(aux .> 10 * eps() * aux[1]))
    V_g = V_g[:, 1:n_Bpg] .* aux[1:n_Bpg]'

    ## Arnoldi Iteration
    hor = reduceDist_hor
    A(x::Matrix{Float64}) = B_gg' * x - B_pg' * (B_gp'*x)
    # [V_g,~] = block_arnoldi_func(A,V_g,hor)
    V_g, ~ = deflated_block_arnoldi(A, V_g, hor)
    n_g = size(V_g, 2)

    ## Build State-Space Reduction transform
    state_red = spzeros(Float64, n_v + n_g + n_Z, n_total + n_p + n_Z)
    state_red[1:n_v,1:n_v] = my_speye(n_v)
    state_red[n_v+1:n_v+n_g,n_v+1:n_total] = V_g'
    state_red[:,n_total+1:n_total+n_p] = 0
    state_red[n_v+n_g+1:n_v+n_g+n_Z,n_total+n_p+1:n_total+n_p+n_Z] = my_speye(n_Z)

    ## Build inverse transform
    inv_state_red = spzeros(Float64, n_total+n_p+n_Z,n_v+n_g+n_Z)
    inv_state_red[1:n_v,1:n_v] = my_speye(n_v)
    inv_state_red[n_v+1:n_total,n_v+1:n_g+n_v] = V_g
    inv_state_red[n_total+1:n_total+n_p,n_v+1:n_v+n_g] = -B_pg*V_g
    inv_state_red[n_total+1:n_total+n_p,1:n_v] = -B_pv
    inv_state_red[n_total+1:n_total+n_p,n_v+n_g+1:n_v+n_g+n_Z] = -B_pZ
    inv_state_red[n_total+n_p+1:n_total+n_p+n_Z,n_v+n_g+1:n_v+n_g+n_Z] = my_speye(n_Z)

    return state_red, inv_state_red, n_g
end

"""
```
function cleanG0Sparse(g0, g1, c, pi, psi)
```
Inverts g0. For certain form on non-invertible g0 coming from constraints,
it will compute the proper inverse.

(by SeHyoun, June 2016)
"""
function cleanG0Sparse(g0, g1, c, pi, psi)
    n = size(g0,1)

    tmp = (max(abs(g0),[],2)==0)

    redundant = findall(tmp)
    keep=findall(1 .- tmp)
    n_keep = length(keep)

    base = sparse(n,n_keep)
    base[keep,:] = my_speye(n_keep)

    base[redundant,:] = -g1[redundant,redundant] \ g1[redundant,keep]

    inv_base = sparse(n_keep,n)
    inv_base[:,keep] = my_speye(n_keep)

    g0  = inv_base * g0 * base
    g1  = inv_base * g1 * base
    g1  = g0 \ g1
    psi = g0 \ inv_base * psi
    pi  = g0 \ inv_base * pi
    c   = g0 \ inv_base * c

    return base, inv_base, g1, c, pi, psi
end

function covv(x,y)
    cov_xy = cov(x,y)
    cov_xy = cov_xy[1,2]
    return cov_xy
end

function Gini_fn(gg, c, dab_aux)
    S_c  = cumsum(gg .* dab_aux .* c[:])/sum(gg .* dab_aux .* c[:])
    trapez_c = 0.5*(S_c[1]*gg[1]*dab_aux[1] +
                      sum((S_c[2:end] + S_c[1:end-1]).*gg[2:end].*dab_aux[2:end]))
    C_Gini = 1 - 2*trapez_c
    return C_Gini
end

"""
```
function gstate(G1, impact; pick = Matrix{Float64}(undef, 0, 0))
```
#function  [GS,pickn,ok,uis,vs]=gstate(G1,impact,pick)
# Arguments
- `G1`:     the coefficient on  lagged y from the output of gensys
- `impact`: the coefficient on exogenous shocks from the output of gensys
- `pick`:   an optional guess at a matrix of coefficients that extracts the state vector from y

# Output
- `ok`: 0 if pick matrix is just no good.
    1 if pick matrix is usable forward, but loses initial date information
    2 if pick matrix is not usable forward, is part of a correct state vector, but is not complete
    3 if pick matrix is     usable forward, is part of a correct state vector, but is not complete
    4 if pick matrix is not usable forward, summarizes past, but is redundant
    5 if pick matrix is     usable forward, summarizes past, but is redundant
    6 if pick matrix summarizes past and is not redundant, but is not usable forward
    7 if pick matrix is perfect, both forward and as history summary

- `pickn`: a matrix of coefficients that extracts the state vector from y.  Equal
    to pick if pick is supplied and ok=1; otherwise, an ok-maximizing pick matrix that
    is somewhat similar to pick.

- `GS`: the matrix of coefficients on the lagged state variable
- `uis,vs`: If uis'*vs is square and full rank, ok=7 is possible, otherwise not.
    If ok<2, an ok>2 can be found by trying a pick in the row space of vs'.
    Any pick with pick*uis full column rank will provide a foward state (i.e. ok an odd number).
- `uiu`:   uiu'y(t)=0, and this is interpretable as the decision rule setting controls
    as functions of the states.

The solution was in the form y(t)=G1*y(t-1)+impact*z(t).
Now it's in the form y(t)=GS*pickn*y(t-1)+impact*z(t).
In general, `pickn` is mxn with m<n.
"""
function gstate(G1, impact; pick = Matrix{Float64}(undef, 0, 0))
    REALSMALL = 1e-9
    nr, nc = size(G1)
    F = svd(G1)
    u, d, v = F.U, F.S, F.V #[u d v] = svd(G1)
    top = findall(d .> REALSMALL) #diag(d)>REALSMALL)
    nd = top[end]
    us = u[:,top]
    uu = u[:,top+1:end]
    d  = d[top,top]
    vs = v[:,top]
    if length(pick) == 0 #nargin <= 2
        pick = vs'
        vp   = v
        vps  = vs
        dp   = Matrix{Float64}(I,nd)
        ups  = Matrix{Float64}(I,nr,nd)
        ndp  = nd
    else
        Fp = svd(pick)
        up, dp, vp = Fp.U, Fp.S, Fp.V #[up dp vp] = svd(pick)
        topp = findall(dp .> REALSMALL)
        ndp  = topp[end]
        vps  = vp[:,topp]
        dp   = dp[topp,topp]
        ups  = up[:,topp]
    end

    # If we were worried about efficiency, we'd skip some of this when pick = vs.
    #Does pick summarize history?  (pick in v' row space, v in vp' row space).
    pinv = all(all((pick'-vs*vs'*pick').^2<REALSMALL))
    vinp = all(all((vs-vps*vps'*vs).^2<REALSMALL))
    okpast = pinv+vinp
    #  Does pick summarize all current info?  (impact in us column space, pick*uu full rank)
    Fi = svd(hcat(impact, us)) # RECA: [impact us]
    ui, di, vi = Fi.U, Fi.S, Fi.V
    topi = findall(di .> REALSMALL)
    ndi = length(topi)
    uis = ui[:,topi]
    uiu = ui[:,ndi+1:size(ui,2)]
    if ndi < size(G1,1)
        if size(pick,1) < size(uis,2)
            oknow = 0
        else
            Ft = svd(pick*uis)
            ut, dt, vt = Ft.U, Ft.S, Ft.V
            toppu = findall(dt .> REALSMALL)
            if length(toppu) < size(us,2)
                oknow = 0
            else
                oknow = 1
            end
        end
    else
        oknow = 0
    end
    if vinp
        GS = G1/pick
        pickn = pick
    elseif pinv
        r = vs-vps*vps'*vs
        Fr = svd(r)
        ur, dr, vr = Fr.U, Fr.S, Fr.V
        topr = findall(dr .> REALSMALL)
        p2 = ur[:,topr]
        pickn = [pick;p2']
        GS = G1/pickn
    elseif oknow
        GS = G1/[pick;uiu']
        GS = GS[:,1:size(pick,1)]
        pickn = pick
    else
        pickn = vs'
        GS = us*d
    end
    ok = oknow + 2*pinv + 4*vinp
    return GS, pickn, ok, uis, uiu, vs
end


text:
using StaticArrays

include("partialbridge!.jl")

T = 2.0
dt = 1/1000

tt = 0.:dt:T
#struct PBIntegratedDiffusion <: ContinuousTimeProcess{ℝ{2}}
#    γ::Float64
#end

Bridge.b(t::Float64, x, P::PBIntegratedDiffusion) = ℝ{2}(x[2], βu(t, x[2], P))

@test Bridge.b!(0.0, X1.yy[1], copy( X1.yy[1]),  P) == Bridge.b(0.0, X1.yy[1], P)

Bridge.σ(t, x, P::PBIntegratedDiffusion) = ℝ{2}(0.0, P.γ)

@test Bridge.σ!(0.0, X1.yy[1], 1.0, copy( X1.yy[1]),  P) == Bridge.σ(0.0, X1.yy[1], P)

Bridge.constdiff(::PBIntegratedDiffusion) = true

#struct PBIntegratedDiffusionAux <: ContinuousTimeProcess{ℝ{2}}
#    γ::Float64
#end

Bridge.b(t::Float64, x, P::PBIntegratedDiffusionAux) = ℝ{2}(x[2], βu(t, x[2], P))

@test Bridge.b!(0.0, X1.yy[1], copy( X1.yy[1]),  Pt) == Bridge.b(0.0, X1.yy[1], Pt)

Bridge.σ(t, P::PBIntegratedDiffusionAux) =  ℝ{2}(0.0, P.γ)
Bridge.σ(t, x, P::PBIntegratedDiffusionAux) = Bridge.σ(t, P)

Bridge.B(t, P::PBIntegratedDiffusionAux) = @SMatrix [0.0 1.0; 0.0 -1.0]
Bridge.β(t, P::PBIntegratedDiffusionAux) = ℝ{2}(0, 1/2)
Bridge.a(t, P::PBIntegratedDiffusionAux) = @SMatrix [0.0 0.0; 0.0 P.γ^2]

# Generate Data
Random.seed!(1)

W2 = sample(tt, Wiener())

@test W2.yy == W1.yy

X2 = solve(Euler(), ℝ{2}(x0), W2, P)

@test X2.yy == X1.yy

L = @SMatrix [1. 0.]
Σ = @SMatrix [0.01]
v = ℝ{1}(2.5)

ϵ = 0.02
# Solve Backward Recursion

S2 = typeof(L)
S = typeof(L*L')
T = typeof(diag(L*L'))

N = length(tt)
Lt = zeros(S2, N)
Ht = zeros(S, N)
νt = zeros(T, N)

Po2 = Bridge.PartialBridgeνH(tt, P, Pt, L, v, ϵ, Σ)

Xo2 = copy(X2)
solve!(Euler(), Xo2, ℝ{2}(x0), W2, Po2)

@test norm(Xo1.yy - Xo2.yy) < 500*eps()

# Likelihood

ll2 = llikelihood(Bridge.LeftRule(), Xo2, Po2)

@test ll1 ≈ ll2

lls2 = Float64[]

    # MCMC parameter
function mcmc2(x0, tt, Po2)
    iterations = 1000
    subsamples = 0:100:iterations
    ρ = 0.9

    # initalization
    W = sample(tt, Wiener())
    X = solve(Euler(), x0, W, Po2)
    ll = llikelihood(Bridge.LeftRule(), X, Po2)

    acc = 0

    Wo = copy(W)
    Xo = copy(X)
    Wrho = copy(W)

    XX = Any[]
    if 0 in subsamples
        push!(XX, copy(X))
    end
    for iter in 1:iterations
        # Proposal
        sample!(Wrho, Wiener())
        Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*Wrho.yy

        solve!(Euler(), Xo, x0, Wo, Po2)
        llo = llikelihood(Bridge.LeftRule(), Xo, Po2)
        push!(lls2, llo)

        if log(rand()) <= llo - ll
            X, Xo = Xo, X
            W, Wo = Wo, W

            ll = llo
            acc += 1
        end
        if iter in subsamples
            push!(XX, copy(X))
        end
    end

end

@time @testset "MCMC" begin
    mcmc2(ℝ{2}(x0), tt, Po2)

    @test norm(lls1 - lls2) < sqrt(eps())

end


text:
module CartesianDDM

using Base.Iterators
using LinearAlgebra
using OffsetArrays

import Base: OneTo
import Base: parent, size, getindex, setindex!, similar
import Base: *, \

import LinearAlgebra: ldiv!

export CartesianDDMDimension
export CartesianDDMContext
export CartesianDDMVector
export CartesianDDMRAS

export decompose

include("decomposition.jl")
include("arrays.jl")
include("vectors.jl")
include("matrices.jl")
include("coherence.jl")
include("preconditioner.jl")
include("ras.jl")

end


text:
# N stochastic A->B reactions at different rates
using DiffEqBase, DiffEqJump, DiffEqJumpExtensions
using Base.Test

doprint = true

function A_to_B_mean(N)
    Nsims    = 32000
    tf       = .1
    baserate = .1
    A0       = 100

    rates    = ones(Float64, N) * baserate;
    cumsum!(rates, rates)
    exactmean = (t) -> A0*exp(-sum(rates) * t)

    # jump reactions
    jumpvec = []
    for i in 1:N
        ratefunc = (u,p,t) -> rates[i] * u[1]
        affect!  = function (integrator)
            integrator.u[1] -= 1
            integrator.u[2] += 1
        end
        push!(jumpvec, ConstantRateJump(ratefunc, affect!))
    end

    # convert jumpvec to tuple to send to JumpProblem...
    jumps = ((jump for jump in jumpvec)...)

    prob = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, DiffEqJumpExtensions.Direct(), save_positions=(false,false), jumps...)

    Asamp = zeros(Int64,Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Asamp[i] = sol[1,end]
    end

    if doprint
        println("samp mean: ", mean(Asamp), ", act mean = ", exactmean(tf))
    end

    @test abs(mean(Asamp) - exactmean(tf)) < 1.

end

# tuples
A_to_B_mean(5)

# function wrappers
A_to_B_mean(15)

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
Base.@kwdef mutable struct CollectedDiagnostics
    onetime_done::Bool = false
    ΣMH_z::Union{Nothing, Array} = nothing
end
const Collected = CollectedDiagnostics()

function collect_onetime(mpicomm, dg, Q)
    if !Collected.onetime_done
        FT = eltype(Q)
        grid = dg.grid
        grid_info = basic_grid_info(dg)
        topl_info = basic_topology_info(grid.topology)
        topology = grid.topology
        Nq1 = grid_info.Nq[1]
        Nq2 = grid_info.Nq[2]
        Nq3 = grid_info.Nqk
        npoints = prod(grid_info.Nq)
        nrealelem = topl_info.nrealelem
        nvertelem = topl_info.nvertelem
        nhorzelem = topl_info.nhorzrealelem

        vgeo = array_device(Q) isa CPU ? grid.vgeo : Array(grid.vgeo)

        Collected.ΣMH_z = zeros(FT, Nq3, nvertelem)

        for eh in 1:nhorzelem, ev in 1:nvertelem
            e = ev + (eh - 1) * nvertelem
            for k in 1:Nq3, j in 1:Nq2, i in 1:Nq1
                ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
                MH = vgeo[ijk, grid.MHid, e]
                Collected.ΣMH_z[k, ev] += MH
            end
        end

        # compute the full number of points on a slab on rank 0
        MPI.Reduce!(Collected.ΣMH_z, +, 0, mpicomm)

        Collected.onetime_done = true
    end

    return nothing
end


text:
#=
Sudoku is a puzzle where you're given a partially-filled 9 by 9 grid with digits. The objective is to fill the grid with the constraint that every row, column, and box (3 by 3 subgrid) must contain all of the digits from 1 to 9.

Implement an efficient sudoku solver.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
##########
# Addition
##########

"""
    Base.:+(a::AbstractUncertainValue, b::AbstractUncertainValue) -> UncertainValue

Addition operator for pairs of uncertain values.

Computes the element-wise sum between for a default of `n = 10000` realizations of `a` and
`b`, then returns an uncertain value based on a kernel density estimate to the distribution
of the element-wise sums.

Use the `+(a, b, n)` syntax to tune the number (`n`) of draws.
"""
function Base.:+(a::AbstractUncertainValue, b::AbstractUncertainValue; n::Int = 30000)
    UncertainValue(resample(a, n) .+ resample(b, n))
end

"""
    Base.:+(a::Real, b::AbstractUncertainValue) -> UncertainValue

Addition operator for between scalars and uncertain values.

Computes the element-wise sum between `a` and `b` for a default of `n = 10000` realizations
of `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise sums.

Use the `+(a, b, n)` syntax to tune the number (`n`) of draws.
"""
Base.:+(a::Real, b::AbstractUncertainValue; n::Int = 30000) =
    UncertainValue(a .+ resample(b, n))

"""
    Base.:+(a::AbstractUncertainValue, b::Real) -> UncertainValue

Addition operator for between uncertain values and scalars.

Computes the element-wise sum between `a` and `b` for a default of `n = 10000` realizations
of `a`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise sums.

Use the `+(a, b, n)` syntax to tune the number (`n`) of draws.
"""
Base.:+(a::AbstractUncertainValue, b::Real; n::Int = 30000) =
    UncertainValue(resample(a, n) .+ b)

"""
    Base.:+(a::AbstractUncertainValue, b::AbstractUncertainValue, n::Int) -> UncertainValue

Addition operator for pairs of uncertain values.

Computes the element-wise sum between `a` and `b` for `n` realizations
of `a` and `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise sums.

Call this function using the `+(a, b, n)` syntax.
"""
function Base.:+(a::AbstractUncertainValue, b::AbstractUncertainValue, n::Int)
    UncertainValue(resample(a, n) .+ resample(b, n))
end
"""
    Base.:+(a::Real, b::AbstractUncertainValue, n::Int) -> UncertainValue

Addition operator for scalar-uncertain value pairs.

Computes the element-wise sum between `a` and `b` for `n` realizations
of `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise sums.

Call this function using the `+(a, b, n)` syntax.
"""
Base.:+(a::Real, b::AbstractUncertainValue, n::Int) =
    UncertainValue(a .+ resample(b, n))

"""
    Base.:+(a::AbstractUncertainValue, b::Real, n::Int) -> UncertainValue

Addition operator for scalar-uncertain value pairs.

Computes the element-wise sum between `a` and `b` for `n` realizations
of `a`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise sums.

Call this function using the `+(a, b, n)` syntax.
"""
Base.:+(a::AbstractUncertainValue, b::Real, n::Int) =
    UncertainValue(resample(a, n) .+ b)

#######################################
# Adding vectors of uncertain values
#######################################
function Base.:+(a::AbstractVector{AbstractUncertainValue},
        b::AbstractVector{AbstractUncertainValue})
    [a[i] + b[i] for i = 1:length(a)]
end

function Base.:+(a::AbstractVector{AbstractUncertainValue}, b::AbstractUncertainValue)
    [a[i] + b for i = 1:length(a)]
end

function Base.:+(a::AbstractUncertainValue, b::AbstractVector{AbstractUncertainValue})
    [a + b[i] for i = 1:length(b)]
end

function Base.:+(a::AbstractVector{AbstractUncertainValue},
        b::AbstractVector{AbstractUncertainValue}, n::Int)

    [+(a[i], b[i], n) for i = 1:length(a)]
end

function Base.:+(a::AbstractVector{AbstractUncertainValue},
        b::AbstractUncertainValue, n::Int)

    [+(a[i], b, n) for i = 1:length(a)]
end

function Base.:+(a::AbstractUncertainValue,
        b::AbstractVector{AbstractUncertainValue}, n::Int)

    [+(a, b[i], n) for i = 1:length(b)]
end

#####################################################################################
# Special cases
#####################################################################################

import ..UncertainValues: CertainValue

##################
# `CertainValue`s
#################
"""
    Base.:+(a::Union{CertainValue, Real}, b::Union{CertainValue, Real})

Addition of certain values with themselves or scalars acts as regular addition, but
returns the result wrapped in a `CertainValue` instance.
"""
Base.:+(a::Union{CertainValue, Real}, b::Union{CertainValue, Real})

Base.:+(a::CertainValue, b::CertainValue) = CertainValue(a.value + b.value)
Base.:+(a::CertainValue, b::Real) = CertainValue(a.value + b)
Base.:+(a::Real, b::CertainValue) = CertainValue(a + b.value)


text:
#=
Reader for cell/param files
=#
using Unitful
using UnitfulAtomic

"Clean lines by removing new line symbols and comments"
function clean_cell_lines(lines_in::Vector{String})
    lines_out = String[]
    for line in lines_in
        tmp = strip(line)
        length(tmp) == 0 && continue
        # Location of the comment
        matched = match(r"[#!]", tmp, 1)
        # No comment - keep the whole line
        if matched === nothing
            push!(lines_out, string(tmp))
        elseif matched.offset == 1
            continue
        else
            # Remove the trailing part
            push!(lines_out, string(strip(tmp[1:matched.offset-1])))
        end
    end
    lines_out
end

"Read in CASTEP's cell/param files"
function read_cell_castep(fname)
    f = open(fname)
    lines = clean_cell_lines(readlines(f))
    close(f)
    kws = Dict{Symbol, String}()
    blocks = Dict{Symbol, Vector{String}}()
    inblock = false
    block_name = ""
    current_block = String[]
    for il in eachindex(lines)
        line = lowercase(lines[il])
        # Empty line - do nothing
        length(il) == 0 && continue

        # Are we hitting a block
        mblock = match(r"%block (\w+)", line)
        if mblock !== nothing
            block_name = mblock.captures[1]
            inblock = true
            continue
        end

        mblock = match(r"%endblock (\w+)", line)
        if mblock !== nothing
            block_name = mblock.captures[1]
            inblock = false
            # Record the block
            blocks[Symbol(block_name)] = current_block
            # Reset current block
            current_block = String[]
            block_name = ""
            continue
        end

        # In inside the block - push the lines
        if inblock == true
            push!(current_block, string(strip(line)))
            continue
        end

        # Process keyword argumetns
        tokens = split(line, r"[:= ]+", limit=2)
        if length(tokens) == 1
            kws[Symbol(tokens[1])] = ""
        else
            kws[Symbol(tokens[1])] = string(tokens[2])
        end
    end
    CastepInFile(kws, blocks)
end

"Type for CASTEP input file"
struct CastepInFile
    kws::Dict
    blocks::Dict
end

function Base.show(io::IO, c::CastepInFile)
    println("CastepInFile with keys: ", join(keys(c.kws), ", "))
    println("Blocks: ", join(keys(c.blocks), ", "))
end

"Allow access fields by c[:xx] notation"
function Base.getindex(c::CastepInFile, s::Symbol)
    if (s in keys(c.kws)) & (s in keys(c.blocks))
        throw(ErrorException("Key $(s) exists in both key and blocks"))
    end
    if s in keys(c.kws)
        return c.kws[s]
    end
    if s in keys(c.kws)
        return c.blocks[s]
    end
    throw(KeyError(s))
end

"""Read out the lattice from the input file"""
function read_lattice(ci::CastepInFile)
    cellmat = zeros(Float64, 3, 3)
    unit = ""
    if :lattice_cart in keys(ci.blocks)
        i = 1
        for line in ci.blocks[:lattice_cart]
            tokens = split(line)
            if length(tokens) == 1
                unit = tokens[1]
                continue
            end
            cellmat[:, i] = [parse(Float64, s) for s in tokens]
            i += 1
        end
    elseif :lattice_abc in keys(ci.blocks)
        a, b, c = [parse(Float64, x) for x in split(ci.blocks[:lattice_abc][1])]
        α, β, γ = [parse(Float64, x) for x in split(ci.blocks[:lattice_abc][2])]
        cellmat = cellpar2mat(a, b, c, α, β, γ, degree=true)
    else
        throw(ErrorException("No lattice related block is found"))
    end
    cellmat * u"Å"
end

"""
Read in the symmetry operations

Returns matrix of the rotations and that for the displacements
"""
function read_symm_disp(ci::CastepInFile)
    lines = ci.blocks[:symmetry_ops]
    vecs = []
    for line in lines
        push!(vecs, [parse(Float64, x) for x in split(line)])
    end
    nvec = length(vecs)
    data = reshape(hcat(vecs...), 3, 4, div(nvec, 4))
    # Separate rotational symmetries and the displacement
    data[:, 1:3, :], data[:, 4, :]
end


text:
function predict_new_rbf(a_all,a_layer1,x_new,obs,neurons,n_per_part,inds_all,x_old,vars,cc)
    println("hello")
    layer1=zeros(obs,neurons)
    for i=1:neurons
        if i in 1:1:neurons println(Dates.format(now(), "HH:MM:SS")," Predicting responce of ",i," neuron, of ",neurons," Total.") end
        x1=copy(x_old[inds_all[n_per_part[i]+1:n_per_part[i+1]],:])
        obs1=size(x1,1)
        phi1 = calc_phi(x1,x_new,obs1,obs,cc,vars)
        # layer1[:,i].=phi1*a_all[i,1:obs1]
        layer1[:,i].=phi1*a_all[i]

    end

    predl1=[layer1 ones(obs)]*a_layer1

    return predl1,layer1

end


text:
"""
    Neighborhood

Neighborhoods define the pattern of surrounding cells in the "neighborhood"
of the current cell. The `neighbors` function returns the surrounding
cells as an iterable.

The main kinds of neighborhood are demonstrated below:

![Neighborhoods](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/Neighborhoods.png)

Neighborhoods can be used in [`NeighborhoodRule`](@ref) and [`SetNeighborhoodRule`](@ref) -
the same shapes with different purposes. In a `NeighborhoodRule` the neighborhood specifies
which cells around the current cell are returned as an iterable from the `neighbors` function.
These can be counted, summed, compared, or multiplied with a kernel in an
`AbstractKernelNeighborhood`, using [`kernelproduct`](@ref).

In `SetNeighborhoodRule` neighborhoods give the locations of cells around the central cell,
as [`offsets`] and absolute [`positions`](@ref) around the index of each neighbor. These
can then be written to manually.
"""
abstract type Neighborhood{R,L} end

ConstructionBase.constructorof(::Type{<:T}) where T <: Neighborhood{R,L} where {R,L} =
    T.name.wrapper{R,L}

radius(hood::Neighborhood{R}) where R = R
neighbors(hood::Neighborhood{R}) where R = map(i -> _buffer(hood)[i], bufindices(hood))
_buffer(hood::Neighborhood) = hood._buffer
@inline positions(hood::Neighborhood, I) = map(o -> o .+ I, offsets(hood))
@inline function bufindices(hood::Neighborhood{R}) where R
    map(offsets(hood)) do o
        i = o[1] + (R + 1)
        j = o[2] + (R + 1)
        i + (j - 1) * (2R + 1)
    end
end

Base.eltype(hood::Neighborhood) = eltype(_buffer(hood))
Base.length(hood::Neighborhood{<:Any,L}) where L = L
Base.iterate(hood::Neighborhood, args...) = iterate(neighbors(hood), args...)
Base.getindex(hood::Neighborhood, i) = getindex(_buffer(hood), bufindices(hood)[i])

"""
    RadialNeighborhood <: Neighborhood

Square neighborhoods with radius `R`, and side length `2R + 1`
"""
abstract type RadialNeighborhood{R,L} <: Neighborhood{R,L} end

"""
    Moore <: RadialNeighborhood

    Moore(radius::Int=1)

Moore neighborhoods define the neighborhood as all cells within a horizontal or
vertical distance of the central cell. The central cell is omitted.
"""
struct Moore{R,L,B} <: RadialNeighborhood{R,L}
    _buffer::B
end
# Buffer is updated later during the simulation.
Moore(radius::Int=1) = Moore{radius}()
Moore{R}(_buffer=nothing) where R = Moore{R,(2R+1)^2-1}(_buffer)
Moore{R,L}(_buffer::B=nothing) where {R,L,B} = Moore{R,L,B}(_buffer)

@generated function bufindices(hood::Moore{R}) where R
    buflen = (2R + 1)^2 # Use linear indexing
    centerpoint = buflen ÷ 2 + 1
    exp = Expr(:tuple)
    for i in 1:buflen
        if i != centerpoint
            push!(exp.args, :($i))
        end
    end
    return exp
end
@generated function offsets(hood::Moore{R}) where R
    exp = Expr(:tuple)
    for j in -R:R, i in -R:R
        if i != (0, 0)
            push!(exp.args, :($i, $j))
        end
    end
    return exp
end
@inline _setbuffer(n::Moore{R,L}, buf::B2) where {R,L,B2} = Moore{R,L,B2}(buf)

# Neighborhood specific `sum` for performance
Base.sum(hood::Moore) = sum(_buffer(hood)) - _centerval(hood)

_centerval(hood::Neighborhood{R}) where R = _buffer(hood)[R + 1, R + 1]

"""
    Window <: RadialNeighborhood

    Window{R}()

A neighboorhood of radius R that includes the central cell.
`R = 1` gives a 3x3 matrix.
"""
struct Window{R,L,B} <: RadialNeighborhood{R,L}
    _buffer::B
end
Window(R::Int) = Window{R}()
Window{R}(_buffer=nothing) where R = Window{R,(2R+1)^2}(_buffer)
Window{R,L}(_buffer::B=nothing) where {R,L,B} = Window{R,L,B}(_buffer)
Window(A::AbstractArray) = Window{(size(A, 1) - 1) ÷ 2}()

@inline _setbuffer(::Window{R,L}, buf::B2) where {R,L,B2} = Window{R,L,B2}(buf)

# The central cell is included
@inline function offsets(hood::Window{R}) where R
    D = 2R + 1
    ntuple(i -> (rem(i-1, D)-R, (i-1) ÷ D - R), D^2)
end
bufindices(hood::Window{R}) where R = Base.OneTo((2R+1)^2)

neighbors(hood::Window) = _buffer(hood)

"""
    AbstractKernelNeighborhood <: Neighborhood

Abstract supertype for kernel neighborhoods.

These can wrap any other neighborhood object, and include a kernel of
the same length and positions as the neighborhood.
"""
abstract type AbstractKernelNeighborhood{R,L} <: Neighborhood{R,L} end

neighborhood(hood::AbstractKernelNeighborhood) = hood.neighborhood
neighbors(hood::AbstractKernelNeighborhood) = neighbors(neighborhood(hood))
offsets(hood::AbstractKernelNeighborhood) = offsets(neighborhood(hood))
positions(hood::AbstractKernelNeighborhood, I) = positions(neighborhood(hood), I)
kernel(hood::AbstractKernelNeighborhood) = hood.kernel

kernelproduct(hood::AbstractKernelNeighborhood) =
    kernelproduct(neighborhood(hood), kernel(hood))
function kernelproduct(hood::Neighborhood{<:Any,L}, kernel) where L
    sum = zero(first(hood))
    @simd for i in 1:L
        sum += hood[i] * kernel[i]
    end
    return sum
end
function kernelproduct(hood::Window{<:Any,L}, kernel) where L
    sum = zero(first(hood))
    @simd for i in 1:L
        @inbounds sum += _buffer(hood)[i] * kernel[i]
    end
    return sum
end

"""
    Kernel <: AbstractKernelNeighborhood

    Kernel(neighborhood, kernel)

Wrap any other neighborhood object, and includes a kernel of
the same length and positions as the neighborhood.

`R = 1` gives 3x3 matrices.
"""
struct Kernel{R,L,N,K} <: AbstractKernelNeighborhood{R,L}
    neighborhood::N
    kernel::K
end
Kernel(A::AbstractMatrix) = Kernel(Window(A), A)
function Kernel(hood::N, kernel::K) where {N<:Neighborhood{R,L},K} where {R,L}
    length(hood) == length(kernel) || _kernel_length_error(hood, kernel)
    Kernel{R,L,N,K}(hood, kernel)
end
function Kernel{R,L}(hood::N, kernel::K) where {R,L,N<:Neighborhood{R,L},K}
    Kernel{R,L,N,K}(hood, kernel)
end

@noinline _kernel_length_error(hood, kernel) =
    throw(ArgumentError("Neighborhood length $(length(hood)) does not match kernel length $(length(kernel))"))

@inline function _setbuffer(n::Kernel{R,L,<:Any,K}, buf) where {R,L,K}
    hood = _setbuffer(neighborhood(n), buf)
    return Kernel{R,L,typeof(hood),K}(hood, kernel(n))
end

"""
    AbstractPositionalNeighborhood <: Neighborhood

Neighborhoods are tuples or vectors of custom coordinates tuples
that are specified in relation to the central point of the current cell.
They can be any arbitrary shape or size, but should be listed in column-major
order for performance.
"""
abstract type AbstractPositionalNeighborhood{R,L} <: Neighborhood{R,L} end

const CustomOffset = Tuple{Vararg{Int}}
const CustomOffsets = Union{AbstractArray{<:CustomOffset},Tuple{Vararg{<:CustomOffset}}}

"""
    Positional <: AbstractPositionalNeighborhood

    Positional(coord::Tuple{Vararg{Int}}...)
    Positional(offsets::Tuple{Tuple{Vararg{Int}}})

Neighborhoods that can take arbitrary shapes by specifying each coordinate,
as `Tuple{Int,Int}` of the row/column distance (positive and negative)
from the central point.

The neighborhood radius is calculated from the most distance coordinate.
For simplicity the buffer read from the main grid is a square with sides
`2r + 1` around the central point.
"""
struct Positional{R,L,O<:CustomOffsets,B} <: AbstractPositionalNeighborhood{R,L}
    "A tuple of tuples of Int, containing 2-D coordinates relative to the central point"
    offsets::O
    _buffer::B
end
Positional(args::CustomOffset...) = Positional(args)
Positional(offsets::CustomOffsets, _buffer=nothing) =
    Positional{_absmaxcoord(offsets),length(offsets)}(offsets, _buffer)
Positional{R,L}(offsets::O, _buffer::B=nothing) where {R,L,O<:CustomOffsets,B} =
    Positional{R,L,O,B}(offsets, _buffer)

# Calculate the maximum absolute value in the offsets to use as the radius
_absmaxcoord(offsets::Union{AbstractArray,Tuple}) = maximum(map(x -> maximum(map(abs, x)), offsets))

ConstructionBase.constructorof(::Type{Positional{R,L,C,B}}) where {R,L,C,B} =
    Positional{R,L}

Base.length(hood::Positional) = length(offsets(hood))

offsets(hood::Positional) = hood.offsets
@inline _setbuffer(n::Positional{R,L,O}, buf::B2) where {R,L,O,B2} =
    Positional{R,L,O,B2}(offsets(n), buf)

"""
    LayeredPositional <: AbstractPositional

    LayeredPositional(layers::Positional...)

Sets of [`Positional`](@ref) neighborhoods that can have separate rules for each set.

`neighbors` for `LayeredPositional` returns a tuple of iterators
for each neighborhood layer.
"""
struct LayeredPositional{R,L,La,B} <: AbstractPositionalNeighborhood{R,L}
    "A tuple of custom neighborhoods"
    layers::La
    _buffer::B
end
LayeredPositional(layers::Positional...) = LayeredPositional(layers)
LayeredPositional(layers::Tuple{Vararg{<:Positional}}, _buffer=nothing) =
    LayeredPositional{maximum(map(radius, layers)),map(length, layers)}(layers, _buffer)
LayeredPositional{R,L}(layers, _buffer::B) where {R,L,B} = begin
    # Child layers must have the same _buffer, and the same R and L parameters
    layers = map(l -> Positional{R,L}(offsets(l), _buffer), layers)
    LayeredPositional{R,L,typeof(layers),B}(layers, _buffer)
end

@inline neighbors(hood::LayeredPositional) = map(l -> neighbors(l), hood.layers)
@inline offsets(hood::LayeredPositional) = map(l -> offsets(l), hood.layers)
@inline positions(hood::LayeredPositional, args...) = map(l -> positions(l, args...), hood.layers)
@inline _setbuffer(n::LayeredPositional{R,L}, buf) where {R,L} =
    LayeredPositional{R,L}(n.layers, buf)

@inline Base.sum(hood::LayeredPositional) = map(sum, neighbors(hood))

"""
    VonNeumann(radius=1) -> Positional

A convenience wrapper to build Von-Neumann neighborhoods as
a [`Positional`](@ref) neighborhood.
"""
function VonNeumann(radius=1, _buffer=nothing)
    offsets = Tuple{Int,Int}[]
    rng = -radius:radius
    for j in rng, i in rng
        distance = abs(i) + abs(j)
        if distance <= radius && distance > 0
            push!(offsets, (i, j))
        end
    end
    return Positional(Tuple(offsets), _buffer)
end

# Get the size of a neighborhood dimension from its radius,
# which is always 2r + 1.
@inline hoodsize(hood::Neighborhood{R}) where R = hoodsize(R)
@inline hoodsize(radius::Integer) = 2radius + 1


text:
@doc raw"""
    Elliptope{N,K} <: AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType}

The Elliptope manifold, also known as the set of correlation matrices, consists of all
symmetric positive semidefinite matrices of rank $k$ with unit diagonal, i.e.,

````math
\begin{aligned}
\mathcal E(n,k) =
\bigl\{p ∈ ℝ^{n × n}\ \big|\ &a^\mathrm{T}pa \geq 0 \text{ for all } a ∈ ℝ^{n},\\
&p_{ii} = 1 \text{ for all } i=1,\ldots,n,\\
&\text{and } p = qq^{\mathrm{T}} \text{ for } q \in  ℝ^{n × k} \text{ with } \operatorname{rank}(p) = \operatorname{rank}(q) = k
\bigr\}.
\end{aligned}
````

And this manifold is working solely on the matrices $q$. Note that this $q$ is not unique,
indeed for any orthogonal matrix $A$ we have $(qA)(qA)^{\mathrm{T}} = qq^{\mathrm{T}} = p$,
so the manifold implemented here is the quotient manifold. The unit diagonal translates to
unit norm columns of $q$.

The tangent space at $p$, denoted $T_p\mathcal E(n,k)$, is also represented by matrices
$Y\in ℝ^{n × k}$ and reads as

````math
T_p\mathcal E(n,k) = \bigl\{
X ∈ ℝ^{n × n}\,|\,X = qY^{\mathrm{T}} + Yq^{\mathrm{T}} \text{ with } X_{ii} = 0 \text{ for } i=1,\ldots,n
\bigr\}
````
endowed with the [`Euclidean`](@ref) metric from the embedding, i.e. from the $ℝ^{n × k}$

This manifold was for example
investigated in[^JourneeBachAbsilSepulchre2010].

# Constructor

    Elliptope(n,k)

generates the manifold $\mathcal E(n,k) \subset ℝ^{n × n}$.

[^JourneeBachAbsilSepulchre2010]:
    > Journée, M., Bach, F., Absil, P.-A., and Sepulchre, R.:
    > “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices”,
    > SIAM Journal on Optimization (20)5, pp. 2327–2351, 2010.
    > doi: [10.1137/080731359](https://doi.org/10.1137/080731359),
    > arXiv: [0807.4423](http://arxiv.org/abs/0807.4423).
"""
struct Elliptope{N,K} <: AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType} end

Elliptope(n::Int, k::Int) = Elliptope{n,k}()

@doc raw"""
    check_point(M::Elliptope, q; kwargs...)

checks, whether `q` is a valid reprsentation of a point $p=qq^{\mathrm{T}}$ on the
[`Elliptope`](@ref) `M`, i.e. is a matrix
of size `(N,K)`, such that $p$ is symmetric positive semidefinite and has unit trace.
Since by construction $p$ is symmetric, this is not explicitly checked.
Since $p$ is by construction positive semidefinite, this is not checked.
The tolerances for positive semidefiniteness and unit trace can be set using the `kwargs...`.
"""
function check_point(M::Elliptope{N,K}, q; kwargs...) where {N,K}
    mpv = invoke(check_point, Tuple{supertype(typeof(M)),typeof(q)}, M, q; kwargs...)
    mpv === nothing || return mpv
    row_norms_sq = sum(abs2, q; dims=2)
    if !all(isapprox.(row_norms_sq, 1.0; kwargs...))
        return DomainError(
            row_norms_sq,
            "The point $(q) does not represent a point p=qq^T on $(M) diagonal is not only ones.",
        )
    end
    return nothing
end

@doc raw"""
    check_vector(M::Elliptope, q, Y; kwargs... )

Check whether $X = qY^{\mathrm{T}} + Yq^{\mathrm{T}}$ is a tangent vector to
$p=qq^{\mathrm{T}}$ on the [`Elliptope`](@ref) `M`,
i.e. `Y` has to be of same dimension as `q` and a $X$ has to be a symmetric matrix with
zero diagonal.

The tolerance for the base point check and zero diagonal can be set using the `kwargs...`.
Note that symmetric of $X$ holds by construction an is not explicitly checked.
"""
function check_vector(M::Elliptope{N,K}, q, Y; kwargs...) where {N,K}
    mpv = invoke(
        check_vector,
        Tuple{supertype(typeof(M)),typeof(q),typeof(Y)},
        M,
        q,
        Y;
        kwargs...,
    )
    mpv === nothing || return mpv
    X = q * Y' + Y * q'
    n = diag(X)
    if !all(isapprox.(n, 0.0; kwargs...))
        return DomainError(
            n,
            "The vector $(X) is not a tangent to a point on $(M) (represented py $(q) and $(Y), since its diagonal is nonzero.",
        )
    end
    return nothing
end

function decorated_manifold(M::Elliptope)
    return Euclidean(representation_size(M)...; field=ℝ)
end

@doc raw"""
    manifold_dimension(M::Elliptope)

returns the dimension of
[`Elliptope`](@ref) `M`$=\mathcal E(n,k), n,k ∈ ℕ$, i.e.
````math
\dim \mathcal E(n,k) = n(k-1) - \frac{k(k-1)}{2}.
````
"""
@generated function manifold_dimension(::Elliptope{N,K}) where {N,K}
    return N * (K - 1) - div(K * (K - 1), 2)
end

"""
    project(M::Elliptope, q)

project `q` onto the manifold [`Elliptope`](@ref) `M`, by normalizing the rows of `q`.
"""
project(::Elliptope, ::Any)

project!(::Elliptope, r, q) = copyto!(r, q ./ (sqrt.(sum(abs2, q, dims=2))))

"""
    project(M::Elliptope, q, Y)

Project `Y` onto the tangent space at `q`, i.e. row-wise onto the oblique manifold.
"""
project(::Elliptope, ::Any...)

function project!(::Elliptope, Z, q, Y)
    Y2 = (Y' - q' .* sum(q' .* Y', dims=1))'
    Z .= Y2 - q * lyap(q' * q, q' * Y2 - Y2' * q)
    return Z
end

@doc raw"""
    retract(M::Elliptope, q, Y, ::ProjectionRetraction)

compute a projection based retraction by projecting $q+Y$ back onto the manifold.
"""
retract(::Elliptope, ::Any, ::Any, ::ProjectionRetraction)

retract!(M::Elliptope, r, q, Y, ::ProjectionRetraction) = project!(M, r, q + Y)

@doc raw"""
    representation_size(M::Elliptope)

Return the size of an array representing an element on the
[`Elliptope`](@ref) manifold `M`, i.e. $n × k$, the size of such factor of $p=qq^{\mathrm{T}}$
on $\mathcal M = \mathcal E(n,k)$.
"""
@generated representation_size(::Elliptope{N,K}) where {N,K} = (N, K)

function Base.show(io::IO, ::Elliptope{N,K}) where {N,K}
    return print(io, "Elliptope($(N), $(K))")
end

"""
    vector_transport_to(M::Elliptope, p, X, q)

transport the tangent vector `X` at `p` to `q` by projecting it onto the tangent space
at `q`.
"""
vector_transport_to(::Elliptope, ::Any, ::Any, ::Any, ::ProjectionTransport)

function vector_transport_to!(M::Elliptope, Y, p, X, q, ::ProjectionTransport)
    project!(M, Y, q, X)
    return Y
end

@doc raw"""
    zero_vector(M::Elliptope,p)

returns the zero tangent vector in the tangent space of the symmetric positive
definite matrix `p` on the [`Elliptope`](@ref) manifold `M`.
"""
zero_vector(::Elliptope, ::Any...)

zero_vector!(::Elliptope{N,K}, v, ::Any) where {N,K} = fill!(v, 0)


text:
function TimeComponents(arg0::jdouble)
    return TimeComponents((jdouble,), arg0)
end

function TimeComponents(arg0::jint, arg1::jdouble)
    return TimeComponents((jint, jdouble), arg0, arg1)
end

function TimeComponents(arg0::jint, arg1::jint, arg2::jdouble)
    return TimeComponents((jint, jint, jdouble), arg0, arg1, arg2)
end

function TimeComponents(arg0::jint, arg1::jint, arg2::jdouble, arg3::jint)
    return TimeComponents((jint, jint, jdouble, jint), arg0, arg1, arg2, arg3)
end

function compare_to(obj::TimeComponents, arg0::TimeComponents)
    return jcall(obj, "compareTo", jint, (TimeComponents,), arg0)
end

function equals(obj::TimeComponents, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function from_seconds(::Type{TimeComponents}, arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jint)
    return jcall(TimeComponents, "fromSeconds", TimeComponents, (jint, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function get_hour(obj::TimeComponents)
    return jcall(obj, "getHour", jint, ())
end

function get_minute(obj::TimeComponents)
    return jcall(obj, "getMinute", jint, ())
end

function get_minutes_from_utc(obj::TimeComponents)
    return jcall(obj, "getMinutesFromUTC", jint, ())
end

function get_second(obj::TimeComponents)
    return jcall(obj, "getSecond", jdouble, ())
end

function get_seconds_in_local_day(obj::TimeComponents)
    return jcall(obj, "getSecondsInLocalDay", jdouble, ())
end

function get_seconds_in_utc_day(obj::TimeComponents)
    return jcall(obj, "getSecondsInUTCDay", jdouble, ())
end

function hash_code(obj::TimeComponents)
    return jcall(obj, "hashCode", jint, ())
end

function parse_time(::Type{TimeComponents}, arg0::JString)
    return jcall(TimeComponents, "parseTime", TimeComponents, (JString,), arg0)
end

function to_string(obj::TimeComponents)
    return jcall(obj, "toString", JString, ())
end


text:
using Test
using JSExpr

@testset "objects" begin
    @testset "braces syntax" begin
        @test string(@js(
            {
                "foo" => "bar",
                "spam" => 1 + 2,
            }
        )) == """{"foo": "bar", "spam": 1 + 2}"""

        @test string(@js(
            {
                :foo => "bar",
                :spam => 1 + 2,
            }
        )) == """{"foo": "bar", "spam": 1 + 2}"""

        @test string(@js(
            {
                foo = "bar",
                "spam" => 1 + 2,
            }
        )) == """{"foo": "bar", "spam": 1 + 2}"""

        @test string(@js(
            { foo }
        )) == """{"foo": foo}"""

        @test string(@js(
            {
                a,
                b = "foo",
                :c => "bar",
            }
        )) == """{"a": a, "b": "foo", "c": "bar"}"""

        @test_throws ErrorException JSExpr.crawl(:({ foo => "bar" }))
        @test_throws ErrorException JSExpr.crawl(:({ "foo" }))
        @test_throws ErrorException JSExpr.crawl(:({ :foo=foo }))
    end

    @testset "namedtuple syntax" begin
        @test string(@js(
            (
                foo="foo",
                bar="bar",
            )
        )) == """{"foo": "foo", "bar": "bar"}"""

        @test_throws ErrorException JSExpr.crawl(:(
            ("foo"=3, foo="bar")
        ))
        @test_throws ErrorException JSExpr.crawl(:(
            (foo="foo", bar)
        ))
    end
end


text:
@testset "AlexNet" begin
  model = alexnet()
  @test size(model(rand(256, 256, 3, 50))) == (1000, 50)
end


text:
__precompile__()

module ImagineHardware

using Unitful, UnitAliases
import Unitful: V, μm, m, A, mA, N, s, Hz, g, μF
import Base: position, clamp

include("amplifiers.jl")
include("registered_amplifiers.jl")
include("pids.jl")
include("registered_pids.jl")
include("piezos.jl")
include("registered_positioners.jl")
include("piezosystems.jl")
include("registered_piezosystems.jl")

export  AbstractAmplifier,
        SimpleAmplifier,
        retrieve_amplifier, register_amplifier,
        amplify, gain, offset, in_min, in_max, out_min, out_max, current_max, name,

        PidController,
        kp, ki, kd,
        set_kp!, set_ki!, set_kd!,
        samprate,
        retrieve_pid, register_pid,

        AbstractPositioner,
        Piezo,
        retrieve_positioner, register_positioner,
	pmin, pmax, pspan, vmin_in, vmax_in, vspan_in, vmin_out, vmax_out, vspan_out, pos2mon, mon2pos,
        closed_loop_pad,
	max_displacement,
	capacitance,
	motion_axis,
	stiffness, stiffness_all,
        clamp,
	effective_mass,
	resonance, resonance_all,
	max_force,
        loaded_resonance,
        closed2open, open2closed,
        pos2mon, mon2pos,
        pos2mod, mod2pos,
        pos2displacement, displacement2pos,
        input2displacement, displacement2input,
        mon2displacement, displacement2mon,
        mod2displacement, displacement2mod,
        input2pos, pos2input,

        PiezoSystem,
        amplifier, piezo, pid_controller,
        retrieve_piezosystem, register_piezosystem

end


text:
using Revise, PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
const PDMP = PiecewiseDeterministicMarkovProcesses

function F_fd!(ẋ, xc, xd, parms, t)
	# vector field used for the continuous variable
	if mod(xd[1], 2) == 0
		ẋ[1] = 1 + xd[1]
	else
		ẋ[1] = -xc[1]
	end
	nothing
end

rate_tcp(x) = 1/x

function R_fd!(rate, xc, xd, parms, t, issum::Bool)
	rate[1] = 1.0 + rate_tcp(xd[1]) * xc[1]
	if issum == false
		return 0., 0.
	else
		return sum(rate), 0.
	end
end

Dummy! = PDMP.Delta_dummy

xc0 = [ 1.0 ]
xd0 = [ 1  ]

nu_fd = [[1 0];[0 -1]]
parms = [0.0]

# works:
Random.seed!(12)
	problem = PDMP.PDMPProblem(F_fd!, R_fd!, nu_fd, xc0, xd0, parms, (0.0, 10000.))
	res =  @time PDMP.solve(problem, CHV(CVODE_Adams()); save_positions = (false, false), n_jumps = 3000)
	# res =  @time PDMP.pdmp!(xc0, xd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0; algo = :chv, ode = CVODE_Adams()) #.967ms 4.38k allocations

Random.seed!(12)
	problem = PDMP.PDMPProblem(F_fd!, R_fd!, nu_fd, xc0, xd0, parms, (0.0, 10000.))
	res =  @time PDMP.solve(problem, CHV(Tsit5()); save_positions = (false, false), n_jumps = 3000)
	# res =  @time PDMP.pdmp!(xc0, xd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0; algo = :chv, n_jumps = 3000,   ode = Tsit5(), save_positions=(false,false)) #1.037ms 466 allocations

# Random.seed!(12)
# 	res =  @time  PDMP.chv_diffeq!(xc0, xd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0,false; n_jumps = 3000, ode = Tsit5() ,save_positions = (false, false), rate = zeros(2), xc0_extended = zeros(2))

Random.seed!(12)
	problem = PDMP.PDMPProblem(F_fd!, R_fd!, nu_fd, xc0, xd0, parms, (0.0, 10000.))
	res =  @time PDMP.solve(problem, CHV(AutoTsit5(Rosenbrock23(autodiff=true))); save_positions = (false, false), n_jumps = 3000)
	# res =  @time PDMP.pdmp!(xc0, xd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0; algo = :chv, n_jumps = 3000,   ode = AutoTsit5(Rosenbrock23(autodiff=true)), save_positions=(false,false)) #9ms

# fail because of autodiff
Random.seed!(12)
	problem = PDMP.PDMPProblem(F_fd!, R_fd!, nu_fd, xc0, xd0, parms, (0.0, 10000.))
	res =  @time PDMP.solve(problem, CHV(TRBDF2(autodiff=true)); save_positions = (false, false), n_jumps = 3000)

using StaticArrays
sxc0 = @MVector [ 1.0 ]
sxd0 = @MVector [1]
ratevec = similar(sxc0, Size(2))
sxc0_e = similar(sxc0, Size(2))

problem = PDMP.PDMPProblem(F_fd!, R_fd!, nu_fd, sxc0, sxd0, parms, (0.0, 10000.))
res =  @time PDMP.solve(problem, CHV(Tsit5()); save_positions = (false, false), n_jumps = 3000)
# ress =  @time  PDMP.chv_diffeq!(sxc0, sxd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0,false; n_jumps = 3000, ode = Tsit5(),save_positions = (false, false), rate = ratevec, xc0_extended = sxc0_e)

length(sxc0)


text:
struct ComponentDescription
    name::Symbol
    componenttype::Type{<:Component}
    attributes::DataFrame
end

const componentdescriptions = ComponentDescription[]

function componentdescription(name::Symbol)
    j = findfirst(cd -> cd.name == name, componentdescriptions)
    !isnothing(j) ? componentdescriptions[j] : error("Component $componentname has not been registered")
end
function componentdescription(::Type{T}) where T <: Component
    j = findfirst(cd -> T <: cd.componenttype, componentdescriptions)
    !isnothing(j) ? componentdescriptions[j] : error("Component type $T has not been registered")
end

function getdefault(c::Component, attr::Symbol)
    attrs = attributes(typeof(c))
    i = findfirst(isequal(attr), attrs.attribute)
    isnothing(i) && error("Attribute $attr for component $c is not known")
    ret = attrs.default[i]
    ismissing(ret) && error("Attribute $attr for component $c does not have a default")
    ret
end

resolve(::Type{Component}, name::Symbol) = componentdescription(name).componenttype
attributes(name) = componentdescription(name).attributes

addcomponent(cd::ComponentDescription) where T<:Component = push!(componentdescriptions, cd)
addcomponent(::Type{T}) where T<:Component = addcomponent(T, naming(T))
addcomponent(::Type{T}, name::Symbol) where T<:Component = addcomponent(ComponentDescription(name, T, DataFrame()))

function addcomponent(::Type{T}, name::Symbol, axes, filename) where T<:Component
    axes = (first(axes)=>name, Base.tail(axes)...)
    df = DataFrame(CSV.read(filename, truestrings=["t"], falsestrings=["f"]))
    df.attribute = Symbol.(df.attribute)
    df.default = map(r->astype(r.dtype, r.default), eachrow(df))
    rename_dimensions(x) = tuple(recode(Symbol.(split(x, ',')), axes...)...)
    df.dimensions = map(x->ismissing(x) ? () : rename_dimensions(x), df.dimensions)

    addcomponent(ComponentDescription(name, T, df))
end

function naming(::Type{T}) where T <: Component
    j = findfirst(cd -> cd.componenttype == T, componentdescriptions)
    !isnothing(j) && return componentdescriptions[j].name

    S = T
    while S isa UnionAll
        S = S.body
    end

    s = lowercase(string(S.name.name))
    Symbol(s, in(s[end], ('s', 'h')) ? "es" : "s")
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
using Thermodynamics, Documenter
using DocumenterCitations

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

pages = Any[
    "Home" => "index.md",
    "Installation" => "Installation.md",
    "API" => "API.md",
    "How-to-guide" => "HowToGuide.md",
    "Tested profiles" => "TestedProfiles.md",
    "Temperature profiles" => "TemperatureProfiles.md",
    "Developer docs" => "DevDocs.md",
    "Thermodynamics overview" => "Formulation.md",
    "References" => "References.md",
]

mathengine = MathJax(
    Dict(
        :TeX => Dict(
            :equationNumbers => Dict(:autoNumber => "AMS"),
            :Macros => Dict(),
        ),
    ),
)

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
)

makedocs(
    bib,
    sitename = "Thermodynamics.jl",
    strict = true,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [Thermodynamics],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/Thermodynamics.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
#
using Test
using SpectralElements

dir = "../examples"
files = [
         "p2d.jl",
        ]
for file in files
    @time @testset "$file" begin include("$dir/$file") end # @safetestset ?
end
#


text:
# Our colleague's project depends on a useless package
# (https://github.com/crstnbr/Useless.jl)
using Useless

life_universe_everything(question) = 42

println(life_universe_everything("What is the meaning of life?"))

# Step 1: instantiate the environment
# Step 2: run the code: julia --project=. code.jl


text:
include("Include.jl")

# path to file -
path_model_file = "Coagulation.net"

# load the file -
model_buffer = read_model_file(path_model_file)

# build the default dictionary -
dd = build_default_model_dictionary(model_buffer)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module CodedComputing

using LinearAlgebra, SparseArrays, HDF5, DataFrames
using Distributions
using DataStructures

include("Linalg.jl")
include("concurrentcircularbuffer.jl")
include("mul.jl")
include("partitioning.jl")
include("Datasets.jl")
include("profiling.jl")
include("eventdriven.jl")
include("loadbalancing.jl")
include("treegradient.jl")

end


text:
using MomentClosure
using MomentClosure: expected_coeff
using ModelingToolkit: value, simplify
using SymbolicUtils: expand
using Test

@parameters b
b = b.val
@test expected_coeff(2, 3) == 8
@test isequal(expected_coeff(b, 3), b*(1 + (6b) + (6(b^2))))
expr = ((b*((1 + b)^-1)) + ((b^6)*((1 + b)^-6)) + (57(b^2)*((1 + b)^-2)) + (302(b^3)*((1 + b)^-3)) + (302(b^4)*((1 + b)^-4)) + (57(b^5)*((1 + b)^-5)))*((1 + b)^6)
@test isequal(expected_coeff(b, 6), expr)

@parameters t, c₁, c₂, c₃, c₄, Ω
@variables X(t), Y(t)

S_mat = [ 1 -1  1 -1;
         -1  1  0  0]
a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X]
rn = ReactionSystemMod(t, [X, Y], [c₁, c₂, c₃, c₄, Ω], a, S_mat)

sys = generate_central_moment_eqs(rn, 2, 4)
expr1 = sys.odes.eqs[2].rhs
@test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs)
μ = sys.μ
M = sys.M
# simplify(value()) is needed due to diffs between ModelingToolkit Num type and SymbolicUtils types
expr2 = c₂*μ[1,0] + c₁*M[1,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*Ω^-2 - c₁*M[2,1]*Ω^-2 -
        2*c₁*M[1,1]*μ[1,0]*Ω^-2 - c₁*M[2,0]*μ[0,1]*Ω^-2 - c₁*μ[0,1]*Ω^-2*μ[1,0]^2
expr2 = simplify(value.(expr2))
@test isequal(expand(expr1), expr2)

sys = generate_central_moment_eqs(rn, 2)
expr1 = sys.odes.eqs[2].rhs
@test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs)
@test isequal(expand(expr1), expr2)

sys = generate_raw_moment_eqs(rn, 2)
μ = sys.μ
expr1 = sys.odes.eqs[4].rhs
expr2 = c₂*μ[2,0] + c₁*μ[1,1]/Ω^2 - c₁*μ[3,1]/Ω^2 -c₂*(μ[1,0] + μ[1,1]) +
        c₃*μ[0,1]*Ω - c₄*μ[1,1] - c₁*μ[1,2]/Ω^2 + c₁*μ[2,2]/Ω^2
expr2 = simplify(value.(expr2))
@test isequal(expand(expr1), expand(expr2))


text:
# Overarching parent type
abstract type GQLItem end

##########
# Values #
##########

struct Loc <: GQLItem
    line::Int
    column::Int
end
Base.show(io::IO, ::MIME"text/plain", loc::Loc) = print(io, "Line $(loc.line) Column $(loc.column)")
Base.show(io::IO, loc::Loc) = print(io, "Line $(loc.line) Column $(loc.column)") # for string interpolation

struct Variable <: GQLItem
    name::String
    loc::Loc
end
Variable(name) = Variable(name, Loc(0,0))
struct Alias # TODO: remove?
    alias::String
    name::String
    loc::Loc
end
Alias(alias, name) = Alias(alias, name, Loc(0,0))

struct Enum <: GQLItem
    name::String
    loc::Loc
end
Enum(name) = Enum(name, Loc(0,0))

struct ObjectField <: GQLItem
    name::String
    value
    loc::Loc
end
ObjectField(name, value) = ObjectField(name, value, Loc(0,0))

struct InputObject <: GQLItem
    object_fields::Vector{ObjectField}
    loc::Loc
end
InputObject(object_fields) = InputObject(object_fields, Loc(0,0))

#################
# SelectionSets #
#################

struct Argument <: GQLItem
    name::String
    value
    loc::Loc
end
Argument(name, value) = Argument(name, value, Loc(0,0))

struct Directive <: GQLItem
    name::String
    arguments::Union{Nothing, Vector{Argument}}
    loc::Loc
end
Directive(name, arguments) = Directive(name, arguments, Loc(0,0))

abstract type Selection  <: GQLItem end

struct Field{T} <: Selection
    alias::Union{Nothing, String}
    name::String
    arguments::Union{Nothing, Vector{Argument}}
    directives
    selection_set::T # avoid circular deps
    loc::Loc
end
Field(alias, name, args, dirs, ss) = Field(alias, name, args, dirs, ss, Loc(0,0))

struct Fragment <: Selection end

struct SelectionSet <: GQLItem
    selections::Vector{Selection}
    loc::Loc
end
SelectionSet(selections) = SelectionSet(selections, Loc(0,0))

struct FragmentSpread <: Selection
    name::String
    directives::Union{Vector{Directive}, Nothing}
    loc::Loc
end
FragmentSpread(name, directives) = FragmentSpread(name, directives, Loc(0,0))

struct InlineFragment <: Selection
    named_type::Union{String, Nothing}
    directives::Union{Vector{Directive}, Nothing}
    selection_set::SelectionSet
    loc::Loc
end
InlineFragment(name, dirs, ss) = InlineFragment(name, dirs, ss, Loc(0,0))

function print_selection(io, field::Field)
    !isnothing(field.alias) && print(io, field.alias, ": ")
    print(io, field.name)
    !isnothing(field.arguments) && print(io, "(args...)")
    !isnothing(field.directives) && print(io, "@...")
    !isnothing(field.selection_set) && print(io, "{selection_set...}")
end

function print_selection_set(io, selection_set::SelectionSet)
    print(io, "{\n")
    for selection in selection_set.selections
        print(io, "    ")
        print_selection(io, selection)
        print(io, "\n")
    end
    print(io, "}")
end
function Base.show(io::IO, ::MIME"text/plain", selection_set::SelectionSet)
    print_selection_set(io, selection_set)
end

struct VariableDefinition <: GQLItem
    name::String
    type::String
    value
    directives::Union{Nothing, Vector{Directive}}
    loc::Loc
end
VariableDefinition(name, type, value, dirs) = VariableDefinition(name, type, value, dirs, Loc(0,0))

###############
# Definitions #
###############

abstract type Definition  <: GQLItem end

struct FragmentDefinition <: Definition
    name::String
    named_type::String
    directives::Union{Nothing, Vector{Directive}}
    selection_set::SelectionSet
    loc::Loc
end
FragmentDefinition(name, named_type, dirs, ss) = FragmentDefinition(name, named_type, dirs, ss, Loc(0,0))

struct Operation <: Definition
    operation_type::String
    name::Union{Nothing, String}
    variable_definitions::Union{Nothing, Vector{VariableDefinition}}
    directives::Union{Nothing, Vector{Directive}}
    selection_set::SelectionSet
    loc::Loc
end
Operation(op_type, named_type, var_defs, dirs, ss) = Operation(op_type, named_type, var_defs, dirs, ss, Loc(0,0))

function Base.show(io::IO, ::MIME"text/plain", op::Operation)
    println(io, "Operation")
    print(io, "\n")
    print(io, "$(op.operation_type)")
    !isnothing(op.name) && print(io, " ", op.name)
    !isnothing(op.variable_definitions) && print(io, "(vars...)")
    !isnothing(op.directives) && print(io, "@...")
    print_selection_set(io, op.selection_set)
end

function Base.:(==)(n1::T, n2:: T) where T <: GQLItem
    for name in fieldnames(T)
        if name != :loc
            !isequal(getproperty(n1, name), getproperty(n2, name)) && return false
        end
    end
    return true
end

function Base.hash(n::T, h::UInt) where T <: GQLItem
    for name in fieldnames(T)
        if name != :loc
            h = hash(getproperty(n, name), h)
        end
    end
    return h
end

############
# Document #
############

struct Document <: GQLItem
    definitions::Vector{Definition}
end


text:
import Neurthino: PREM

@test PREM(0) ≈ 13.0885
@test PREM(1000) ≈ 12.87075725
@test PREM(1234) ≈ 12.15988934
@test PREM(6370) ≈ 1.02


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

using Dates: DateTime, Second

struct AccessToken
    token::String
    scope::String
    expires_in::Int  # seconds
    created_on::DateTime
end

function isexpired(access_token::AccessToken)::Bool
    expires_on = access_token.created_on + Second(access_token.expires_in)
    return expires_on - Second(5) < now() # anticipate token expiration by 5 seconds
end

abstract type Credentials end

mutable struct ClientCredentials <: Credentials
    client_id::String
    client_secret::String
    client_credentials_url::Union{String,Nothing}
    access_token::Union{AccessToken,Nothing}
    ClientCredentials(client_id, client_secret, client_credentials_url = nothing) =
        new(client_id, client_secret, client_credentials_url, nothing)
end

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
using Potree
using FileManager
using Visualization
using BenchmarkTools
using Profile

# potree source folder
potree = "C:/Users/Alessio/Documents/potreeDirectory/pointclouds/CAVA" # replace this path with local potree directory

println("==============================================================================")
println("@btime for FileManager: potree2triee, get_all_values, las2pointcloud")
FMtrie = @btime FileManager.potree2trie($potree)
FMall_files = @btime FileManager.get_all_values($FMtrie)
FMPC = @btime FileManager.las2pointcloud($FMall_files...)
println("==============================================================================")
println("==============================================================================")
println("@btime for Potree: potree2triee, get_all_values, las2pointcloud")
Ptrie = @btime Potree.potree2trie($potree)
PMall_files = @btime Potree.get_all_values($Ptrie)
PPC = @btime Potree.las2pointcloud($PMall_files...)
println("==============================================================================")

# Section used to check if there is any difference
# during the output of the 2 packages

# println("==============================================================================")
# println("==============================================================================")
# println("FMall_files == PMall_files (Next line should be true)")
# println(FMall_files == PMall_files)
# println("==============================================================================")
# println("size(FMPC.coordinates) == size(PPC.coordinates) (Next line should be true)")
# println(size(FMPC.coordinates) == size(PPC.coordinates))
# println("==============================================================================")
# println("==============================================================================")

# Section used to check if there is any difference
# during the Visualization of the results

# point cloud
# GL.VIEW(
#     [
#     Visualization.points_color_from_rgb(PPC.coordinates,PPC.rgbs)
#     ]
# )

# # point cloud
# GL.VIEW(
#     [
#     Visualization.points_color_from_rgb(FMPC.coordinates,FMPC.rgbs)
#     ]
# )


text:
"""
Data structure useful for tracking the state of one or more instances of the
Paxos algorithm.
"""
module Ledgers

export Ledger,
  LedgerEntryState,
  LedgerEntry,
  entryOpen,
  entryPromised,
  entryAccepted,
  entryApplied,
  nextInstance,
  addEntry,
  nextBallotNumber!,
  apply!,
  votes!,
  prepare!,
  promise!,
  accept!

using ..Ballots
using ..Nodes

"""
`open` - The entry is new, and has not been populated by any ballot activity

`promised`- Once promised, only requests with a higher ballot can replace
the request in the entry

`accepted` - An accepted entry will eventually by accepted by all members of the cluster

`applied` - The underlying command in the entry's request has been applied
to application state
"""
@enum LedgerEntryState begin
  entryOpen
  entryPromised
  entryAccepted
  entryApplied
end

"""
A `LedgerEntry` captures the state of a `Command` to apply to a `Ledger`s external
model or state machine. Conveniently, a ledger entry has a structure suitable
for use as a "ledger" for 1 instance of the Paxos algorithm: as ballots
for that instance progress the entry contains the state needed for a
leader, follower, or learner to move towards an outcome or final ballot
containing the chosen command for application.
"""
mutable struct LedgerEntry
  state::LedgerEntryState
  ballot::Ballot
end

LedgerEntry(ballot::Ballot) = LedgerEntry(entryOpen, ballot)

# LedgerEntry(request::Request) = LedgerEntry(entryOpen, SequenceNumber(1), request)

# LedgerEntry(sequenceNumber::SequenceNumber) = LedgerEntry(entryOpen, sequenceNumber, nothing)

"""
A ledger is a record of `Command`s to apply to an external data structure or
state machine. A ledger is structured as a sequence of `LedgerEntry` objects.
"""
mutable struct Ledger
  """
  The entries in the ledger, accessed by their index. The reason for using
  a `Dict` instead of an array (with a base offset) is to allow for a potentially
  sparse list of entries.
  """
  entries::Dict{InstanceID,LedgerEntry}

  """
  Index of earliest entry in ledger
  """
  earliestIndex::Union{InstanceID,Nothing}
  """
  Index of latest entry in ledger
  """
  latestIndex::Union{InstanceID,Nothing}
  """
  Index of latest applied entry, or nothing if none applied yet
  """
  latestApplied::Union{Integer,Nothing}

  """
  Lock to protect concurrent access to the ledger and its entries
  """
  lock::ReentrantLock
end

Ledger() = Ledger(Dict(), InstanceID(0), nothing, nothing, ReentrantLock())

function withLedger(fn::Function, ledger::Ledger)
  try
    lock(ledger.lock)
    fn(ledger)
  finally
    unlock(ledger.lock)
  end
end

function Base.isempty(ledger::Ledger)
  withLedger(ledger) do ledger
    isempty(ledger.entries)
  end
end

"""
Return true if the ledger is not empty and has unapplied entries, false
otherwise
"""
function Base.isready(ledger::Ledger)
  withLedger(ledger) do ledger
    return if isempty(ledger)
      false
    else
      if ledger.latestApplied == nothing
        false
      else
        (ledger.latestApplied < ledger.latestIndex) && ledger.latestApplied.state == entryAccepted
      end
    end
  end
end

function Base.length(ledger::Ledger)
  length(ledger.entries)
end

"""
Return the next unused instance (actually, the next unused index) in the ledger
"""
function nextInstance(ledger::Ledger)
  withLedger(ledger) do ledger
    (ledger.latestIndex === nothing) ? ledger.earliestIndex : InstanceID(ledger.latestIndex.value + 1)
  end
end

"""
Create a ballot number for an instance
"""
function nextBallotNumber!(ledger::Ledger, leaderID::NodeID, instanceID = nextInstance(ledger))
  withLedger(ledger) do ledger
    entry = ledger.entries[instanceID]
    entry.ballot.sequenceNumber += 1
    entry.ballotNumber
  end
end

function addEntry(ledger::Ledger, leaderID::NodeID, request::Request)
  withLedger(ledger) do ledger
    instanceID = nextInstance(ledger)
    ballotNumber = BallotNumber(instanceID, SequenceNumber(1),leaderID)
    ballot = Ballot(ballotNumber, request)
    entry = LedgerEntry(entryOpen, ballot)
    ledger.entries[instanceID] = entry
    ballot
  end
end

function apply!(fn::Function, ledger::Ledger, state)
  withLedger(ledger) do ledger
    while isready(ledger)
      nextIndex = (ledger.latestApplied == nothing) ? 0 : (ledger.latestApplied + 1)
      entry = ledger.entries[nextIndex]
      # note we are holding a lock while calling into arbitrary code
      # ....may be worth a rethink
      fn(state, entry.ballot.request.command)
      ledger.latestApplied = nextIndex
    end
  end
end

"""
Record choices for each instance in `ballots`
"""
function prepare!(ledger::Ledger, ballots::Dict{InstanceID,Ballot})::Dict{InstanceID,Ballot}
  preparedBallots = Dict{InstanceID,Ballot}()
  withLedger(ledger) do ledger
    for instanceID in ballots
      entry = ledger.entries[instanceID]
      ballot = ballots[instanceID]
      sequenceNumber = ballot.number.sequenceNumber
      if entry.state == entryOpen && entry.sequenceNumber < sequenceNumber
          entry.sequenceNumber = sequenceNumber
          preparedBallots[instanceID] = ballot
      end
    end
  end
  preparedBallots
end

function promise!(ledger::Ledger, promises::Vector{BallotNumber})
  withLedger(ledger) do ledger
    for promise in promises
      instanceID = promise.instanceID
      sequenceNumber = promise.sequenceNumber
      entry = ledger.entries[instanceID]
      if entry.state == entryOpen && entry.sequenceNumber < sequenceNumber
        entry.state = entryPromised
        entry.sequenceNumber = sequenceNumber
      end
    end
  end
end

"""
Compute a ballot number corresponding to the highest sequence number seen
for each instance. This is useful for the "propose" phase of the protocol.
"""
function votes!(ledger::Ledger, ballotNumbers::Vector{BallotNumber})
  withLedger(ledger) do ledger
    map(ballotNumbers) do ballotNumber
      instanceID = ballotNumber.instanceID
      entry = get!(ledger.entries, instanceID, LedgerEntry(ballotNumber.sequenceNumber))
      BallotNumber(instanceID, max(entry.sequenceNumber, ballotNumber.sequenceNumber))
    end
  end
end

function promise!(ledger::Ledger, approved::Dict{InstanceID,BallotNumber})::Dict{InstanceID,BallotNumber}
  promises = Dict{InstanceID,BallotNumber}()
  withLedger(ledger) do ledger
    for instanceID in approved
      ballotNumber = approved[instanceID]
      entry = ledger.entries[instanceID]
      if entry.state == entryOpen && entry.sequenceNumber == ballotNumber.sequenceNumber
        # TODO think about configuration changes
        promises[instanceID] = ballotNumber
      end
    end
  end
  promises
end

function accept!(ledger::Ledger, ballotNumbers::Vector{BallotNumber})
  acceptances = Vector{BallotNumber}()
  withLedger(ledger) do ledger
    accepted = map(ballotNumbers) do ballotNumber
      instanceID = ballotNumber.instanceID
      entry = ledger.entries[instanceID]
      # TOOO hmmm, should we check leader ID who generated the ballot too?
      if entry.state == entryOpen && entry.sequenceNumber == ballotNumber.sequenceNumber
        entry.state = entryAccepted
        push!(acceptances, ballotNumber)
      end
    end
  end
  acceptances
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 18782940-4b5a-11eb-1a4c-d5d234c2a7bd
begin
using MLJ
#import RDatasets: dataset
using Gadfly
end

# ╔═╡ 254375d0-4b5a-11eb-2717-e3ba997469a1
begin
	function init()
	X = randn(40, 2)
	y = vcat(-ones(20), ones(20))
	return X,y
	end
	X,y=init()
end

# ╔═╡ e6905d60-4b5b-11eb-1ea3-d708f75a0657
begin
ym1 = y .== -1
ym2 = .!ym1
	plot(
layer(x=X[ym1, 1], y=X[ym1, 2], ),
layer(x=X[ym2, 1], y=X[ym2, 2], Theme(default_color=colorant"red")),
 #Scale.shape_discrete(levels=["+ve","-ve"]),
 Guide.manual_color_key("", ["+ve", "-ve"], ["deepskyblue","red"]),
	point_shapes=[Shape.circle, Shape.star1])
	end

# ╔═╡ 6fec4c40-4b61-11eb-0f31-698d897b7b3c
begin
	tX = MLJ.table(X)
	ty = categorical(y);
end

# ╔═╡ a0100f10-4b61-11eb-2944-2d6da12d683a
begin
	@load SVC pkg=LIBSVM

	svc_mdl = SVC()
	svc = machine(svc_mdl, tX, ty)

	fit!(svc);
end

# ╔═╡ a77e3830-4b61-11eb-34f1-4dde97f2c85f
begin
	ypred = predict(svc, tX)
	misclassification_rate(ypred, ty)
end

# ╔═╡ Cell order:
# ╠═18782940-4b5a-11eb-1a4c-d5d234c2a7bd
# ╠═254375d0-4b5a-11eb-2717-e3ba997469a1
# ╠═e6905d60-4b5b-11eb-1ea3-d708f75a0657
# ╠═6fec4c40-4b61-11eb-0f31-698d897b7b3c
# ╠═a0100f10-4b61-11eb-2944-2d6da12d683a
# ╠═a77e3830-4b61-11eb-34f1-4dde97f2c85f

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
"""
Decs package

Copyright 2019 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
Originally based/inspired by Decimals.jl (by jack@tinybike.net (Jack Peterson), 7/3/2014)
"""
module Decs

macro dec_str(s, flags...) parse(Dec, s) end

# Note: I don't implement the number() function that Decimals had
# (I don't believe it should be part of the API)

import Base: ==, +, -, *, /, <, inv, round

export Dec, dec, @dec_str

const DIGITS = 20

const DFP_MARKER = BigInt(0)

struct Dec <: AbstractFloat
    sgn::Bool   # true = negative
    pow::Int32  # power of 10
    val::BigInt # value
end

Dec(x::Dec) = x
Dec(x::Real) = parse(Dec, string(x))
Dec(x::Integer, scl=0) = Dec(x < 0, scl, abs(x))

# Until I change to use DecFP for speed and to be able to support +-Inf and NaN:
Base.isfinite(x::Dec) = true
Base.isnan(x::Dec) = false

const dectab = Vector{BigInt}(undef, 128)
const rndtab = Vector{BigInt}(undef, 128)

function __init__()
    for p = 1:length(dectab); v = BigInt(10)^p; dectab[p] = v; rndtab[p] = div(v, 2); end
end

# Primitives to help make these functions more generic later on

const IntDec = Union{Integer, Dec}

_getsign(x::Dec) = x.sgn
_getsign(x::Unsigned) = false
_getsign(x::Integer) = x < 0

_getcoeff(x::Dec) = x.val
_getcoeff(x::Integer) = BigInt(abs(x))
_getcoeff(x::BigInt) = x

_getint(x::Dec) = x.val
_getint(x::Integer) = abs(x)

_getscale(x::Dec) = x.pow
_getscale(x::Integer) = 0

_rnd10(d) = d > length(dectab) ? div(BigInt(10)^d, 2) : rndtab[d]
_pow10(d) = d > length(dectab) ? BigInt(10)^d : dectab[d]
_scale(x, d) = _getcoeff(x) * _pow10(d)

_eq(x::IntDec, y::IntDec) = _getint(x) == _getint(y)

_lt(x::IntDec, y::IntDec) = _getint(x) < _getint(y)

# Promotion rules
Base.promote_rule(::Type{Dec}, ::Type{<:Real}) = Dec

# override definitions in Base
Base.promote_rule(::Type{BigFloat}, ::Type{Dec}) = Dec
Base.promote_rule(::Type{BigInt}, ::Type{Dec}) = Dec

# Addition
function +(x::Dec, y::IntDec)
    xscl = _getscale(x)
    yscl = _getscale(y)
    # Quickly deal with zero case, make sure sign is correct in -0.0 + -0.0 case
    iszero(x) && return iszero(y) ? Dec(_getsign(x) & _getsign(y), 0, min(xscl, yscl)) : y
    iszero(y) && return x
    # Make both the same scale
    dscl = xscl - yscl
    if dscl == 0
        xval = _getcoeff(x)
        yval = _getcoeff(y)
    else
        xval, yval = dscl < 0 ? (_getint(x), _scale(y, -dscl)) : (_scale(x, dscl), _getint(y))
        abs(xscl) < abs(yscl) && (xscl = yscl)
    end
    # Simple case where signs are the same
    (xsgn = _getsign(x)) == _getsign(y) && return Dec(xsgn, xscl, xval + yval)
    # Signs are different, we need to subtract, possibly change sign
    (diff = xval - yval) < 0 ? Dec(!xsgn, xscl, -diff) : Dec(xsgn, xscl, diff)
end
+(x::Integer, y::Dec) = y + x

# Negation
-(x::Dec) = Dec(!x.sgn, x.pow, x.val)

# Subtraction
-(x::Dec, y::IntDec)  = x + -y
-(x::Integer, y::Dec) = -x + y

# Multiplication
*(x::Dec, y::IntDec) =
    Dec(xor(_getsign(x), _getsign(y)), _getscale(x) + _getscale(y), _getint(x) * _getint(y))
*(x::Integer, y::Dec) = y * x

# Inversion
function Base.inv(x::Dec)
    str = string(x)
    if str[1] == '-'
        str = str[2:end]
    end
    b = ('.' in str) ? length(split(str, '.')[1]) : 0
    c = round(BigInt(10)^(-x.pow + DIGITS) / x.val)
    q = (x.pow < 0) ? 1 - b - DIGITS : -b - DIGITS
    #normalize(Dec(x.sgn, q, c))
    Dec(x.sgn, q, c)
end

# Division
function /(x::Dec, y::Dec)
    x * inv(y)
end

# TODO exponentiation

@noinline argerr(str) = throw(ArgumentError("cannot parse \"$str\" as Dec"))

function _makedec(sgn, str, begpos, pos)
    #println("_makedec($sgn, \"$str\", $begpos, $pos) => $(pos-begpos) \"", str[begpos:pos-1], '"')
    (pos - begpos < 19
     ? Dec(sgn, 0, parse(Int64, str[begpos:pos-1]))
     : Dec(sgn, 0, parse(BigInt, str[begpos:pos-1])))
end

function _makedec(sgn, str, begpos, frcpos, frcend, scl)
    #println("_makedec($sgn, \"$str\", $begpos, $frcpos, $frcend, $scl) => \"",
    #str[begpos:frcpos-1], '"')
    v = begpos == frcpos ? BigInt(0) : parse(BigInt, str[begpos:frcpos-1])
    frcpos == frcend && return Dec(sgn, scl, v)
    diff = frcend - frcpos - 1
    #println("sgn = $sgn diff = $diff v = $v scl = $scl, \"", str[frcpos+1:frcend-1], '"')
    Dec(sgn, scl - diff, v * _pow10(diff) + parse(BigInt, str[frcpos+1:frcend-1]))
end

# Convert a string to a decimal, e.g. "0.01" -> Dec(0, -2, 1)
function Base.parse(::Type{Dec}, str::AbstractString)
    # start with +, -, 0, 1-9 or .
    (siz = sizeof(str)) == 0 && argerr(str)
    pos = 1
    c = str[pos]
    sgn = false
    if c == '+'
        (pos += 1) > siz && argerr(str)
        c = str[pos]
    elseif c == '-'
        (pos += 1) > siz && argerr(str)
        c = str[pos]
        sgn = true
    end
    # 'e' or 'E' is not allowed before having at least one digit
    (c == 'e' || c == 'E') && argerr(str)
    # may have any number of leading 0s, ignore them
    begpos = pos
    while c == '0'
        (pos += 1) > siz && return Dec(sgn, 0, 0)
        c = str[pos]
    end
    # Look for an immediate '.' (i.e. only fractional part)
    if c == '.'
        # if the string ends here, it must have had at least one leading zero to be well formed
        # "-.", "+." and "." are not valid, "-0.", "+0." and "0." are valid
        if (pos += 1) > siz
            begpos == siz && argerr(str)
            return Dec(sgn, 0, 0)
        end
        begpos = pos
        while (c = str[pos]) == '0'
            (pos += 1) > siz && return Dec(sgn, pos - begpos, 0)
        end
        # Now we can have 'e', 'E', or digits
        (c == 'e' || c == 'E') && return Dec(sgn, parse(Int32, str[pos+1:siz]) - (pos - begpos), 0)
        frcpos = pos
        while isdigit(c)
            # pos == siz && println("sgn = $sgn, zer = $zer str = \"", str[frcpos:siz], '"')
            (pos += 1) > siz && return Dec(sgn, begpos - pos, parse(BigInt, str[frcpos:siz]))
            c = str[pos]
        end
        (c == 'e' || c == 'E') || argerr(str)
        #println("sgn = $sgn v = ", parse(BigInt, str[frcpos:pos-1]), " scl = ",
        #parse(Int32, str[pos+1:siz]) - zer)
        return Dec(sgn, parse(Int32, str[pos+1:siz]) - (pos - begpos),
                   parse(BigInt, str[frcpos:pos-1]))
    end
    begpos = pos
    # We can have 'e', 'E', or digits 1-9 now
    while isdigit(c)
        (pos += 1) > siz && return _makedec(sgn, str, begpos, pos)
        c = str[pos]
    end
    # We can have '.', 'e', 'E' now
    frcpos = pos
    if c == '.'
        (pos += 1) > siz && return _makedec(sgn, str, begpos, pos-1)
        c = str[pos]
        while isdigit(c)
            (pos += 1) > siz && return _makedec(sgn, str, begpos, frcpos, pos, 0)
            c = str[pos]
        end
    end
    frcend = pos
    # We can now have only 'e' or 'E'
    (c == 'e' || c == 'E' || pos < siz) || argerr(str)
    # We can only have '+', '-', or digit(s) now
    c = str[pos += 1]
    (c == '+' || c == '-') && pos < siz && (c = str[pos += 1])
    # Must have at least one digit
    while isdigit(c)
        (pos += 1) > siz &&
            return _makedec(sgn, str, begpos, frcpos, frcend, parse(Int32, str[frcend+1:siz]))
        c = str[pos]
    end
    argerr(str)
end

dec(x::Real) = Dec(x)
dec(str::AbstractString) = parse(Dec, String(str))

const strzeros = repeat('0', 256)

function outzeros(io::IO, cnt::Integer)
    for i = 1:(cnt>>8)
        print(io, strzeros)
    end
    print(io, strzeros[1:(cnt&255)])
end

Base.tostr_sizehint(x::Dec) =
    Base.GMP.MPZ.sizeinbase(x.val, 10) + (x.pow == 0 ? ndigits(x.pow) + 3 : 1)

# Convert a decimal to a string
function Base.print(io::IO, x::Dec)
    x.sgn && print(io, '-')
    scl = x.pow
    if scl < 0
        c = string(x.val)
        len = sizeof(c)
        shift = scl + len
        if shift > 0
            print(io, c[1:shift], '.', c[(shift+1):len])
        elseif shift == 0
            print(io, "0.", c)
        elseif shift > -4
            print(io, "0.") ; outzeros(io, -shift) ; print(io, c)
        else
            print(io, c, 'e', scl)
        end
    else
        print(io, x.val)
        if scl > 2
            print(io, 'e', scl)
        elseif scl == 2
            print(io, "00")
        elseif scl == 1
            print(io, '0')
        end
    end
    # x.sgn ? x.pow == 0 ? print(io, '-', x.val) : print(io, '-', x.val, 'e', x.pow)
end

# Zero/one value
Base.zero(::Type{Dec}) = Dec(0, 0, 0)
Base.one(::Type{Dec})  = Dec(0, 0, 1)

# convert a decimal to any subtype of Real
#Base.convert(::Type{T}, x::Dec) where {T<:Real} =  parse(T, string(x))

(::Type{T})(x::Dec) where {T<:Integer} = parse(T, string(x))
(::Type{T})(x::Dec) where {T<:AbstractFloat} = parse(T, string(x))
(::Type{Rational{T}})(x::Dec) where {T<:Integer} = convert(Rational{T}, x)
(::Type{Rational{BigInt}})(x::Dec) = convert(Rational{BigInt}, x)

# fast case for Rationals
function Base.convert(::Type{T}, x::Dec) where {T<:Rational{<:Signed}}
    scl = _getscale(x)
    val = _getsign(x) ? -x.val : x.val
    scl < 0 ? T(val, _pow10(-scl)) : T(scl == 0 ? val : val * _pow10(scl), 1)
end

@noinline _inexact(T, x) = throw(InexactError(:convert, T, x))

# fast case for Integers
function Base.convert(::Type{T}, x::Dec) where {T<:Signed}
    scl = _getscale(x)
    val = _getsign(x) ? -x.val : x.val
    if scl < 0
        d, r = divrem(val, _pow10(-scl))
        r == 0 || _inexact(T, x)
        T(d)
    else
        T(scl == 0 ? val : val * _pow10(scl))
    end
end

# sign
Base.signbit(x::Dec) = x.sgn

# Equality

Base.iszero(x::Dec) = iszero(x.val)

# equals() now depends on == instead of the other way round.

function ==(x::Dec, y::IntDec)
    # Check if both are zero, regardless of sign
    iszero(x) && return iszero(y)
    iszero(y) && return false
    # Make sure signs are the same
    _getsign(x) == _getsign(y) || return false
    # If scales are the same, don't bother to equalize the scales
    (d = _getscale(x) - _getscale(y)) == 0 && return _eq(x, y)
    # Find out how much we need to multiply by to equalize the scales
    # Note: further optimization would use tables to see if the size (in limbs, or bits) of the two operands
    # could possibly be ==, without even doing the (somewhat expensive) scaling operation.
    d < 0 ? _eq(x, _scale(y, -d)) : _eq(_scale(x, d), y)
end

==(x::Integer, y::Dec) = y == x

function <(x::Dec, y::IntDec)
    # Check for both zeros, regardless of sign
    iszero(x) && return ifelse(iszero(y), false, !_getsign(y))
    iszero(y) && return _getsign(x)
    # Check signs
    xsgn = _getsign(x)
    ysgn = _getsign(y)
    xsgn == ysgn || return xsgn
    # If scales are the same, don't bother to equalize the scales
    (dscl = _getscale(x) - _getscale(y)) == 0 && return xor(_lt(x, y), ysgn)
    # Find out how much we need to multiply by to equalize the scales
    # Note: further optimization would use tables to see if the size (in limbs, or bits) of the two operands
    # are such that one is definitely larger than the other, without even doing the (somewhat expensive) scaling operation.
    xor(dscl < 0 ? _lt(x, _scale(y, -dscl)) : _lt(_scale(x, dscl), y), ysgn)
end

function <(x::Integer, y::Dec)
    # Check for both zeros, regardless of sign
    iszero(x) && return ifelse(iszero(y), false, !_getsign(y))
    iszero(y) && return _getsign(x)
    # Check signs
    xsgn = _getsign(x)
    ysgn = _getsign(y)
    xsgn == ysgn || return xsgn
    # If scales are the same, don't bother to equalize the scales
    xor((dscl = _getscale(y)) == 0 ? _lt(x, y) :
        (dscl < 0 ? _lt(_scale(x, -dscl), y) : _lt(x, _scale(y, dscl))), ysgn)
end

# 12345e-2 rounded to -1, you want 12e1
# 12345e-4 rounded to 0, you want to divrem by _pow10(-scl)
# 12345e-4 rounded to 3, you want to divrem by _pow10(-scl-digits)
# 12345e-4 rounded to >= 4, return x
# Rounding
function Base._round_digits(x::Dec, r::RoundingMode, digits::Integer, base)
    base != 10 && error("base=$base not implemented for type Dec yet")
    (scl = _getscale(x) + digits) >= 0 && return x
    d, r = divrem(_getint(x), _pow10(-scl))
    c = cmp(r, _rnd10(-scl))
    #println("scl = $scl, d = $d, r = $r, c = $c")
    if c != 0
        flg = (c < 0)
    elseif r === RoundNearest
        flg = iseven(d)
    else
        error("Rounding mode $r not implemented for type Dec yet")
    end
    Dec(_getsign(x), -digits, flg ? d : d + 1)
end

end # Decs


text:
function decompose(ln::ReferenceLine,x::Float64)
    @assert x ∈ ln
    a,b = vertices(ln)
    return LagrangeLine(x,a),LagrangeLine(x,b)
end

function decompose(tri::ReferenceTriangle,x)
    @assert x ∈ tri
    a,b,c = vertices(tri)
    return LagrangeTriangle(a,x,b), LagrangeTriangle(b,x,c), LagrangeTriangle(c,x,a)
end

function decompose(sq::ReferenceSquare,x)
    @assert x ∈ sq
    a,b,c,d = vertices(sq)
    return LagrangeTriangle(a,x,b), LagrangeTriangle(b,x,c), LagrangeTriangle(c,x,d), LagrangeTriangle(d,x,a)
end

function decompose(sq::ReferenceSquare,x,target::ReferenceSquare)
    @assert x ∈ sq
    a,b,c,d = vertices(sq)
    return  LagrangeRectangle(x,(0,x[2]),(0,0),(x[1],0)),
            LagrangeRectangle(x,(x[1],0),(1,0),(1,x[2])),
            LagrangeRectangle(x,(1,x[2]),(1,1),(x[1],1)),
            LagrangeRectangle(x,(x[1],1),(0,1),(0,x[2]))
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
@doc raw"""
Simple perturbation in electric field in (x,y,z) coordinates:
```math
E(x,y,z) = E_0 \, \begin{pmatrix}
0 \\
0 \\
\cos (2\pi z) \\
\end{pmatrix}
```

Parameters: `E₀`
"""
module EzCosZ

    import ..ElectromagneticFields
    import ..ElectromagneticFields: CartesianPerturbation, code
    import ..AnalyticCartesianField: X, Y, Z

    export EzCosZPerturbation

    const DEFAULT_E₀ = 1.0

    struct EzCosZPerturbation{T <: Number} <: CartesianPerturbation
        name::String
        E₀::T
        EzCosZPerturbation{T}(E₀::T) where T <: Number = new("EzCosZ", E₀)
    end

    EzCosZPerturbation(E₀::T=DEFAULT_E₀) where T <: Number = EzCosZPerturbation{T}(E₀)

    function init(E₀=DEFAULT_E₀)
        EzCosZPerturbation(E₀)
    end

    macro code(E₀=DEFAULT_E₀)
        code(init(E₀); escape=true)
    end

    function Base.show(io::IO, equ::EzCosZPerturbation)
        print(io, "Simple perturbation in electric field")
    end

    ElectromagneticFields.φ(x::AbstractVector, equ::EzCosZPerturbation) = equ.E₀ / (2π) * sin(2π * Z(x,equ))

    ElectromagneticFields.get_functions(::EzCosZPerturbation) = (X=X, Y=Y, Z=Z)

end