text:
"""
Types and methods related to parity.
"""
module Symmetries

export Parity, parity
export AngularMomentum, angularmomentum
export AngularSymmetry
export absdiff

using AtomicLevels
using HalfIntegers: HalfInteger, HalfInt, half

#
# Parity
# ------------------------------------------------------------------------------------------

"""
Represents the parity (even/odd) of a state.

Can be constructed from:

* `true` or `false` values, corresponding to even and odd, respectively.
* Characters `+` or `-`, corresponding to even and odd, respectively.
* The [`evenp`](@ref) and [`oddp`](@ref) constant instances.
"""
struct Parity
    even_parity :: Bool
end

"""
    even <: Parity

Represents even parity.
"""
const even = Parity(true)

"""
    odd <: Parity

Represents odd parity.
"""
const odd = Parity(false)

function Parity(c::Char)
    if c == '+'
        even
    elseif c == '-'
        odd
    else
        throw(ArgumentError("Invalid parity character '$c'. Must be '+' or '-'."))
    end
end

function Parity(x::Integer)
    if x == 1
        return Parity(true)
    elseif x == -1
        return Parity(false)
    else
        throw(ArgumentError("Invalid numeric parity value $x."))
    end
end

Parity(x::AtomicLevels.Parity) = Parity(x.p)

function Base.parse(::Type{Parity}, s::AbstractString)
    s = strip(s)
    length(s) == 1 || throw(Meta.ParseError("Can't parse '$s' into Parity (bad length)."))
    try
        Parity(first(strip(s)))
    catch e
        isa(e, ArgumentError) && throw(Meta.ParseError(e.msg))
        rethrow()
    end
end

Base.show(io::IO, p::Parity) = print(io, p.even_parity ? "+" : "-")

"""
    parity(object) -> Parity

Returns the parity (represented by an object of type [`Parity`](@ref)) of an object.
"""
function parity end

#
# Angular momentum
# ------------------------------------------------------------------------------------------

"""
Represent integer and half-integer angular momentum values.

Use `Rational` values to define half-integer angular momenta.
"""
struct AngularMomentum
    twoj :: UInt
    function AngularMomentum(j::Union{Integer, Rational, HalfInt})
        if j < zero(j)
            throw(ArgumentError("Angular momentum can't be negative: $j"))
        end

        twoj = if denominator(j) == 1
            2*numerator(j)
        elseif denominator(j) == 2
            numerator(j)
        else
            throw(ArgumentError("Angular momentum needs to be integer or half-integer: $j"))
        end

        return new(twoj)
    end
end

Base.convert(::Type{AngularMomentum}, j::Union{Integer, Rational}) = AngularMomentum(j)
Base.convert(::Type{Rational}, am::AngularMomentum) = convert(Rational{Int}, am)
Base.convert(::Type{Rational{Int}}, am::AngularMomentum) = Rational(Int(am.twoj), 2)
Base.convert(::Type{HalfInteger}, am::AngularMomentum) = half(convert(Int, am.twoj))

import Base: +, *
+(a::AngularMomentum, b::AngularMomentum) = AngularMomentum(convert(Rational, a) + convert(Rational, b))
function *(a::T, b::AngularMomentum) where {T <: Integer}
    a < zero(T) && throw(ArgumentError("Can't multiply AngularMomentum with a negative number."))
    AngularMomentum(a * b.twoj // 2)
end
absdiff(a::AngularMomentum, b::AngularMomentum) = AngularMomentum(abs(convert(Rational, a) - convert(Rational, b)))

Base.isless(a::AngularMomentum, b::AngularMomentum) = (a.twoj < b.twoj)

"""
    parse(::Type{AngularMomentum}, ::AbstractString)

Parses the string `s` into the corresponding `2J`-value. String can either be a number or a
fraction of the form `<x>/2`/ An empty string (also one including just whitespace) parses
into 0.
"""
function Base.parse(::Type{AngularMomentum}, s::AbstractString)
    if in('/', s)
        J_numerator, J_denominator = split(s, '/'; limit=2)
        @assert parse(Int, J_denominator) == 2
        AngularMomentum(parse(Int, J_numerator) // 2)
    elseif !isempty(strip(s))
        AngularMomentum(parse(Int, s))
    else
        AngularMomentum(0)
    end
end

Base.show(io::IO, am::AngularMomentum) = print(io, iseven(am.twoj) ? string(div(am.twoj, 2)) : "$(am.twoj)/2")

"""
    angularmomentum(object) -> AngularMomentum

Returns the angular momentum of an object, represented by an [`AngularMomentum`](@ref)
object.
"""
function angularmomentum end

#
# Angular symmetry (angular momentum + parity)
# ------------------------------------------------------------------------------------------

"""
Represents the angular symmetry (parity and angular momentum) of a state.
"""
struct AngularSymmetry
    angmom :: AngularMomentum
    parity :: Parity
    AngularSymmetry(angmom, parity::Parity) = new(angmom, parity)
end
AngularSymmetry(angmom, parity) = AngularSymmetry(angmom, Parity(parity))

angularmomentum(as::AngularSymmetry) = as.angmom
parity(as::AngularSymmetry) = as.parity

function Base.show(io::IO, as::AngularSymmetry)
    show(io, as.angmom)
    show(io, as.parity)
end

end


text:
module DNS

include("module_common.jl")

const _libcloud_dns_types = PyCall.PyNULL()
const _libcloud_dns_providers = PyCall.PyNULL()
const _libcloud_dns_base = PyCall.PyNULL()

immutable DNSDriver
    o::PyObject

    function DNSDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_dns_providers[:get_driver](provider)
        DNSDriver(cls(args...; kwargs...))
    end
    DNSDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::DNSDriver) = print(io, c.o[:__str__]())

immutable Zone
    o::PyObject

    id::Nullable{Compat.String}
    domain::Compat.String
    typ::Compat.String
    ttl::Nullable{Int}
    driver::DNSDriver
    extra::Dict

    function Zone(o::PyObject)
        new(o, o[:id], o[:domain], o[:type], o[:ttl], DNSDriver(o[:driver]), o[:extra])
    end
end
PyObject(o::Zone) = o.o
convert(::Type{Zone}, o::PyObject) =  Zone(o)
show(io::IO, o::Zone) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_dns_base, "Zone") Zone

immutable Record
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    typ::Compat.String
    data::Any
    zone::Zone
    driver::DNSDriver
    ttl::Nullable{Int}
    extra::Dict

    function Record(o::PyObject)
        new(o, o[:id], o[:name], o[:type], o[:data], Zone(o[:zone]), DNSDriver(o[:driver]), o[:ttl], o[:extra])
    end
end
PyObject(o::Record) = o.o
convert(::Type{Record}, o::PyObject) =  Record(o)
show(io::IO, o::Record) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_dns_base, "Record") Record

const _dnsdriver_fns = [
    :list_record_types, :iterate_zones, :list_zones, :iterate_records, :list_records,
    :get_zone, :get_record, :create_zone, :update_zone, :create_record, :update_record,
    :delete_zone, :delete_record, :export_zone_to_bind_format, :export_zone_to_bind_zone_file
]

const _zone_fns = [
    :list_records, :create_record, :update, :delete, :export_to_bind_format, :export_to_bind_zone_file
]

const _record_fns = [
    :update, :delete
]

for f in union(Set(_dnsdriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_dns_base, "DNSDriver", $sf) $(f)(dns::DNSDriver, args...; kwargs...) = dns.o[$(sf)](args...; kwargs...)
end

for f in _zone_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_dns_base, "Zone", $sf) $(f)(zone::Zone, args...; kwargs...) = zone.o[$(sf)](args...; kwargs...)
end

for f in _record_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_dns_base, "Record", $sf) $(f)(rec::Record, args...; kwargs...) = rec.o[$(sf)](args...; kwargs...)
end

# Initialize cloud dns
function __init__()
    copy!(_libcloud_dns_types, pyimport("libcloud.dns.types"))
    copy!(_libcloud_dns_providers, pyimport("libcloud.dns.providers"))
    copy!(_libcloud_dns_base, pyimport("libcloud.dns.base"))

    global const DNSProvider = pywrap(_libcloud_dns_types[:Provider])
    global const RecordType = pywrap(_libcloud_dns_types[:RecordType])

    _map_types(_libcloud_dns_base, (Zone, Record))
end

# types
export DNSProvider, DNSDriver, RecordType, Zone, Record
# dns driver functions
export list_record_types, iterate_zones, list_zones, iterate_records, list_records,
       get_zone, get_record, create_zone, update_zone, create_record, update_record,
       delete_zone, delete_record, export_zone_to_bind_format, export_zone_to_bind_zone_file
# zone functions
export list_records, create_record, update, delete, export_to_bind_format, export_to_bind_zone_file
# record functions
export update, delete

end # module DNS

ext: jl


text:
using DifferentialEquations, NLsolve

# Entire analysis for one press simulation (1 system)
function press_unit(par, p_strength)
    # Parameters
    u0 = [1.0, 1.5, 1.5]
    t_end = 5000 # way more than needed
    t_span = (0.0, t_end)
    t_start = 220
    p_start = 300.0
    len = 100000
    t_grid = range(t_start, t_end, length = len)
    t_press = range(p_start, t_end, length = len)
    Ω = par.Ω

    # functions
    press_event(u, t, integrator) = t ∈ p_start
    function forcing_press!(integrator)
        if integrator.t == p_start
            # addition for press
            integrator.p.K = p_strength + integrator.p.K_base
        end
        return
    end
    cb_press = DiscreteCallback(press_event, forcing_press!)

    # run fixed case first for active omnivory
    if par.pref == adapt_pref
        # compute ω
        # NB: This may often increases computation time as we may compute
        # results for passive omnivory twice but the code is much easier to
        # read. Here the time is not a major concern as it takes few minutes to
        # run.
        par_fix = deepcopy(par)
        par_fix.pref = fixed_pref
        prob_fix = ODEProblem(model!, u0, t_span, deepcopy(par_fix),
            tstops = p_start)
        sol_fix = solve(prob_fix, reltol = 1e-8, abstol = 1e-8,
            callback = cb_press)
        # sol_fix(t_grid)[1] starts at t_start
        eq = nlsolve(
            (du, u) -> model!(du, u, par_fix, 0.0),
            sol_fix(t_grid)[1]
        ).zero
        ### Solving for ω we have `ω = Ω * C^* / (Ω * C^* + (1 - Ω) * R^*)`
        ω = Ω * eq[2] / (Ω * eq[2] + (1 - Ω) * eq[1])
        par.ω = ω
    end

    par_after = deepcopy(par)
    par_after.K = par.K + p_strength
    prob = ODEProblem(model!, u0, t_span, deepcopy(par), tstops = p_start)
    sol = solve(prob, reltol = 1e-8, abstol = 1e-8, callback = cb_press)
    sol_grid = sol(t_grid)
    t_hit_eq = find_times_hit_equil_press(sol(t_press))

    eq_befor = find_eq(sol_grid[1], par)
    eq_after = find_eq(sol_grid[end], par_after)
    λ1 = λ1_stability(cmat(eq_after, par_after))
    λ1_imag = λ1_stability_imag(cmat(eq_after, par_after))
    react = ν_stability(cmat(eq_after, par_after))

    # Measure of Overshoot
    predator_OS = overshoot(sol, eq_after, 3, t_hit_eq[3], t_end)
    consumer_OS = overshoot(sol, eq_after, 2, t_hit_eq[2], t_end)
    resource_OS = overshoot(sol, eq_after, 1, t_hit_eq[1], t_end)

    # Calculate max-min metric
    predator_mm = min_max(sol, 3, t_hit_eq[3], t_end)
    consumer_mm = min_max(sol, 2, t_hit_eq[2], t_end)
    resource_mm = min_max(sol, 1, t_hit_eq[1], t_end)

    push!(
        [],
        λ1,
        -1 ./ λ1,
        λ1_imag,
        react,
        predator_OS,
        consumer_OS,
        resource_OS,
        predator_mm,
        consumer_mm,
        resource_mm,
        eq_befor,
        eq_after,
        par,
        sol
    )

end

# Perform the pulse for the 3 systems:
# - Food chain: "chain";
# - Passive omnivory: "passive";
# - Active omnivory: "active";
function press(base_par, Ω, p_strength)

    par_chain = deepcopy(base_par)
    par_chain.Ω = 0.0
    par_passive = deepcopy(base_par)
    par_passive.Ω = Ω
    par_passive.pref = fixed_pref
    par_active = deepcopy(par_passive)
    par_active.pref = adapt_pref

    push!(
        [],
        press_unit(par_chain, p_strength),
        press_unit(par_passive, p_strength),
        press_unit(par_active, p_strength)
    )

end

max_stars_repo_name: McCannLab/Labmnivory

max_issues_repo_name: McCannLab/Labmnivory

max_forks_repo_name: McCannLab/Labmnivory


text:
function right_nullspace(M::AbstractMatrix{T}) where {T}
    A, pivots = row_echelon_form(M)
    ncolsA = size(A, 2)
    length(pivots) == ncolsA && return zeros(T, ncolsA, 0)
    W = zeros(T, ncolsA, ncolsA - length(pivots))
    for (i, el) in enumerate(setdiff(1:ncolsA, pivots))
        W[el, i] += 1
        for (j, k) in enumerate(pivots)
            if j < el
                W[k, i] -= A[j, el]
            end
        end
    end
    return W
end

function left_nullspace(M::AbstractMatrix)
    return transpose(right_nullspace(transpose(M)))
end

function left_eigen(M::AbstractMatrix{T}) where {T<:FiniteFields.GF}
    @assert ==(size(M)...)
    eigen = Dict{T,typeof(M)}()
    cumdim = 0
    for i in T # brute force; TODO: use factorization of characteristic/minimal polynomials
        cumdim >= size(M, 1) && break
        #do left eigenspaces!
        basis = first(row_echelon_form!(left_nullspace(M - i * I)))
        if size(basis, 1) > 0 # nullity
            cumdim += size(basis, 1)
            eigen[i] = basis
        end
    end
    return eigen
end

function _find_l(M::AbstractMatrix)
    # this function should be redundant when defining a better structure for echelonized subspaces
    l = Int[]
    for i = 1:size(M, 2)
        j = findfirst(isone, @view M[length(l)+1:end, i])
        if j !== nothing
            push!(l, i)
        end
    end
    return l
end

# EigenSpaceDecomposition

function eigen_decomposition!(M::Matrix{T}) where {T<:FiniteFields.GF}
    eigspace_ptrs = Vector{Int}()
    eigen = left_eigen(M)
    sizehint!(eigspace_ptrs, length(eigen) + 1)
    push!(eigspace_ptrs, 1)
    for val in sort!(collect(keys(eigen))) #to get deterministic behaviour
        basis = eigen[val]
        dim = size(basis, 1)
        cd = eigspace_ptrs[end]
        ran = cd:cd+dim-1
        M[ran, :] = basis
        push!(eigspace_ptrs, cd + dim)
    end
    @assert eigspace_ptrs[end] == size(M, 1) + 1 "Matrix does not split over $T"
    return M, eigspace_ptrs
end

mutable struct EigenSpaceDecomposition{T<:FiniteFields.GF}
    basis::Matrix{T}
    eigspace_ptrs::Vector{Int}

    function EigenSpaceDecomposition(
        basis::Matrix{T},
        eigspace_ptrs::AbstractVector{<:Integer},
    ) where {T<:FiniteFields.GF}

        @assert eigspace_ptrs[1] == 1
        @assert eigspace_ptrs[end] == size(basis, 1) + 1
        return new{T}(basis, eigspace_ptrs)
    end
end

EigenSpaceDecomposition(M::Matrix{T}) where {T<:FiniteFields.GF} =
    EigenSpaceDecomposition(eigen_decomposition!(deepcopy(M))...)

function Base.show(
    io::IO,
    ::MIME"text/plain",
    esd::EigenSpaceDecomposition{T},
) where {T}
    println(io, tuple(diff(esd.eigspace_ptrs)...), "-splitting over ", T)
    print(io, esd.basis)
end

function Base.show(io::IO, esd::EigenSpaceDecomposition{T}) where {T}
    print(io, tuple(diff(esd.eigspace_ptrs)...), "-splitting over ", T)
end

Base.length(esd::EigenSpaceDecomposition) = length(esd.eigspace_ptrs) - 1

function Base.getindex(esd::EigenSpaceDecomposition, i::Int)
    @boundscheck 1 <= i <= length(esd)
    return esd.basis[esd.eigspace_ptrs[i]:esd.eigspace_ptrs[i+1]-1, :]
end

function Base.iterate(esd::EigenSpaceDecomposition, s = 1)
    s > length(esd) && return nothing
    first_last = esd.eigspace_ptrs[s]:esd.eigspace_ptrs[s+1]-1
    return (esd.basis[first_last, :], s + 1)
end

Base.eltype(esd::EigenSpaceDecomposition{T}) where {T} = Matrix{T}

LinearAlgebra.isdiag(esd::EigenSpaceDecomposition) =
    esd.eigspace_ptrs == 1:length(esd)+1

function refine(esd::EigenSpaceDecomposition{T}, M::Matrix{T}) where {T}
    nbasis = Array{T}(undef, 0, size(first(esd), 2))
    nptrs = [1]
    for (i, e) in enumerate(esd)
        if size(e, 1) > 1
            esd2, ptrs = eigen_decomposition!(e * @view M[:, _find_l(e)])
            nbasis = vcat(nbasis, esd2 * e)
            append!(nptrs, ptrs .+ (pop!(nptrs) - 1))
        else
            nbasis = vcat(nbasis, e)
            push!(nptrs, nptrs[end] + 1)
        end
    end
    return EigenSpaceDecomposition(nbasis, nptrs)
end


text:
function mehrmann(L,z;maxiter=10,tol=0., relax=1., x0=[],v=[], output=true)
  if output
      println("Launching Mehrmann...")
  end
  if x0==[]
    x0=ones(ComplexF64,size(L(0))[1])
  end
  if v==[]
    v=ones(ComplexF64,size(L(0))[1])
  end
  #normalize
  x0/=v'*x0

  active=L.active #TODO: better integrate activity into householder

  z0=complex(Inf)
  n=0
  flag=1
  try
    while abs(z-z0)>tol && n<maxiter
      if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
      z0=z
      u=L(z,0)\(L(z,1)*x0)
      #println(size(v)," ",size(x0)," ",size(u))
      z=z0-(v'*x0)/(v'*u)
      x0=u/(v'*u)
      #println("###",v'x0)
      n+=1
    end
  catch excp
    if output
        println("Error occured:")
        println(excp)
        println(typeof(excp))
        println("...aborted Mehrmann!")
    end

    flag=-2
  end

  #convergence checks
  if flag==1
    if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
    if n>=maxiter
      flag=-1
      if output; println("Warning: Maximum number of iterations has been reached!");end
    elseif abs(z-z0)<=tol
      flag=1
      if output; println("Solution has converged!"); end
    else
      if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
      flag=-3
      println(z)
    end
  end

  # if flag==1
  #   #compute left eigenvector
  #   #TODO: degeneracy
  #   lam,y = Arpack.eigs(L(z)',nev=1, sigma = 0,v0=x0)
  #   #normalize
  #   println("size:::$size(y)")
  #   x0./sqrt(x0'*x0)
  #   y./=conj(y'*L(z,1)*x0)
  # else
  #   y=[]
  # end
  y=[]

  return  Solution(L.params,x0,y,L.eigval), n , flag

end
##
function lancaster(L,z;maxiter=10,tol=0., relax=1., x0=[],y0=[], output=true)
  if output
      println("Launching Lancaster...")
  end
  if x0==[]
    x0=ones(ComplexF64,size(L(0))[1])
  end
  if y0==[]
    y0=ones(ComplexF64,size(L(0))[1])
  end

  z0=complex(Inf)
  n=0
  flag=1
  try
    while abs(z-z0)>tol && n<maxiter
            if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
      z0=z
      #solve
      ξ=L(z)\x0
      η=L(z)'\y0
      z=z0-(η'*L(z,0)*ξ)/(η'*L(z,1)*ξ)
      n+=1
    end
  catch excp
    if output
        println("Error occured:")
        println(excp)
        println(typeof(excp))
        println("...aborted Lancaster!")
    end

    flag=-2
  end

  #convergence checks
  if flag==1
    if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
    if n>=maxiter
      flag=-1
      if output; println("Warning: Maximum number of iterations has been reached!");end
    elseif abs(z-z0)<=tol
      flag=1
      if output; println("Solution has converged!"); end
    else
      if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
      flag=-3
      println(z)
    end
  end

  return  Solution(L.params,zeros(ComplexF64,length(x0)),[],L.eigval), n , flag
end

#%%
import LinearAlgebra
function juniper(L,z;maxiter=10,tol=0.,relax=1.,output=true)
    if output
        println("Launching Juniper...")
    end
  z0=complex(Inf)
  n=0
  flag=1
  try
    while abs(z-z0)>tol && n<maxiter
      if output;  println(n,"\t\t",abs(z-z0),"\t", z ); end
      z0=z
      #dz=-1/LinearAlgebra.tr(Array(L(z))\Array(L(z,1))) #TODO: better implementation for sparse types
      #QR=SparseArrays.qr(L(z))
      LU=SparseArrays.lu(L(z),check=false)
      L1=L(z,1)
      dz=0.0+0.0im
      for idx=1:size(LU)[1]
        dz+=(LU\Array(L1[:,idx]))[idx]
      end
      dz=-1/dz
      z=z0+relax*dz
      n+=1
    end
  catch excp
    if output
        println("Error occured:")
        println(excp)
        println(typeof(excp))
        println("...aborted Juniper!")
    end
    flag=-2
  end

  #convergence checks
  if flag==1
    if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
    if n>=maxiter
      flag=-1
      if output; println("Warning: Maximum number of iterations has been reached!");end
    elseif abs(z-z0)<=tol
      flag=1
      if output; println("Solution has converged!"); end
    else
      if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
      flag=-3
      println(z)
    end
  end

    return Solution(L.params,[],[],L.eigval), n, flag
end

##

function guettel(L,z;maxiter=10,tol=0., relax=1., x0=[],y0=[], output=true)
  if output
      println("Launching Guettel...")
  end
  if x0==[]
    x0=zeros(ComplexF64,size(L(0))[1])
    x0[1]=1
  end
  if y0==[]
    y0=zeros(ComplexF64,size(L(0))[1])
    y0[1]=1
  end
  #normalize
  x0/=sqrt(x0'*x0)
  y0/=sqrt(y0'*y0)
  z0=complex(Inf)
  n=0
  flag=1
  try
    while abs(z-z0)>tol && n<maxiter
      if output;  println(n,"\t\t",abs(z-z0),"\t", z ); end
      z0=z
      LU=SparseArrays.lu(L(z),check=false)
      x0=LU\(L(z,1)*x0)
      y0=LU'\(L(z,1)'*y0)
      #normalize
      x0/=sqrt(x0'*x0)
      y0/=sqrt(y0'*y0)
      #inner iteration
      idx=0
      z00=complex(Inf)
      while abs(z-z00)>tol && idx<10
        z00=z
        z=z-(y0'*L(z)*x0)/(y0'*L(z,1)*x0)
        idx+=1
      end
      n+=1
    end
  catch excp
    if output
        println("Error occured:")
        println(excp)
        println(typeof(excp))
        println("...aborted Guettel!")
    end
    flag=-2
  end
  #convergence checks
  if flag==1
    if output;  println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
    if n>=maxiter
      flag=-1
      if output; println("Warning: Maximum number of iterations has been reached!");end
    elseif abs(z-z0)<=tol
      flag=1
      if output; println("Solution has converged!"); end
    elseif isnan(z)
        flag=-5
        if output; println("Warning: computer arithmetics problem. Eigenvalue is NaN"); end
    else
      if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
      flag=-3
      println(z)
    end
  end
      return Solution(L.params,x0,y0,L.eigval), n, flag
end

## count poles and zeros

function count_eigvals(L,Γ,N;output=false)
  function integrand(z)
    z,w=z
    LU=SparseArrays.lu(L(z),check=false)
    L1=L(z,1)
    dz=0.0+0.0im
    for idx=1:size(LU)[1]
      dz+=(LU\Array(L1[:,idx]))[idx]
    end
    return dz*w
  end

  return gauss(0.0+0.0im,integrand,Γ,N,output)/2/pi/1.0im
end


text:
## Reactive signal update performance: o, h, l, c vectors from TimeArray values-component
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")

#s_ohlc = Reactive.Input(ohlc_BA.values[1,:])
s_ohlc = Reactive.Input(ohlc_BA[1])
s_close = lift(s -> s.values[4], s_ohlc, typ=Float64)

vsig = [s_close.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in 2:nmax
        #push!(s_ohlc, ohlc_BA.values[i,:])
        push!(s_ohlc, ohlc_BA[i])
        push!(vsig, s_close.value)
    end
end

println("OHLC updates o, h, l, c row from TimeArray (preserves timestamp of ohlc-row):")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: at least twice as long as ohlc_updates_vector_tavalues.jl

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
@inline issym(ex::Symbol) = true
@inline issym(ex) = false

@inline istuple(ex::Expr) = ex.head == :tuple
@inline istuple(ex) = false

@inline iscolon(ex::Expr) = ex.head == :call && first(ex.args) == :(:)
@inline iscolon(ex) = false

@inline iscollect(ex::Expr) = ex.head == Symbol("'")
@inline iscollect(ex) = false

@inline isleaf(ex::Expr) = istuple(ex) || iscolon(ex) || iscollect(ex)
@inline isleaf(ex::Symbol) = true
@inline isleaf(ex::Integer) = true
@inline isleaf(ex) = false

@inline getleft(ex::Expr) = ex.args[2]
@inline getright(ex::Expr) = ex.args[3]

leftmost(ex) = isleaf(ex) ? ex : leftmost(getleft(ex))
leftmostnode(ex) = isleaf(getleft(ex)) ? ex : leftmostnode(getleft(ex))

@inline isdup(ex::Integer) = true
@inline isdup(ex::Expr) = iscolon(ex) && getleft(ex) isa Integer
@inline isdup(ex) = false

@inline iscollectsym(ex) = iscollect(ex) && issym(first(ex.args))
@inline iscollecttup(ex) = iscollect(ex) && istuple(first(ex.args))

@inline istuplesym(ex) = istuple(ex) && all(issym, ex.args)
@inline istuplesymlike(ex) = istuple(ex) && all(issymlike, ex.args)

@inline issymlike(ex) = issym(ex) || iscollectsym(ex)
@inline istuplelike(ex) = istuplesymlike(ex) || iscollecttup(ex)

@inline hascollect(ex::Expr) = iscollect(ex) || any(iscollect, ex.args)
@inline hascollect(ex) = false

removecollect(ex) = ex
function removecollect(ex::Expr)
  if hascollect(ex)
    if iscollect(ex)
      return first(ex.args)
    else
      ret = copy(ex)
      for (i, e) ∈ enumerate(ex.args)
        if iscollect(e)
          ret.args[i] = first(e.args)
        end
      end
      return ret
    end
  else
    return ex
  end
end

collectcollect(ex) = :(())
function collectcollect(ex::Expr)
  if iscollect(ex)
    return first(ex.args)
  else
    ret = Any[]
    for (i, e) ∈ enumerate(ex.args)
      if iscollect(e)
        push!(ret, first(e.args))
      end
    end
    if length(ret) == 1
      return ret[1]
    else
      return Expr(:tuple, ret...)
    end
  end
end

const legalsym = (:(=>),:→)

@inline islegal(ex::Symbol) = true
@inline islegal(ex::Integer) = ex > 0
@inline islegal(ex) = false
function islegal(ex::Expr)
  global legalsym

  if istuple(ex)
    all(issymlike, ex.args)
  elseif iscolon(ex)
    length(ex.args) == 3 && isdup(ex) ?
      islegal(getleft(ex)) &&
      (issymlike(getright(ex)) ||
       (istuplelike(getright(ex)) && islegal(getright(ex)))) :
      (issymlike(getleft(ex)) ||
       (istuplelike(getleft(ex)) && islegal(getleft(ex)))) &&
      (issymlike(getright(ex)) ||
       (istuplelike(getright(ex)) && islegal(getright(ex))))
  elseif iscollect(ex)
    length(ex.args) == 1 && (issym(first(ex.args)) || istuplesym(first(ex.args)))
  else
    ex.head == :call &&
      first(ex.args) ∈ legalsym &&
      length(ex.args) == 3 &&
      islegal(getleft(ex)) &&
      islegal(getright(ex))
  end
end

struct Code
  in
  out
  lines::Vector{Tuple}
end

Code(in, out) = Code(in, out, [(out, in)])

function Base.show(io::IO, code::Code)
  println(io, "Code: input: $(code.in), output: $(code.out)")
  println(io, "body: ")
  for l in code.lines
    println("\t$(l[1]) = $(l[2])")
  end
  io
end

@inline out(code::Code) = code.out
@inline out(c) = c

function add(code1::Code, code2::Code)
  in = code1.in
  out = code2.out
  lines = cat(code1.lines, code2.lines; dims=1)
  Code(in, out, lines)
end

function add(code::Code, newline)
  in = code.in
  out = newline
  lines = [code.lines..., (out, code.out)]
  Code(in, out, lines)
end

function add(left, right)
  in = left
  out = right
  Code(in, out)
end

function duplicate(code::Code, n::Int)
  in = code.in
  out = code.out

  tp_cls = [(first(code.lines)[1], out), code.lines[2:end]...]
  tp_cls = [tp_cls for i = 1:(n-1)]

  Code(in, out, cat(code.lines, tp_cls...; dims=1))
end

function duplicate(c, n::Int)
  in = c
  out = c
  lines = [(c, c) for i = 1:n]
  Code(in, out, lines)
end

@inline getint(ex::Int) = ex
@inline getint(ex::Expr) = getleft(ex) isa Int ? getleft(ex) : getright(ex)
@inline getsym(ex::Expr) = getleft(ex) isa Int ? getright(ex) : getleft(ex)

@inline getin(ex::Expr) = iscolon(ex) ? getleft(ex) : ex
@inline getin(x::Symbol) = x
@inline getout(ex::Expr) = iscolon(ex) ? getright(ex) : ex
@inline getout(x::Symbol) = x

function _to_code(node)
  if !isleaf(getleft(node))
    pre_code = _to_code(getleft(node))
  else
    pre_code = getleft(node)
  end

  rL = leftmost(getright(node)) #Symbol(or tuple) or Int
  if isdup(rL)#rL isa Int
    rightnode = getright(node)
    if isdup(rightnode)#getright(node) isa Int
      node.args[3] = rL isa Int ? out(pre_code) : getsym(rL)
    else
      leftmostnode(rightnode).args[2] = rL isa Int ? out(pre_code) : getsym(rL)
    end
    if getint(rL) == 1 #don't need duplicate
      code = pre_code
    else
      code = duplicate(pre_code, getint(rL))
    end
  else
    code = add(pre_code, rL)
  end

  if !isleaf(getright(node))
    next_code = _to_code(getright(node))
    code = add(code, next_code)
  end

  code
end

function _postcode(code::Code)
  in = getin(code.in)
  out = getin(code.out)
  lines = map(code.lines) do (vab, farg)
    (getin(vab), getout(farg))
  end
  Code(in, out, lines)
end

function to_code(ex::Expr)
  code = _to_code(ex)
  _postcode(code)
end


text:
function dummy_project(f; keep_files=false)
    folder = joinpath(@__DIR__,"dummy_project")
    isdir(folder) && rm(folder, recursive=true)
    mkdir(folder)
    initialize_project(folder)
    try
        f(folder)
    catch e
        keep_files|| rm(folder,recursive=true)
        rethrow(e)
    end
    keep_files|| rm(folder,recursive=true)
    @quickactivate
end

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
mutable struct NestedTimeIntervalDisplayProperty
    unix_starttime::Float64

end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,U,p)
    rhoU = rho.*U
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

const TOL = 1e-16
"Approximation parameters"
N = 3 # The order of approximation
K = 50
T = 1.0
T = 6.0
#T = 0.0039

# # Sod shocktube
# const γ = 1.4
# const Bl = -0.5
# const Br = 0.5
# const rhoL = 1.0
# const rhoR = 0.125
# const pL = 1.0
# const pR = 0.1
# const xC = 0.0
# const GIFINTERVAL = 20
# T = 0.2

# Leblanc shocktube
const γ = 5/3
const Bl = 0.0
const Br = 9.0
const rhoL = 1.0
const rhoR = 0.001
const pL = 0.1
const pR = 1e-7
const xC = 3.0
const GIFINTERVAL = 600
T = 6.0

"Mesh related variables"
VX = LinRange(Bl,Br,K+1)
EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K))

"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
Vq = vandermonde_1D(N,rq)/VDM
M = Vq'*diagm(wq)*Vq
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# operators
Qr = M*Dr
B = zeros(N+1,N+1)
B[1,1] = -1
B[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Qr0 = S0+1/2*B

# Drop zeros
Qr = Matrix(droptol!(sparse(Qr),TOL))
Qr0 = Matrix(droptol!(sparse(Qr0),TOL))
B = Matrix(droptol!(sparse(B),TOL))
L = Matrix(droptol!(sparse(L),TOL))
S0 = Matrix(droptol!(sparse(S0),TOL))
S = Matrix(droptol!(sparse((Qr-Qr')/2),TOL))

rf = [-1.0;1.0]
nrJ = [-1.0;1.0]
Vf = vandermonde_1D(N,rf)/VDM

"""High order mesh"""
x = V1*VX[transpose(EToV)]
xf = Vf*x
mapM = reshape(1:2*K,2,K)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

# """Periodic"""
# mapP[1] = mapM[end]
# mapP[end] = mapP[1]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

# """Low order mesh"""
# VX0 = [x[1:end-1,:][:]; x[end]]
# #EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K)
# EToV0 = transpose(reshape(sort([1:N*K; 2:N*K+1]),2,N*K))
# x0 = VX0[transpose(EToV0)]
# xf0 = x0
# mapM0 = reshape(1:2*N*K,2,N*K)
# mapP0 = copy(mapM0)
# mapP0[1,2:end] .= mapM0[2,1:end-1]
# mapP0[2,1:end-1] .= mapM0[1,2:end]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Initial condition"""
rho_x(x) = (x <= xC) ? rhoL : rhoR
u_x(x) = 0.0
p_x(x) = (x <= xC) ? pL : pR

rho = @. rho_x(x)
u = @. u_x(x)
p = @. p_x(x)
U = primitive_to_conservative_hardcode(rho,u,p)

function limiting_param(U_low, P_ij)
    l = 1.0
    # Limit density
    if U_low[1] + P_ij[1] < -TOL
        l = min(abs(U_low[1])/(abs(P_ij[1])+1e-14), 1.0)
    end

    # limiting internal energy (via quadratic function)
    a = P_ij[1]*P_ij[3]-1.0/2.0*P_ij[2]^2
    b = U_low[3]*P_ij[1]+U_low[1]*P_ij[3]-U_low[2]*P_ij[2]
    c = U_low[3]*U_low[1]-1.0/2.0*U_low[2]^2

    l_eps_ij = 1.0
    if b^2-4*a*c >= 0
        r1 = (-b+sqrt(b^2-4*a*c))/(2*a)
        r2 = (-b-sqrt(b^2-4*a*c))/(2*a)
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

function flux_lowIDP(U_i,U_j,f_i,f_j,c_ij,wavespd)
    return c_ij*(f_i+f_j)-abs(c_ij)*wavespd*(U_j-U_i)
end

function flux_high(f_i,f_j,c_ij)
    return c_ij*(f_i+f_j)
end

function rhs_IDP(U,K,N,wq,S,S0,dt,Mlump_inv)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_low       = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_high      = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_low_P     = [zeros(2),zeros(2),zeros(2)] # 1: left boundary, 2: right boundary
    F_high_P    = [zeros(2),zeros(2),zeros(2)]
    L           = zeros(N+1,N+1) # Array of limiting params
    L_P         = zeros(K-1) # limiting params at left/right boundary
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    # Determine dt
    wavespd_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
        end
    end

    for k = 1:K
        # Assemble matrix of low and high order algebraic fluxes
        # interior of the element
        for i = 1:N+1
            for j = 1:N+1
                if i != j # skip diagonal
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k])
                    for c = 1:3
                        F_low[c][i,j]  = flux_lowIDP(U[c][i,k],U[c][j,k],flux[c][i,k],flux[c][j,k],S0[i,j],wavespd)
                        F_high[c][i,j] = flux_high(flux[c][i,k],flux[c][j,k],S[i,j])
                    end
                end
            end
        end

        # Assemble matrix of low and high order algebraic fluxes
        # interface of the element
        U_left = (k == 1)  ? [rhoL; 0.0; pL/(γ-1)] : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        f_left = (k == 1)  ? [0.0; pL; 0.0]        : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right = (k == K) ? [rhoR; 0.0; pR/(γ-1)] : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        f_right = (k == K) ? [0.0; pR; 0.0]        : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]))

        for c = 1:3
            F_low_P[c][1] = flux_lowIDP(U[c][1,k],U_left[c],flux[c][1,k],f_left[c],-0.5,wavespd_l)
            F_low_P[c][2] = flux_lowIDP(U[c][end,k],U_right[c],flux[c][end,k],f_right[c],0.5,wavespd_r)

            F_high_P[c][1] = flux_high(flux[c][1,k],f_left[c],-0.5)
            F_high_P[c][2] = flux_high(flux[c][end,k],f_right[c],0.5)
        end

        # Calculate limiting parameters over interior of the element
        P_ij = zeros(3,1)
        # TODO: redundant
        U_low =  [zeros(N+1,1),zeros(N+1,1),zeros(N+1,1)]
        for c = 1:3
            U_low[c] .= sum(-F_low[c],dims=2)
            U_low[c][1] -= F_low_P[c][1]
            U_low[c][end] -= F_low_P[c][2]
            U_low[c] .= U[c][:,k]+dt/J*Mlump_inv*U_low[c]
        end
        for i = 1:N+1
            lambda_j = (i >= 2 && i <= N) ? 1/N : 1/(N+1)
            m_i = J*wq[i]
            for j = 1:N+1
                if i != j
                    for c = 1:3
                        P_ij[c] = dt/(m_i*lambda_j)*(F_low[c][i,j]-F_high[c][i,j])
                    end
                    L[i,j] = limiting_param([U_low[1][i]; U_low[2][i]; U_low[3][i]],P_ij)
                end
            end
        end

        # Symmetrize limiting parameters
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    l_ij = min(L[i,j],L[j,i])
                    L[i,j] = l_ij
                    L[j,i] = l_ij
                end
            end
        end

        # Limiting parameters at the interface
        if k <= K-1
            lambda_j = 1/(N+1)
            rhs_low = zeros(3,1)
            U_lowP = [U_low[1][end]; U_low[2][end]; U_low[3][end]]
            # current element
            m_i = wq[end]*J
            for c = 1:3
                P_ij[c] = dt/(m_i*lambda_j)*(F_low_P[c][2]-F_high_P[c][2])
            end

            L_P[k] = limiting_param(U_lowP, P_ij)

            # next element
            U_lowP = zeros(3,1)
            P_ij = zeros(3,1)
            m_i = wq[1]*J
            lambda_j = 1/(N+1)

            # TODO: Very awkward... need to repeat the calculation for the next element
            F_low_tmp = [zeros(1,N+1),zeros(1,N+1),zeros(1,N+1)]
            F_high_tmp = [zeros(1,N+1),zeros(1,N+1),zeros(1,N+1)]
            F_lowP_tmp = zeros(3,1)
            F_highP_tmp = zeros(3,1)
            i = 1
            for j = 1:N+1
                if i != j # skip diagonal
                    wavespd = max(wavespd_arr[i,k+1],wavespd_arr[j,k+1])
                    for c = 1:3
                        F_low_tmp[c][i,j] = flux_lowIDP(U[c][i,k+1],U[c][j,k+1],flux[c][i,k+1],flux[c][j,k+1],S0[i,j],wavespd)
                        F_high_tmp[c][i,j] = flux_high(flux[c][i,k+1],flux[c][j,k+1],S[i,j])
                    end
                end
            end

            U_left = [U[1][end,k]; U[2][end,k]; U[3][end,k]]
            f_left = [flux[1][end,k]; flux[2][end,k]; flux[3][end,k]]
            wavespd_l = max(wavespd_arr[1,k+1],wavespeed_1D(U_left[1],U_left[2],U_left[3]))
            for c = 1:3
                F_lowP_tmp[c] = flux_lowIDP(U[c][1,k+1],U_left[c],flux[c][1,k+1],f_left[c],-0.5,wavespd_l)
                F_highP_tmp[c] = flux_high(flux[c][1,k+1],f_left[c],-0.5)
            end

            U_low = zeros(3,1)
            for c = 1:3
                U_low[c] = sum(-F_low_tmp[c])-F_lowP_tmp[c][1]
                U_low[c] = U[c][1,k+1]+dt*1/J*Mlump_inv[1]*U_low[c]
            end

            for c = 1:3
                P_ij[c] = dt/(m_i*lambda_j)*(F_lowP_tmp[c]-F_highP_tmp[c])
            end

            Ltmp = limiting_param([U_low[1]; U_low[2]; U_low[3]],P_ij)
            L_P[k] = min(Ltmp,L_P[k])

        end

        # construct rhs
        for c = 1:3
            # With limiting
            rhsU[c][:,k] = sum((L.-1).*F_low[c] - L.*F_high[c],dims=2)

            if k > 1
                rhsU[c][1,k] += (L_P[k-1]-1)*F_low_P[c][1] - L_P[k-1]*F_high_P[c][1]
            else
                rhsU[c][1,k] += -F_low_P[c][1]
            end
            if k < K
                rhsU[c][N+1,k] += (L_P[k]-1)*F_low_P[c][2] - L_P[k]*F_high_P[c][2]
            else
                rhsU[c][N+1,k] += -F_low_P[c][2]
            end

            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end

    end

    return rhsU
end

function rhs_IDPlow(U,K,N,Mlump_inv,p,flux,J)
    dfdx = (zeros(N+1,K),zeros(N+1,K),zeros(N+1,K))
    for i = 2:K*(N+1)-1
        for c = 1:3
            dfdx[c][i] = 1/2*(flux[c][mod1(i+1,K*(N+1))] - flux[c][mod1(i-1,K*(N+1))])
        end
    end
    dfdx[1][1] = 1/2*(flux[1][2] - 0.0)
    dfdx[2][1] = 1/2*(flux[2][2] - pL)
    dfdx[3][1] = 1/2*(flux[3][2] - 0.0)
    dfdx[1][end] = 1/2*(0.0 - flux[1][end-1])
    dfdx[2][end] = 1/2*(pR - flux[2][end-1])
    dfdx[3][end] = 1/2*(0.0 - flux[3][end-1])

    visc = (zeros(N+1,K),zeros(N+1,K),zeros(N+1,K))
    for i = 2:K*(N+1)-1
        wavespd_curr = wavespeed_1D(U[1][i],U[2][i],U[3][i])
        wavespd_R = wavespeed_1D(U[1][mod1(i+1,K*(N+1))],U[2][mod1(i+1,K*(N+1))],U[3][mod1(i+1,K*(N+1))])
        wavespd_L = wavespeed_1D(U[1][mod1(i-1,K*(N+1))],U[2][mod1(i-1,K*(N+1))],U[3][mod1(i-1,K*(N+1))])
        dL = 1/2*max(wavespd_curr,wavespd_L)
        dR = 1/2*max(wavespd_curr,wavespd_R)
        for c = 1:3
            visc[c][i] = dL*(U[c][mod1(i-1,K*(N+1))]-U[c][i]) + dR*(U[c][mod1(i+1,K*(N+1))]-U[c][i])
        end
    end

    # i = 1
    wavespd_curr = wavespeed_1D(U[1][1],U[2][1],U[3][1])
    wavespd_R = wavespeed_1D(U[1][2],U[2][2],U[3][2])
    wavespd_L = wavespeed_1D(rhoL,0.0,pL/(γ-1))
    dL = 1/2*max(wavespd_curr,wavespd_L)
    dR = 1/2*max(wavespd_curr,wavespd_R)
    visc[1][1] = dL*(rhoL-U[1][1]) + dR*(U[1][2]-U[1][1])
    visc[2][1] = dL*(0.0-U[2][1]) + dR*(U[2][2]-U[2][1])
    visc[3][1] = dL*(pL/(γ-1)-U[3][1]) + dR*(U[3][2]-U[3][1])

    # i = end
    wavespd_curr = wavespeed_1D(U[1][end],U[2][end],U[3][end])
    wavespd_R = wavespeed_1D(rhoR,0.0,pR/(γ-1))
    wavespd_L = wavespeed_1D(U[1][end-1],U[2][end-1],U[3][end-1])
    dL = 1/2*max(wavespd_curr,wavespd_L)
    dR = 1/2*max(wavespd_curr,wavespd_R)
    visc[1][end] = dL*(U[1][end-1]-U[1][end]) + dR*(rhoR-U[1][end])
    visc[2][end] = dL*(U[2][end-1]-U[2][end]) + dR*(0.0-U[2][end])
    visc[3][end] = dL*(U[3][end-1]-U[3][end]) + dR*(pR/(γ-1)-U[3][end])

    rhsU = (x->1/J*Mlump_inv*x).(.-dfdx.+visc)
    return rhsU
end

function rhs_high(U,K,N,Mlump_inv,S,p,flux,J)
    dfdx = (zeros(N+1,K),zeros(N+1,K),zeros(N+1,K))
    for k = 1:K
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    for c = 1:3
                        dfdx[c][i,k] += S[i,j]*(flux[c][i,k]+flux[c][j,k])
                    end
                end
            end
        end

        fLB = [0.0;pL;0.0]
        fRB = [0.0;pR;0.0]
        # At interfaces
        for c = 1:3
            fL = (k > 1) ? flux[c][end,k-1] : fLB[c]
            fR = (k < K) ? flux[c][1,k+1]   : fRB[c]
            dfdx[c][1,k] += -1.0/2.0*(flux[c][1,k]+fL)
            dfdx[c][end,k] += 1.0/2.0*(flux[c][end,k]+fR)
        end
    end

    rhsU = (x->-1/J*Mlump_inv*x).(dfdx)
    return rhsU
end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x))]

# # Forward Euler
# while t < T
#     dt = 0.0001
#     rhsU = rhs_IDP(U,K,N,wq,S,S0,dt,Mlump_inv)
#     @. U = U + dt*rhsU
#     global t = t + dt
#     println("Current time $t with time step size $dt, and final time $T")
# end

# Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM
# gr(size=(300,300),ylims=(0,1.2),legend=false,markerstrokewidth=1,markersize=2)
# plt = plot(Vp*x,Vp*U[1])

Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM
gr(size=(300,300),ylims=(0,1.2),legend=false,markerstrokewidth=1,markersize=2)
plot()

dt = 0.0001
Nsteps = Int(T/dt)
@gif for i = 1:Nsteps
    rhsU = rhs_IDP(U,K,N,wq,S,S0,dt,Mlump_inv)
    @. U = U + dt*rhsU
    global t = t + dt
    println("Current time $t with time step size $dt, and final time $T")
    if i % GIFINTERVAL == 0
        plot(Vp*x,Vp*U[1])
    end
end every GIFINTERVAL

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
#Local function for checking a BootInput is valid for a call to dbootvecinds
function check_bi_for_dbootvecinds(bi::BootInput)::Bool
    bi.numobs < 1 && error("Number of observations must be strictly positive")
    bi.numresample < 1 && error("Number of resamples must be strictly positive")
    bi.numobsperresample < 1 && error("Number of observations per resample must be strictly positive")
    bi.blocklength <= 0.0 && error("Block length must be strictly positive")
    isnan(bi.blocklength) && error("Block length is set to NaN")
    isinf(bi.blocklength) && error("Block length is infinite")
    return true
end

"""
dbootinds_one(bi::BootInput)::Vector{Int}
dbootinds_one(data::T; kwargs...)::Vector{Int}

Returns a single resampling index that, when used to index the original dataset,
will provide a single resampled dataset.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.
"""
(dbootinds_one(bi::BootInput{BootIID})::Vector{Int}) = rand(1:bi.numobs, bi.numobsperresample)
function dbootinds_one(bi::BootInput{BootStationary})::Vector{Int}
    bi.blocklength <= 1.0 && return rand(1:bi.numobs, bi.numobsperresample)
    inds = zeros(Int, bi.numobsperresample)
    geo1 = Geometric(1 / bi.blocklength)
    (c, geodraw) = (1, 1)
    for n = 1:bi.numobsperresample
        if c == geodraw #Start a new block
            inds[n] = rand(1:bi.numobs)
            geodraw = rand(geo1) + 1
            c = 1
        else #Next obs in existing block
            inds[n-1] == bi.numobs ? (inds[n] = 1) : (inds[n] = inds[n-1] + 1)
            c += 1
        end
    end
    return inds
end
function dbootinds_one(bi::BootInput{BootMoving})::Vector{Int}
    bl = ceil(Int, bi.blocklength)
    bl == 1 && return rand(1:bi.numobs, bi.numobsperresample)
    inds = zeros(Int, bi.numobsperresample)
    blockstart_ub = max(1, bi.numobs-bl+1)
    for n = 1:bl:bi.numobsperresample
        inds[n] = rand(1:blockstart_ub) #Start of block
        for s = n+1:min(n+bl-1, bi.numobsperresample) #Iterate through block (use of min avoids bounds error)
            inds[s] = inds[s-1] + 1
        end
    end
	return inds
end
function dbootinds_one(bi::BootInput{BootCircular})::Vector{Int}
    bl = ceil(Int, bi.blocklength)
    bl == 1 && return rand(1:bi.numobs, bi.numobsperresample)
    inds = zeros(Int, bi.numobsperresample)
    for n = 1:bl:bi.numobsperresample
        inds[n] = rand(1:bi.numobs) #Start of block
        for s = n+1:min(n+bl-1, bi.numobsperresample) #Iterate through block (use of min avoids bounds error)
            inds[s-1] == bi.numobs ? (inds[s] = 1) : (inds[s] = inds[s-1] + 1)
        end
    end
	return inds
end
function dbootinds_one(bi::BootInput{BootNoOverlap})::Vector{Int}
    bl = ceil(Int, bi.blocklength)
    bl == 1 && return rand(1:bi.numobs, bi.numobsperresample)
    inds = zeros(Int, bi.numobsperresample)
	blockstart_ub = max(1, bi.numobs-bl+1)
	blockstartvalues = collect(1:bl:blockstart_ub) #Build valid set of start indices for any block
    length(blockstartvalues) == 1 && error("Not enough observations to perform non-overlapping block bootstrap given block length: num obs = $(bi.numobs), block length = $(bl)")
    (blockstartvalues[end] + bl - 1 > bi.numobs) && pop!(blockstartvalues)
    for n = 1:bl:bi.numobsperresample
        inds[n] = blockstartvalues[rand(1:length(blockstartvalues))] #Start of block
        for s = n+1:min(n+bl-1, bi.numobsperresample) #Iterate through block (use of min avoids bounds error)
            inds[s] = inds[s-1] + 1
        end
    end
	return inds
end
#(dbootinds_one(bi::BootInput{BootTapered})::Vector{Int}) = error("Routines for the tapered block bootstrap are currently not completed. Users interested in contributing should check the package github page.")
(dbootinds_one(bi::BootInput{BootDummy})::Vector{Int}) = error("Logic fail. It should not have been possible to call this method. Please file an issue with full stacktrace.")
(dbootinds_one(data, bi::BootInput)::Vector{Int}) = dbootinds_one(bi)
(dbootinds_one(data ; kwargs...)::Vector{Int}) = dbootinds_one(data, BootInput(data ; kwargs...))

"""
    dbootinds(data::T ; bi::BootInput)::Vector{Vector{Int}}
    dbootinds(data::T ; kwargs...)::Vector{Vector{Int}}

Each inner vector of the returned Vector{Vector{Int}} provides indices that, when used to
index the original dataset, will provide a single resampled dataset.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.

Please use dbootinds_one if you only want to obtain a single Vector{Int} resampling index.
"""
(dbootinds(bi::BootInput)::Vector{Vector{Int}}) = check_bi_for_dbootvecinds(bi) ? [ dbootinds_one(bi) for n = 1:bi.numresample ] : error("Logic fail in check_bi_for_dbootvecinds. Please file an issue.")
(dbootinds(data, bi::BootInput)::Vector{Vector{Int}}) = dbootinds(bi)
(dbootinds(data ; kwargs...)::Vector{Vector{Int}}) = dbootinds(data, BootInput(data ; kwargs...))


text:
include("common.jl")

cxx"""
#include "lldb/Target/Process.h"
#include "lldb/lldb-enumerations.h"
#include "/Users/kfischer/julia/deps/llvm-svn/tools/lldb/source/Plugins/Process/Utility/RegisterContextMemory.h"
#include "/Users/kfischer/julia/deps/llvm-svn/tools/lldb/source/Plugins/Process/Utility/ThreadMemory.h"
"""

lldb_exec(dbg,"log enable lldb")

tlist = icxx"$target->m_process_sp->GetThreadList();"
process = icxx"$target->m_process_sp.get();";
thread = tlist[0]
reg_ctx = icxx"$thread->GetRegisterContext().get();"
tf=0x80193f68

mreg_ctx = icxx"""
new RegisterContextMemory(*$thread,0,$reg_ctx,$tf,38*sizeof(uint32_t));
"""

thread = icxx"""
auto T2 = new ThreadMemory (*$process,
                  0x20,
                  "user level",
                  "",
                  $tf);
T2->SetRegisterContext($mreg_ctx);
return T2;
"""
#=
TSP = icxx" lldb::ThreadSP T($thread); return T; "
icxx" $thread->GetStatus (*$dbg->GetOutputFile().get(),0,10,0); ";

=#

icxx"""
lldb_private::RegisterValue reg_value;
reg_value.SetType(lldb::eRegisterValueGCC);
$mreg_ctx->ReadRegister($(mreg_ctx[29]),reg_value);
"""

data = Array(UInt32,47)
icxx"
lldb_private::Error error;
$process->ReadMemory($tf,$(pointer(data)),$(sizeof(data)),error);
"


text:
module BoxTest

using Projections, Test

@testset "box" begin
y = [-2.5, -1.5, -0.5, 0.0, 0.5, 1.5, 2.5]
a = fill(-1.0, length(y))
b = fill(1.0, length(y))
@test all(project(Box(a, b), y) .== [clamp(y[i], a[i], b[i]) for i in 1:length(y)])
end

end


text:
module LOBPCG

#export locg, LOCGBuffersSimple, LOCGBuffersGeneral
export lobpcg

include("LOBPCG3.jl")

end # module


text:
using InteractionNetworkModels, StructuredDistances, Distributions
using Plots, Random

E = [
    [1,2,1,2],
    [1,2,1],
    [3,4,3]
]
d = MatchingDistance(FastLCS(101))
K_inner, K_outer = (DimensionRange(1,10), DimensionRange(1,25))
model = SIM(
    E, 3.5,
    d,
    1:10,
    K_inner, K_outer
)

# Auxiliary sampler
β = 0.7
mcmc_move = InvMcmcMixtureMove(
    (
        EditAllocationMove(ν=3),
        InsertDeleteMove(ν=3, len_dist=TrGeometric(0.8, 1, model.K_inner.u))
    ),
    (β, 1-β)
)
mcmc_sampler = InvMcmcSampler(
    mcmc_move,
    lag=30, burn_in=3000
)

@time x = mcmc_sampler(
    model, desired_samples=5000, burn_in=0, lag=1
)
plot(x)

@time x = mcmc_sampler(
    model,
    desired_samples=50,
    lag=500,
    burn_in=10000
)
acceptance_prob(mcmc_sampler)

plot(x)
summaryplot(x)

E_prior = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,7.0)

data = x.sample
posterior = SimPosterior(data, E_prior, γ_prior)

# Construct posterior sampler
β = 0.8
mode_move = InvMcmcMixtureMove(
    (
        EditAllocationMove(ν=1),
        PathPermutationMove(ν=1),
        InsertDeleteMove(ν=1, len_dist=TrGeometric(0.8, 1, model.K_inner.u)),
        SplitMergeMove(ν=1)
    ),
    (β/2, β/2, (1-β)/2, (1-β)/2)
    # (β, 1-β)
)

posterior_sampler = IexMcmcSampler(
    mode_move, mcmc_sampler,
    ε=0.3,
    aux_init_at_prev=true
)

# E_init = shuffle.(model.mode)
E_init = [vcat(model.mode...)]
@time x = posterior_sampler(
    posterior,
    desired_samples=1_000,
    S_init=E_init,
    γ_init=3.5
)
acceptance_prob(posterior_sampler)

# With fast matching distance
d_fast = FastMatchingDistance(FastLCS(101), 1000)
model_fast = SIM(
    E, 3.5,
    d_fast,
    1:10,
    K_inner, K_outer
)

E_prior_fast = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
posterior_fast = SimPosterior(data, E_prior_fast, γ_prior)

E_init = [vcat(model.mode...)]
@time x = posterior_sampler(
    posterior_fast,
    desired_samples=1_000,
    S_init=E_init,
    γ_init=3.5
)
acceptance_prob(posterior_sampler)

plot(x.suff_stat_trace)
plot(x,E)
x.S_sample


text:
cd(@__DIR__)

@testset "Notes" begin
    midi = load("doxy.mid")

    notes = getnotes(midi.tracks[4])

    @test notes.tpq == 960
    @test typeof(notes[1]) == Note
    @test typeof(notes[1:3]) == Notes{Note}
    @test typeof(notes[1:3]) <: Notes
    @test notes[1:3].notes == notes.notes[1:3]
end

@testset "pitch names" begin
    a = ["C5", "D12", "D#3", "G#52","A23"]
    for n in a
        @test n == replace(pitch_to_name(name_to_pitch(n)), "♯" => "#")
    end
    @test pitch_to_name(name_to_pitch("E#7")) == "F7"
    @test pitch_to_name(name_to_pitch("B#")) == "C5"
    @test pitch_to_name(name_to_pitch("F♯")) == "F♯4"
    @test name_to_pitch("C4") == 60
    @test pitch_to_name(name_to_pitch("C")) == "C4"
    @test pitch_to_name(name_to_pitch("Cb")) == "B3"
    @test pitch_to_name(name_to_pitch("G♭")) == "F♯4"
    @test pitch_to_name(name_to_pitch("G♭"); flat=true) == "G♭4"
    @test pitch_to_name(name_to_pitch("Bb7")) == "A♯7"

    n = Note(0, 1, 1, 1)
    @test pitch_to_name(n.pitch) == "C-1"
end

@testset "copying notes" begin
    midi = load("doxy.mid")
    notes = getnotes(midi.tracks[4])
    n2 = copy(notes)
    notes[1].pitch = 1
    @test n2[1].pitch ≠ 1
end


text:
# sampling functions

"""
    sample_trajectory!(x, sampler; options=MDOptions())

In place applciation of the `sampler` to `x`.  Number of iterations are
set using the `options` argument.

### Fields

* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
### Optional Fields
* `options`   - Sampling options, including number of iteration
"""
function sample_trajectory!(x::Tx, sampler::S; options = MDOptions()) where {Tx,S<:AbstractSampler}

    state = InitState!(x, sampler)
    for _ = 1:options.n_iters
        UpdateState!(state, sampler)
    end
    x
end

"""
    sample_trajectory(x₀, sampler; options=MDOptions())

Run the `sampler` starting at `x₀`.  Number of iterations and interval between
saves are set using the `options` argument.  For Metropolis samplers, the
running acceptance rates are also resturned.

### Fields
* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
### Optional Fields
* `options`   - Sampling options, including number of iteration
* `constraints` - Constraints on the trajectory
"""
function sample_trajectory(x₀::Tx, sampler::S; options = MDOptions()) where {Tx,S<:MetropolisSampler}

    n_accept = Int(0)

    state = InitState(x₀, sampler)

    # allocate memory for samples
    samples = Tx[similar(x₀) for i = 1:options.n_save]
    acceptance_rates = zeros(options.n_save)
    save_index = 1
    for i = 1:options.n_iters
        UpdateState!(state, sampler)
        n_accept += state.accept
        if (mod(i, options.n_save_iters) == 0)
            @. samples[save_index] = deepcopy(state.x)
            acceptance_rates[save_index] = n_accept / i
            save_index += 1
        end
    end
    return samples, acceptance_rates
end

function sample_trajectory(x₀::Tx, sampler::S; options = MDOptions()) where {Tx,S<:NonMetropolisSampler}

    state = InitState(x₀, sampler)
    samples = Tx[similar(x₀) for i = 1:options.n_save]
    save_index = 1
    for i = 1:options.n_iters
        UpdateState!(state, sampler)
        if (mod(i, options.n_save_iters) == 0)
            @. samples[save_index] = deepcopy(state.x)
            save_index += 1
        end
    end
    return samples
end


text:
@testset "Utils" begin
    orthants = Orthants(3)

    @test length(orthants) == 8
    eltype(Orthants) == Vector{Int}
    @test first(orthants) == [1, 1, 1]
    @test last(orthants) == [-1, -1, -1]
    @test orthants[4] == [-1, -1, 1]
end


text:
# An example of the ScalarWaveEquation code for the propagation
# of a gaussian wave in time

using ScalarWaveEquation
using Printf
using PyPlot

function get_gaussian(N::Int)
	# returns a normalized Gaussian distribution with (mean, std) = (0.5, 0.1)
	mean = 0.5
	std  = 0.1

	xdata      = range(0.0, stop = 1.0, length = N)
	exponent   = -0.5 * (xdata.-mean) .* (xdata.-mean) / (std*std)
	norm_coeff = std*sqrt(2*pi)

	return exp.(exponent)/norm_coeff
end

Ttot = 0.3  # total duration
Nt   = 1000  # number of time divisions
Nx   = 400  # number of divisions on position axis

# initial conditions are accepted in the form of FunArr instances
# we create them from 1D Arrays
# the x-range is assumed to be [0, 1]
# both initial conditions must have the same size
init_positions  = get_gaussian(Nx)
init_velocities = zeros(Nx)
ufun  = FunArr{Float64, Float64}(init_positions)
utfun = FunArr{Float64, Float64}(init_velocities)

# use the method of lines to calculate u(x, t) over [0, 1]x[0, Ttot]
# returns a FunArr2D instance with the 2D Array
# also returns the energy at each of the Nt time steps
u2fun, energies = scalarwave_and_energy(ufun, utfun, Nt, Ttot)

# the energy tends to increase over time
# a small energy increase indicates good convergence
@printf("The lowest energy calculated is %.8f\n", minimum(energies))
@printf("The highest energy calculated is %.8f\n", maximum(energies))

# u(x, t) can be plotted as a colour map
xlabel("t")
ylabel("x")
pcolormesh(u2fun.tdat, u2fun.xdat, u2fun.u2dat)
colorbar()
savefig("images/figure1.png", dpi = 300)


text:
int_rules_5_1_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.1.5*u*(a+b*arcsin(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * ~x), ~x) => Subst(Antiderivative(((~a + ~b * ~x) ^ ~n * cos(~x)) / (~c * ~d + ~e * sin(~x)), ~x), ~x, asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * ~x), ~x) => -(Subst(Antiderivative(((~a + ~b * ~x) ^ ~n * sin(~x)) / (~c * ~d + ~e * cos(~x)), ~x), ~x, acos(~c * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ ~n) / (~e * (~m + 1)) - ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ ~n) / (~e * (~m + 1)) + ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~m, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~m, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * cos(~x) * (~c * ~d + ~e * sin(~x)) ^ ~m, ~x), ~x, asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sin(~x) * (~c * ~d + ~e * cos(~x)) ^ ~m, ~x), ~x, acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide(ExpandExpression(~Px, ~x), ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide(ExpandExpression(~Px, ~x), ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && PolynomialQ(~Px, ~x)
    #= Antiderivative((~Px)*((~a')+(~b')*asin((~c')*(~x)))^(~n'),~x) :=  With([u=IntHide(Px,x)], Dist((a+b*asin(c*x))^n,u,x) -  b*c*n*Antiderivative(SimplifyIntegrand(u*(a+b*asin(c*x))^(n-1)/sqrt(1-c^2*x^2), x),x)) <-- FreeQ([a,b,c],x) && PolynomialQ(Px,x) && IGtQ(n,0) =#
    #= Antiderivative((~Px)*((~a')+(~b')*acos((~c')*(~x)))^(~n'),~x) :=  With([u=IntHide(Px,x)], Dist((a+b*acos(c*x))^n,u,x) +  b*c*n*Antiderivative(SimplifyIntegrand(u*(a+b*acos(c*x))^(n-1)/sqrt(1-c^2*x^2), x),x)) <-- FreeQ([a,b,c],x) && PolynomialQ(Px,x) && IGtQ(n,0) =#
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide(~Px * (~d + ~e * ~x) ^ ~m, ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide(~Px * (~d + ~e * ~x) ^ ~m, ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(p') * (~d + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~p * (~d + ~e * ~x) ^ ~m, ~x)], Dist((~a + ~b * asin(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && (ILtQ(~m, 0) && LtQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(p') * (~d + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~p * (~d + ~e * ~x) ^ ~m, ~x)], Dist((~a + ~b * acos(~c * ~x)) ^ ~n, u, ~x) + ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && (ILtQ(~m, 0) && LtQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative(((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n) / (~d + ~(e') * ~x) ^ 2, ~x) => With([u = IntHide((~f + ~g * ~x + ~h * (~x) ^ 2) ^ ~p / (~d + ~e * ~x) ^ 2, ~x)], Dist((~a + ~b * asin(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && EqQ(~e * ~g - 2 * ~d * ~h, 0)))
    @apply_utils Antiderivative(((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n) / (~d + ~(e') * ~x) ^ 2, ~x) => With([u = IntHide((~f + ~g * ~x + ~h * (~x) ^ 2) ^ ~p / (~d + ~e * ~x) ^ 2, ~x)], Dist((~a + ~b * acos(~c * ~x)) ^ ~n, u, ~x) + ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && EqQ(~e * ~g - 2 * ~d * ~h, 0)))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~d + ~e * ~x) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolynomialQ(~Px, ~x) && (IGtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~d + ~e * ~x) ^ ~m * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolynomialQ(~Px, ~x) && (IGtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && (LtQ(~m, -2 * ~p - 1) || GtQ(~m, 3))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(Dist(1 / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && (LtQ(~m, -2 * ~p - 1) || GtQ(~m, 3))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~d, 0) && (IGtQ(~n, 0) && (~m == 1 || (~p > 0 || (~n == 1 && ~p > -1 || ~m == 2 && ~p < -2))))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~d, 0) && (IGtQ(~n, 0) && (~m == 1 || (~p > 0 || (~n == 1 && ~p > -1 || ~m == 2 && ~p < -2))))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~m * sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~d * ~g * ~m + 2 * ~e * ~f * ~x + ~e * ~g * (~m + 2) * (~x) ^ 2) * (~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (ILtQ(~m, 0) && (GtQ(~d, 0) && IGtQ(~n, 0))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~m * sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (-((~f + ~g * ~x) ^ ~m) * (~d + ~e * (~x) ^ 2) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) + (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~d * ~g * ~m + 2 * ~e * ~f * ~x + ~e * ~g * (~m + 2) * (~x) ^ 2) * (~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (ILtQ(~m, 0) && (GtQ(~d, 0) && IGtQ(~n, 0))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IGtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * acos(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IGtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative(ExpandIntegrand((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1), (~d * ~g * ~m + ~e * ~f * (2 * ~p + 1) * ~x + ~e * ~g * (~m + 2 * ~p + 1) * (~x) ^ 2) * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (ILtQ(~m, 0) && (IGtQ(~p - 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (-((~f + ~g * ~x) ^ ~m) * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) + (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative(ExpandIntegrand((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1), (~d * ~g * ~m + ~e * ~f * (2 * ~p + 1) * ~x + ~e * ~g * (~m + 2 * ~p + 1) * (~x) ^ 2) * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (ILtQ(~m, 0) && (IGtQ(~p - 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => ((~f + ~g * ~x) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (GtQ(~d, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (-((~f + ~g * ~x) ^ ~m) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) + ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~m, 0) && (GtQ(~d, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (1 / ((~c) ^ (~m + 1) * sqrt(~d))) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * (~c * ~f + ~g * sin(~x)) ^ ~m, ~x), ~x, asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (GtQ(~d, 0) && (GtQ(~m, 0) || IGtQ(~n, 0)))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (-1 / ((~c) ^ (~m + 1) * sqrt(~d))) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * (~c * ~f + ~g * cos(~x)) ^ ~m, ~x), ~x, acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (GtQ(~d, 0) && (GtQ(~m, 0) || IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * asin(~c * ~x)) ^ ~n / sqrt(~d + ~e * (~x) ^ 2), (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acos(~c * ~x)) ^ ~n / sqrt(~d + ~e * (~x) ^ 2), (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative((~f + ~g * ~x) ^ ~m * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative((~f + ~g * ~x) ^ ~m * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0)))))
    @apply_utils Antiderivative((log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (log(~h * (~f + ~g * ~x) ^ ~m) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~a + ~b * asin(~c * ~x)) ^ (~n + 1) / (~f + ~g * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (-(log(~h * (~f + ~g * ~x) ^ ~m)) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) + ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~a + ~b * acos(~c * ~x)) ^ (~n + 1) / (~f + ~g * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative(log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(log(~h * (~f + ~g * ~x) ^ ~m) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0))))
    @apply_utils Antiderivative(log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(log(~h * (~f + ~g * ~x) ^ ~m) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * (~f + ~(g') * ~x) ^ ~m * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && ILtQ(~m + 1 / 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * (~f + ~(g') * ~x) ^ ~m * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(Dist(1 / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && ILtQ(~m + 1 / 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * asin(~c * ~x), v, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(v / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * acos(~c * ~x), v, ~x) + ~b * ~c * Antiderivative(SimplifyIntegrand(v / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~Px * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (PolynomialQ(~Px, ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2)))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (PolynomialQ(~Px, ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2)))
    @apply_utils Antiderivative(~(Px') * (~f + ~(g') * (~d + ~(e') * (~x) ^ 2) ^ ~p) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~f + ~g * (~d + ~e * (~x) ^ 2) ^ ~p) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (PolynomialQ(~Px, ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~p + 1 / 2, 0) && IntegersQ(~m, ~n))))
    @apply_utils Antiderivative(~(Px') * (~f + ~(g') * (~d + ~(e') * (~x) ^ 2) ^ ~p) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~f + ~g * (~d + ~e * (~x) ^ 2) ^ ~p) ^ ~m * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (PolynomialQ(~Px, ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~p + 1 / 2, 0) && IntegersQ(~m, ~n))))
    @apply_utils Antiderivative(~RFx * asin(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand(asin(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ(~c, ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * acos(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand(acos(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ(~c, ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~a + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~RFx * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~a + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~RFx * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * asin(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * asin(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * acos(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * acos(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~a + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p, ~RFx * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~a + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p, ~RFx * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/5.1.5 u (a+b arcsin(c x))^n.jl


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3

# Short description of this file:This script simulates the Fisher-KPP equation
#and fits a neural PDE to the data with the growth (aka reaction) term replaced
# by a feed-forward neural network and the diffusion term with a CNN.
#The resulting equation is solved with a RNN.

#Warning: Reset Julia Kernel (pre-compilation required) to avoid errors in
#redefinitions of constants.

cd(@__DIR__)
using Pkg#; Pkg.activate("."); Pkg.instantiate()
Pkg.add("BSON")
Pkg.add("DifferentialEquations")

using PyPlot, Printf
using LinearAlgebra
using Flux, DiffEqFlux, Optim, DiffEqSensitivity
using BSON: @save, @load
using Flux: @epochs
using DifferentialEquations

#parameter
D = 0.01; #diffusion
r = 1.0; #reaction rate

#domain
X = 1.0; T = 5.0;
dx = 0.04; dt = T/10;
x = collect(0:dx:X);
t = collect(0:dt:T);
tTrain  = collect(0:dt:T/2);
tValid  = collect(T/2:dt:T);
Nx = Int64(X/dx+1);
Nt = Int64(T/dt+1);

#initial conditions
Amp = 1.0;
Delta = 0.2
#IC-1
rho0 = Amp*(tanh.((x .- (0.5 - Delta/2))/(Delta/10)) - tanh.((x .- (0.5 + Delta/2))/(Delta/10)))/2
#IC-2
#rho0 = Amp*(1 .- tanh.((x .- 0.2)/(Delta/6)))/2.

save_folder = "data"

if isdir(save_folder)
    rm(save_folder, recursive=true)
end
mkdir(save_folder)

close("all")
figure()
plot(x, rho0)
title("Initial Condition")
gcf()

########################
# Generate training data
########################
reaction(u) = r * u .* (1 .- u)
lap = diagm(0 => -2.0 * ones(Nx), 1=> ones(Nx-1), -1 => ones(Nx-1)) ./ dx^2
#Periodic BC
lap[1,end] = 1.0/dx^2
lap[end,1] = 1.0/dx^2
#Neumann BC
#lap[1,2] = 2.0/dx^2
#lap[end,end-1] = 2.0/dx^2

function rc_ode(rho, p, t)
    #finite difference
    D * lap * rho + reaction.(rho)
end

prob = ODEProblem(rc_ode, rho0, (0.0, T), saveat=dt)
sol = DifferentialEquations.solve(prob, RK4()) #Tsit5()

ode_data = Array(sol);

figure(figsize=(8,3))

subplot(121)
pcolor(x,t,ode_data')
xlabel("x"); ylabel("t");
colorbar()

subplot(122)
for i in 1:2:Nt
    plot(x, ode_data[:,i], label="t=$(sol.t[i])")
end
xlabel("x"); ylabel(L"$\rho$")
legend(frameon=false, fontsize=7, bbox_to_anchor=(1, 1), loc="upper left", ncol=1)
tight_layout()
savefig(@sprintf("%s/training_data.pdf", save_folder))
gcf()

########################
# Define the neural PDE
########################
n_weights = 10

#for the reaction term
rx_nn = Chain(Dense(1, n_weights, tanh),
                Dense(n_weights, 2*n_weights, tanh),
                Dense(2*n_weights, n_weights, tanh),
                Dense(n_weights, 1),
                x -> x[1])

#conv with bias with initial values as 1/dx^2
w_err = 0.0
init_w = reshape([1.1 -2.5 1.0], (3, 1, 1, 1))
diff_cnn_ = Conv(init_w, [0.], pad=(0,0,0,0))

#initialize D0 close to D/dx^2
D0 = [6.5]

p1,re1 = Flux.destructure(rx_nn)
p2,re2 = Flux.destructure(diff_cnn_)
p = [p1;p2;D0]
full_restructure(p) = re1(p[1:length(p1)]), re2(p[(length(p1)+1):end-1]), p[end]

function nn_ode(u,p,t)
    rx_nn = re1(p[1:length(p1)])

    u_cnn_1   = [p[end-4] * u[end] + p[end-3] * u[1] + p[end-2] * u[2]]
    u_cnn     = [p[end-4] * u[i-1] + p[end-3] * u[i] + p[end-2] * u[i+1] for i in 2:Nx-1]
    u_cnn_end = [p[end-4] * u[end-1] + p[end-3] * u[end] + p[end-2] * u[1]]

    # Equivalent using Flux, but slower!
    #CNN term with periodic BC
    #diff_cnn_ = Conv(reshape(p[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
    #u_cnn = reshape(diff_cnn_(reshape(u, (Nx, 1, 1, 1))), (Nx-2,))
    #u_cnn_1 = reshape(diff_cnn_(reshape(vcat(u[end:end], u[1:1], u[2:2]), (3, 1, 1, 1))), (1,))
    #u_cnn_end = reshape(diff_cnn_(reshape(vcat(u[end-1:end-1], u[end:end], u[1:1]), (3, 1, 1, 1))), (1,))

    [rx_nn([u[i]])[1] for i in 1:Nx] + p[end] * vcat(u_cnn_1, u_cnn, u_cnn_end)
end

########################
# Soving the neural PDE and setting up loss function
########################
prob_nn = ODEProblem(nn_ode, rho0, (0.0, T), p)
sol_nn = concrete_solve(prob_nn,Tsit5(), rho0, p)

function predict_rd(θ)
  # No ReverseDiff if using Flux
  Array(concrete_solve(prob_nn,Tsit5(),rho0,θ,saveat=dt,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

#match data and force the weights of the CNN to add up to zero
function loss_rd(p)
    pred = predict_rd(p)
    sum(abs2, ode_data .- pred) + 10^2 * abs(sum(p[end-4 : end-2])), pred
end

########################
# Training
########################
#Optimizer
opt = ADAM(0.001)

global count = 0
global save_count = 0
save_freq = 50

train_arr = Float64[]
diff_arr = Float64[]
w1_arr = Float64[]
w2_arr = Float64[]
w3_arr = Float64[]

#callback function to observe training
cb = function (p,l,pred)
    rx_nn, diff_cnn_, D0 = full_restructure(p)
    push!(train_arr, l)
    push!(diff_arr, p[end])

    weight = diff_cnn_.weight[:]
    push!(w1_arr, weight[1])
    push!(w2_arr, weight[2])
    push!(w3_arr, weight[3])

    println(@sprintf("Loss: %0.4f\tD0: %0.4f Weights:(%0.4f,\t %0.4f, \t%0.4f) \t Sum: %0.4f"
            ,l, D0[1], weight[1], weight[2], weight[3], sum(weight)))

    global count

    if count==0
        fig = figure(figsize=(8,2.5));
        ttl = fig.suptitle(@sprintf("Epoch = %d", count), y=1.05)
        global ttl
        subplot(131)
        pcolormesh(x,t,ode_data')
        xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
        colorbar()

        subplot(132)
        img = pcolormesh(x,t,pred')
        global img
        xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
        colorbar(); clim([0, 1]);

        ax = subplot(133); global ax
        u = collect(0:0.01:1)
        rx_line = plot(u, rx_nn.([[elem] for elem in u]), label="NN")[1];
        global rx_line
        plot(u, reaction.(u), label="True")
        title("Reaction Term")
        legend(loc="upper right", frameon=false, fontsize=8);
        ylim([0, r*0.25+0.2])

        subplots_adjust(top=0.8)
        tight_layout()
    end

    if count>0
        println("updating figure")
        img.set_array(pred[1:end-1, 1:end-1][:])
        ttl.set_text(@sprintf("Epoch = %d", count))

        u = collect(0:0.01:1)
        rx_pred = rx_nn.([[elem] for elem in u])
        rx_line.set_ydata(rx_pred)
        u = collect(0:0.01:1)

        min_lim = min(minimum(rx_pred), minimum(reaction.(u)))-0.1
        max_lim = max(maximum(rx_pred), maximum(reaction.(u)))+0.1

        ax.set_ylim([min_lim, max_lim])
    end

    global save_count
    if count%save_freq == 0
        println("saved figure")
        savefig(@sprintf("%s/pred_%05d.png", save_folder, save_count), dpi=200, bbox_inches="tight")
        save_count += 1
    end

    display(gcf())
    count += 1

    false
end

#train
res1 = DiffEqFlux.sciml_train(loss_rd, p, ADAM(0.01), cb=cb, maxiters = 10)
res2 = DiffEqFlux.sciml_train(loss_rd, res1.minimizer, ADAM(0.01), cb=cb, maxiters = 50)
res3 = DiffEqFlux.sciml_train(loss_rd, res2.minimizer, BFGS(), cb=cb, maxiters = 100)

pstar = res3.minimizer

### Recurrent Neural Network

master_rnn = Flux.RNN(26, 1, (x -> x))

function predict_rnn(x)
  out = master_rnn.(x)[end]
  Flux.reset!(master_rnn)
  out
end

function loss_rnn(p)
    #pred = predict_rd(p)
    pred_rnn = predict_rnn(p)
    sum(abs2, ode_data .- pred_rnn), pred_rnn
end

ps = Flux.params(master_rnn)

train_data = [nn_ode(predict_rd(p)[:,k], p, tTrain) for k in 1:length(tTrain)]
test_data = [nn_ode(predict_rd(p)[:,k], p, tTrain) for k in length(tTrain):length(tTrain)+ length(tValid) -1]

train_labels = [predict_rd(p)[:,k] for k in 1:length(tTrain)]
test_labels = [predict_rd(p)[:,k] for k in length(tTrain):length(tTrain) + length(tValid) - 1]

loss(x, y) = abs(sum((predict_rnn(x) .- y)))

n_epochs = 50

@epochs n_epochs Flux.train!(loss, ps, zip(train_data, train_labels), opt)

using BenchmarkTools
@btime pred_test = [predict_rnn(test_data[k]) for k in 1:length(test_data)]

train_lab = reshape(train_labels[1], : , 1)
for k in 2:6
    train_lab = cat(train_lab, reshape(train_labels[k], : , 1), dims= (2))

end

fig = figure(figsize=(8,2.5));
subplot(131)
pcolormesh(x,tTrain,train_lab')
xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
colorbar()

pred_test_fig = reshape(pred_test[1], : , 1)
for k in 2:6
    pred_test_fig = cat(pred_test_fig, reshape(pred_test[k], : , 1), dims= (2))

end

subplot(132)
img = pcolormesh(x,tValid,pred_test_fig')
global img
xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
colorbar(); clim([0, 1]);

display(gcf())

## Save trained model
"""
function pipeline()
    u = u0
    for i in t
        uderiv = nn_ode(u, p, t)
        u = predict_rnn(zip(u, uderiv))
    end
end
"""


text:
module TestVisualizationData

using Test
using Gridap.Geometry
using Gridap.CellData
using Gridap.FESpaces
using Gridap.TensorValues
using Gridap.Algebra
using Gridap.Visualization
using Gridap.ReferenceFEs
using Gridap.Geometry: get_reffes

domain = (0,2pi, 0, 4pi)
partition = (2, 3)
model = simplexify(CartesianDiscreteModel(domain,partition))
f(pt) = sin(pt[1])
V0 = TestFESpace(
  reffe=:Lagrangian, order=1, valuetype=Float64,
  conformity=:H1, model=model)
U = TrialFESpace(V0)
trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)
A(u, v) = u ⊙ v
b(v) = v⊙f
t_Ω = AffineFETerm(A,b,trian,quad)
op = AffineFEOperator(U,V0,t_Ω)
u = solve(op)

visdata = visualization_data(trian, cellfields=Dict("u" =>u))

ncells = num_cells(visdata.grid)
nodes_per_cell = num_nodes(first(get_reffes(visdata.grid)))
@test size(visdata.nodaldata["u"]) == (ncells * nodes_per_cell,)
@test size(visdata.celldata["cell"]) == (ncells,)
@test visdata.grid isa Visualization.VisualizationGrid
@test visdata isa Visualization.VisualizationData

end#module


text:
# ---
# title: 1029. Two City Scheduling
# id: problem1029
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/two-city-scheduling/description/>
# hidden: true
# ---
#
# A company is planning to interview `2n` people. Given the array `costs` where
# `costs[i] = [aCosti, bCosti]`, the cost of flying the `ith` person to city `a`
# is `aCosti`, and the cost of flying the `ith` person to city `b` is `bCosti`.
#
# Return _the minimum cost to fly every person to a city_ such that exactly `n`
# people arrive in each city.
#
#
#
# **Example 1:**
#
#
#
#     Input: costs = [[10,20],[30,200],[400,50],[30,20]]
#     Output: 110
#     Explanation:
#     The first person goes to city A for a cost of 10.
#     The second person goes to city A for a cost of 30.
#     The third person goes to city B for a cost of 50.
#     The fourth person goes to city B for a cost of 20.
#
#     The total minimum cost is 10 + 30 + 50 + 20 = 110 to have half the people interviewing in each city.
#
#
# **Example 2:**
#
#
#
#     Input: costs = [[259,770],[448,54],[926,667],[184,139],[840,118],[577,469]]
#     Output: 1859
#
#
# **Example 3:**
#
#
#
#     Input: costs = [[515,563],[451,713],[537,709],[343,819],[855,779],[457,60],[650,359],[631,42]]
#     Output: 3086
#
#
#
#
# **Constraints:**
#
#   * `2 * n == costs.length`
#   * `2 <= costs.length <= 100`
#   * `costs.length` is even.
#   * `1 <= aCosti, bCosti <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CSTUtils
    export OverlayNode
    export function_def_call, isexpr, text, trailing_trivia, leading_trivia,
        span_text, fullspan_text, id_name

    include("treewalking.jl")

    function text(o::OverlayNode, leading_trivia = true, trailing_trivia = true)
        o.buffer[((leading_trivia ? first(o.fullspan) : first(o.span)):(trailing_trivia ? last(o.fullspan) : last(o.span)))]
    end

    function trailing_trivia(o::OverlayNode)
        o.buffer[(1+last(o.span):last(o.fullspan))]
    end

    function leading_trivia(o::OverlayNode)
        o.buffer[(first(o.fullspan):first(o.fullspan)-1)]
    end

    span_text(o::OverlayNode) = text(o, false, false)
    fullspan_text(o::OverlayNode) = text(o, true, true)

    isexpr(x::EXPR{T}, ::Type{S}) where {T, S} = T == S

    isexpr(x::T, ::Type{T}) where {T} = true
    isexpr(x::T1, ::Type{T2}) where {T1, T2} = false
    # Operator
    isexpr(x::CSTParser.OPERATOR, o::Type{CSTParser.OPERATOR}, op::Tokens.Kind) = x.kind == op && x.dot == false
    isexpr(x, o::Type{CSTParser.OPERATOR}, op::Tokens.Kind) = false
    # Keyword
    isexpr(x::CSTParser.KEYWORD, k::Type{CSTParser.KEYWORD}, kw::Tokens.Kind) = x.kind == kw
    isexpr(x, k::Type{CSTParser.KEYWORD}, kw::Tokens.Kind) = false
    # Literal
    isexpr(x::CSTParser.LITERAL, k::Type{CSTParser.LITERAL}, lit::Tokens.Kind) = x.kind == lit
    isexpr(x, k::Type{CSTParser.LITERAL}, lit::Tokens.Kind) = false

    isexpr(x::CSTParser.PUNCTUATION, k::Type{CSTParser.PUNCTUATION}, lit::Tokens.Kind) = x.kind == lit
    isexpr(x, k::Type{CSTParser.PUNCTUATION}, lit::Tokens.Kind) = false

    isexpr(x::OverlayNode, ::Type{S}, kind::Tokens.Kind) where {S} = isexpr(x.expr, S, kind)
    isexpr(x::OverlayNode, ::Type{CSTParser.OPERATOR}, kind::Tokens.Kind)  = isexpr(x.expr, CSTParser.OPERATOR, kind)
    isexpr(x::OverlayNode, ::Type{CSTParser.KEYWORD}, kind::Tokens.Kind)  = isexpr(x.expr, CSTParser.KEYWORD, kind)
    isexpr(x::OverlayNode, ::Type{CSTParser.LITERAL}, kind::Tokens.Kind) = isexpr(x.expr, CSTParser.LITERAL, kind)
    isexpr(x::OverlayNode, ::Type{CSTParser.PUNCTUATION}, kind::Tokens.Kind) = isexpr(x.expr, CSTParser.PUNCTUATION, kind)
    isexpr(x::OverlayNode, ::Type{S}) where {S} = isexpr(x.expr, S)

    id_name(x::CSTParser.IDENTIFIER) = x.val
    id_name(x::OverlayNode) = id_name(x.expr)
end