text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SDL2_image"
version = v"2.0.5"

# Collection of sources required to complete build
sources = [
    "http://www.libsdl.org/projects/SDL_image/release/SDL2_image-2.0.5.zip" =>
    "eee0927d1e7819d57c623fe3e2b3c6761c77c474fe9bc425e8674d30ac049b1c",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SDL2_image-*/
export CPPFLAGS="-I${prefix}/include"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-pic
make -j${nproc}
make install
if [[ "${target}" == *-freebsd* ]]; then
    # We need to manually build the shared library for FreeBSD
    cd "${libdir}"
    ar x libSDL2_image.a
    cc -shared -o libSDL2_image.${dlext} *.o
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libSDL2_image", "SDL2_image"], :libsdl2_image)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    PackageSpec(name="SDL2_jll", uuid="ab825dc5-c88e-5901-9575-1e5e20358fcf")
    PackageSpec(name="libwebp_jll", uuid="c5f90fcd-3b7e-5836-afba-fc50a0988cb2")
    PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8")
    PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f")
    PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828")
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""

# basis_set

Determine basis_set

$(SIGNATURES)

Part of API, exported.
"""
function basis_set(dag::DAG; debug=false)
  as = transpose(topological_sort(dag.a))
  nod = dag.vars[topological_order(dag.a)]
  dv = 1:length(nod)

  # Correct if adj matixi is upper triangular
  @assert istriu(as) "Sorted adj matrix not upper triangular."

  ind = Vector{Vector{Symbol}}()
  for r in dv
    for s in r:length(nod)
      ((as[r, s] !== 0) || (s == r)) && continue
      ed = nod[[r, s]]
      pa_r = filter(i -> as[i, r] == 1, dv)
      pa_s = filter(i -> as[i, s] == 1, dv)
      dsep = union(nod[pa_r], nod[pa_s])
      debug && println("r = $(nod[r]), s= $(nod[s]), pa_r = $(nod[pa_r]), pa_s = $(nod[pa_s]), dsep = $dsep")
      dsep = setdiff(dsep, ed)
      append!(ed, dsep)
      append!(ind, [ed])
    end
  end
  BasisSet(dag, ind)
end

export
  basis_set

ext: jl


text:
import Gaugefields:staggered_U

struct StaggeredFermion_2D_wing{NC} <: AbstractFermionfields_2D{NC}
    NC::Int64
    NX::Int64
    #NY::Int64
    #NZ::Int64
    NT::Int64
    NDW::Int64
    NG::Int64 #size of the Gamma matrix. In Staggered fermion, this is one.
    NV::Int64
    f::Array{ComplexF64,4}
    Dirac_operator::String

    function StaggeredFermion_2D_wing(NC,NX,NT)
        NG = 1
        NDW = 1
        NY = 1
        NZ = 1
        NV = NC*NX*NY*NZ*NT*NG
        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NX+2NDW,NT+2NDW,NG)
        Dirac_operator = "Staggered"
        return new{NC}(NC,NX,NT,NDW,NG,NV,f,Dirac_operator)
    end
end

function Base.size(x::StaggeredFermion_2D_wing{NC}) where NC
    return (x.NC,x.NX,x.NT,x.NG)
    #return (x.NV,)
end

function Base.length(x::StaggeredFermion_2D_wing{NC}) where {NC}
    return NC*x.NX*x.NT*x.NG
end

function Base.similar(x::T) where T <: StaggeredFermion_2D_wing
    return StaggeredFermion_2D_wing(x.NC,x.NX,x.NT)
end

function Dx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1},boundarycondition) where  {T <: StaggeredFermion_2D_wing,G <:AbstractGaugefields}
    #temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    #clear!(temp)
    set_wing_fermion!(x,boundarycondition)
    clear_fermion!(xout)
    for ν=1:2
        xplus = shift_fermion(x,ν)
        Us = staggered_U(U[ν],ν)
        mul!(temp1,Us,xplus)

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        Uminus_s = staggered_U(Uminus,ν)
        mul!(temp2,Uminus_s',xminus)

        add_fermion!(xout,0.5,temp1,-0.5,temp2)

        #fermion_shift!(temp1,U,ν,x)
        #fermion_shift!(temp2,U,-ν,x)
        #add!(xout,0.5,temp1,-0.5,temp2)

    end

    set_wing_fermion!(xout,boundarycondition)

    return
end

function clear_fermion!(x::StaggeredFermion_2D_wing{NC},evensite) where NC
    ibush = ifelse(evensite,0,1)
    for it=1:x.NT
        #for iz=1:x.NZ
            #for iy=1:x.NY
                xran =1+(1+ibush+it)%2:2:x.NX
                for ix in xran
                    @simd for ic=1:NC
                        x[ic,ix, it,1] = 0
                    end
                end
            #end
        #end
    end
    return
end

function set_wing_fermion!(a::StaggeredFermion_2D_wing{NC},boundarycondition) where NC
    NT = a.NT
    #NZ = a.NZ
    #NY = a.NY
    NX = a.NX

    #!  X-direction
    for ialpha=1:1
        for it=1:NT
            #for iz = 1:NZ
                #for iy=1:NY
                    @simd for k=1:NC
                        a[k,0,it,ialpha] = boundarycondition[1]*a[k,NX,it,ialpha]
                    end
                #end
            #end
        end
    end

    for ialpha=1:1
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for k=1:NC
                        a[k,NX+1,it,ialpha] = boundarycondition[1]*a[k,1, it,ialpha]
                    end
                #end
            #end
        end
    end

    for ialpha=1:1

        #T-direction
        #for iz=1:NZ
            #for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix, 0,ialpha] = boundarycondition[2]*a[k,ix, NT,ialpha]
                        a[k,ix, NT+1,ialpha] = boundarycondition[2]*a[k,ix, 1,ialpha]
                    end
                end
            #end
        #end

    end

end


text:
struct MOSResult
    structure::MOSStructure
    result::Semiconductor1D
    VGB::Number
    ψs::Number
end

tox(MOS::MOSResult) = tox(MOS.structure)

C(MOS::MOSStructure) = C(MOS.semiconductor)
tox(MOS::MOSStructure) = MOS.tox
Temperature(MOS::MOSStructure) = Temperature(MOS.semiconductor)
VFB(MOS::MOSStructure) = φₘ(MOS.metal,Temperature(MOS)) - χ(MOS.semiconductor) - Ec(MOS.semiconductor)
ϕb(MOS::MOSStructure) = abs(Ef(MOS.semiconductor) - Ei(MOS.semiconductor))
Cox(MOS::MOSStructure) = MOS.oxide.ϵr*ϵ₀/tox(MOS)
n(MOS::MOSStructure,ψ=0.0) = n(MOS.semiconductor,ψ)
p(MOS::MOSStructure,ψ=0.0) = p(MOS.semiconductor,ψ)
ξ(MOS::MOSStructure,ψ=0.0) = ξ(MOS.semiconductor,ψ)
ξ(MOS::MOSStructure,ψ,ϕ) = ξ(MOS.semiconductor,ψ,ϕ)
Ei(MOS::MOSStructure) = Ei(MOS.semiconductor)
Ec(MOS::MOSStructure) = Ec(MOS.semiconductor)
Ev(MOS::MOSStructure) = Ev(MOS.semiconductor)
Nv(MOS::MOSStructure) = Nv(MOS.semiconductor)
Nc(MOS::MOSStructure) = Nc(MOS.semiconductor)
Ef(MOS::MOSStructure) = Ef(MOS.semiconductor)
ni(MOS::MOSStructure) = ni(MOS.semiconductor)
#Temperature(MOS::MOSStructure) = Temperature(MOS.semiconductor)
Rₗ(MOS::MOSStructure) = MOS.tsi
γ(MOS::MOSStructure) = γ(abs(C(MOS)),MOS.tox)
DistFun(MOS::MOSStructure,E) = DistFun(MOS.semiconductor,E)

function ψs₀(Vg,Na,tox,Vfb,T)
    Cox = ϵ₀*3.9/tox
    ϕt = kb*T
    Gamma = Cox^-1*sqrt(2*q*11.7*ϵ₀*Na)
    α = -Gamma*sqrt(ϕt)
    ϕf = ϕt*log(Na/ni(T,BoltzmanDist(),Silicon()))
    ψs=0
    ug = Vg-Vfb
    if ug == 0.0
        return 0.0
    end
    uc = 2*ϕf
    ec = exp(-uc/ϕt)
    K = 1 - ec
    if uc > 0
        ϕsat = ug+Gamma^2/2.0*(1- sqrt(max(0,1+4*(ug-ϕt)/(Gamma^2))))
        if ϕsat >= uc && ug > 0
            ψs = uc + ϕt*log(max(1,(ug-uc)^2/Gamma^2/ϕt))
        else
            if ug < α
                ψs = -2*ϕt*log(ug/α)
            elseif ug > 0 && abs(ug) > ϕt
                ψs = ϕsat
            end
        end
    end
    ψs = ug > 0 ? min(ug,ψs) : max(ug,ψs)
    return ψs
end

function Cap(Vg,MOS::MOSStructure)
    ψ = ψs(Vg,MOS)
    Xs = 1.0/Cs(ψ,MOS)
    return 1/(1/Cox(MOS)+Xs)
end

function Cs(ψ,MOS::MOSStructure)
    if ψ ∈ -0.03..0.03
        x0 = -0.03
        x1 = 0.03
        C0 = abs(Zygote.gradient(x->-ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,x),x0)[1])
        C1 = abs(Zygote.gradient(x->-ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,x),x1)[1])
        α = (ψ-x0)/(x1-x0)
        return C0*(C1/C0)^α    #exp interp
    end
    return abs(Zygote.gradient(x->-ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,x),ψ)[1])
end

function Qs(Vg,MOS::MOSStructure)
    return -ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,ψs(Vg,MOS))
end

ψs_PSP(Vg,MOS::MOSStructure) = PSP_VaricapPsi(Vg,abs(C(MOS)),MOS.tox,VFB(MOS),Temperature(MOS),ϕb(MOS))

ψs_PSPs(Vg,V,MOS::MOSStructure) = PSP_psiS(Vg,abs(C(MOS)),MOS.tox,VFB(MOS),Temperature(MOS),ϕb(MOS),V)[1]

function ψs(Vg,MOS::MOSStructure)
    ψ₀, = PSP_VaricapPsi(Vg,abs(C(MOS)),MOS.tox,VFB(MOS),Temperature(MOS),ϕb(MOS))
    QF(psi) = -ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,psi)
    fpsiB(psi) = psi + VFB(MOS) - QF(psi)/Cox(MOS)
    return nlsolve(x->fpsiB(x[1])-Vg,[ψ₀]).zero[1]
end
function ψs(Vg,ϕ,MOS::MOSStructure)
    ψ₀, = ψs_PSPs(Vg,ϕ,MOS)
    QF(psi) = -ϵ₀*ϵᵣ(MOS.semiconductor)*ξ(MOS.semiconductor,psi,ϕ)
    fpsiB(psi) = psi + VFB(MOS) - QF(psi)/Cox(MOS)
    return nlsolve(x->fpsiB(x[1])-Vg,[ψ₀]).zero[1]
end

function Poisson(Vg,MOS::MOSStructure,N=501)
    #println("New Version Build PN")
    Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(MOS.semiconductor,Rₗ(MOS))
    Ld = Ld0/lscal
    if Ld > 1.0
        println("Silicon Too small")
    end
    F(x) = DistFun(MOS,x)
    @variables x ψ(..)
    D = Differential(x)^2
    C0 = C(MOS)/Cscal
    ul = (Ef(MOS)-Ei(MOS)+Vg)/Vscal
    eqs = [λ*D(ψ(x)) ~ δ*F(ψ(x))-δ*F(-ψ(x))-C0]
    bcs = [ψ(0) ~ ψs(Vg,MOS)/Vscal+ul,
            ψ(Ld) ~ ul]
    domain = [x ∈ (0.0,Ld)]
    pde_system = PDESystem(eqs,bcs,domain,[x],[ψ(x)])
    discretization = MOLFiniteDifference([x=>Ld/(N-1)],nothing;centered_order=2)
    prob = discretize(pde_system,discretization)
    sol = solve(prob)*Vscal .+ Ei(MOS)
    x = collect(LinRange(0,Ld0,length(sol)))
    push!(x,Rₗ(MOS))
    push!(sol,0)
    return MOSResult(MOS,Semiconductor1D(MOS.semiconductor,x,sol),Vg,sol[1])
end

function PoissonVoronoi(Vg,MOS::MOSStructure,N=20)
    Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(MOS.semiconductor,Rₗ(MOS))
    x,sol = PoissonVoronoi(Ec(MOS),Ev(MOS),Nc(MOS),Nv(MOS),0.0,ψs(Vg,MOS),0.0,lscal,Ld0,MOS.tsi,λ,N,N,x->C(MOS),Vscal,Cscal)
    return MOSResult(MOS,Semiconductor1D(MOS.semiconductor,x,sol),Vg,sol[1])
end

function BandDiagram(MOS::MOSResult)
    tsi = MOS.structure.tsi
    #print(tsi)
    BDplot = BandDiagramJ(MOS.result,tsi)
    BDplot = plot!(BDplot,[-tox(MOS)*1e4,-tox(MOS)*1e4],[-10,10],label=nothing,c=:black)
    BDplot = plot!(BDplot,[0,0],[-10,10],label=nothing,c=:black,ylims=(-2.5,2.5))
    BDplot = plot!(BDplot,[tox(MOS)*1e4-0.1,-tox(MOS)*1e4],[-MOS.VGB+VFB(MOS.structure),-MOS.VGB+VFB(MOS.structure)],label=nothing,c=:black,xlims=(-2.0*tox(MOS)*1e4,tsi*1e4))
    BDplot = plot!(BDplot,[-tox(MOS)*1e4,0],[-MOS.VGB+VFB(MOS.structure),Ef(MOS.result[1])-MOS.ψs],label=L"\psi_{ox}",c=:red,xlabel=L"y \ [\mu m]",legend=:outertopright)
    return BDplot
end


text:
# Autogenerated wrapper script for TempestModel_jll for x86_64-linux-gnu-cxx11
export BaldaufGravityWaveTest_exe, Baroclinic3DCartesianRidgeTest_exe, Baroclinic3DCartesianTest_exe, BaroclinicWaveJWTest_exe, BaroclinicWaveUMJSTest_exe, BarotropicInstabilityTest_exe, DensityCurrentCartesianTest_exe, HeldSuarezTest_exe, HydrostaticMountainCartesianTest_exe, InertiaGravityWaveTest_exe, InertialGravityCartesianXZTest_exe, MountainRossby3DTest_exe, MountainRossbyTest_exe, MountainWaveSphereTest_exe, NonHydroMountainCartesianTest_exe, RobertBubbleCartesianTest_exe, RossbyHaurwitzWaveTest_exe, ScharMountainCartesianTest_exe, ScharMountainSphereTest_exe, ShearJetMtnWave2DCartesianTest_exe, StationaryMountainFlowTest_exe, ThermalBubbleCartesian3DTest_exe, ThermalBubbleCartesianTest_exe

using MKL_jll
using MPICH_jll
using NetCDF_jll
using HDF5_jll
JLLWrappers.@generate_wrapper_header("TempestModel")
JLLWrappers.@declare_executable_product(BaldaufGravityWaveTest_exe)
JLLWrappers.@declare_executable_product(Baroclinic3DCartesianRidgeTest_exe)
JLLWrappers.@declare_executable_product(Baroclinic3DCartesianTest_exe)
JLLWrappers.@declare_executable_product(BaroclinicWaveJWTest_exe)
JLLWrappers.@declare_executable_product(BaroclinicWaveUMJSTest_exe)
JLLWrappers.@declare_executable_product(BarotropicInstabilityTest_exe)
JLLWrappers.@declare_executable_product(DensityCurrentCartesianTest_exe)
JLLWrappers.@declare_executable_product(HeldSuarezTest_exe)
JLLWrappers.@declare_executable_product(HydrostaticMountainCartesianTest_exe)
JLLWrappers.@declare_executable_product(InertiaGravityWaveTest_exe)
JLLWrappers.@declare_executable_product(InertialGravityCartesianXZTest_exe)
JLLWrappers.@declare_executable_product(MountainRossby3DTest_exe)
JLLWrappers.@declare_executable_product(MountainRossbyTest_exe)
JLLWrappers.@declare_executable_product(MountainWaveSphereTest_exe)
JLLWrappers.@declare_executable_product(NonHydroMountainCartesianTest_exe)
JLLWrappers.@declare_executable_product(RobertBubbleCartesianTest_exe)
JLLWrappers.@declare_executable_product(RossbyHaurwitzWaveTest_exe)
JLLWrappers.@declare_executable_product(ScharMountainCartesianTest_exe)
JLLWrappers.@declare_executable_product(ScharMountainSphereTest_exe)
JLLWrappers.@declare_executable_product(ShearJetMtnWave2DCartesianTest_exe)
JLLWrappers.@declare_executable_product(StationaryMountainFlowTest_exe)
JLLWrappers.@declare_executable_product(ThermalBubbleCartesian3DTest_exe)
JLLWrappers.@declare_executable_product(ThermalBubbleCartesianTest_exe)
function __init__()
    JLLWrappers.@generate_init_header(MKL_jll, MPICH_jll, NetCDF_jll, HDF5_jll)
    JLLWrappers.@init_executable_product(
        BaldaufGravityWaveTest_exe,
        "bin/BaldaufGravityWaveTest",
    )

    JLLWrappers.@init_executable_product(
        Baroclinic3DCartesianRidgeTest_exe,
        "bin/Baroclinic3DCartesianRidgeTest",
    )

    JLLWrappers.@init_executable_product(
        Baroclinic3DCartesianTest_exe,
        "bin/Baroclinic3DCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        BaroclinicWaveJWTest_exe,
        "bin/BaroclinicWaveJWTest",
    )

    JLLWrappers.@init_executable_product(
        BaroclinicWaveUMJSTest_exe,
        "bin/BaroclinicWaveUMJSTest",
    )

    JLLWrappers.@init_executable_product(
        BarotropicInstabilityTest_exe,
        "bin/BarotropicInstabilityTest",
    )

    JLLWrappers.@init_executable_product(
        DensityCurrentCartesianTest_exe,
        "bin/DensityCurrentCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        HeldSuarezTest_exe,
        "bin/HeldSuarezTest",
    )

    JLLWrappers.@init_executable_product(
        HydrostaticMountainCartesianTest_exe,
        "bin/HydrostaticMountainCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        InertiaGravityWaveTest_exe,
        "bin/InertiaGravityWaveTest",
    )

    JLLWrappers.@init_executable_product(
        InertialGravityCartesianXZTest_exe,
        "bin/InertialGravityCartesianXZTest",
    )

    JLLWrappers.@init_executable_product(
        MountainRossby3DTest_exe,
        "bin/MountainRossby3DTest",
    )

    JLLWrappers.@init_executable_product(
        MountainRossbyTest_exe,
        "bin/MountainRossbyTest",
    )

    JLLWrappers.@init_executable_product(
        MountainWaveSphereTest_exe,
        "bin/MountainWaveSphereTest",
    )

    JLLWrappers.@init_executable_product(
        NonHydroMountainCartesianTest_exe,
        "bin/NonHydroMountainCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        RobertBubbleCartesianTest_exe,
        "bin/RobertBubbleCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        RossbyHaurwitzWaveTest_exe,
        "bin/RossbyHaurwitzWaveTest",
    )

    JLLWrappers.@init_executable_product(
        ScharMountainCartesianTest_exe,
        "bin/ScharMountainCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        ScharMountainSphereTest_exe,
        "bin/ScharMountainSphereTest",
    )

    JLLWrappers.@init_executable_product(
        ShearJetMtnWave2DCartesianTest_exe,
        "bin/ShearJetMtnWave2DCartesianTest",
    )

    JLLWrappers.@init_executable_product(
        StationaryMountainFlowTest_exe,
        "bin/StationaryMountainFlowTest",
    )

    JLLWrappers.@init_executable_product(
        ThermalBubbleCartesian3DTest_exe,
        "bin/ThermalBubbleCartesian3DTest",
    )

    JLLWrappers.@init_executable_product(
        ThermalBubbleCartesianTest_exe,
        "bin/ThermalBubbleCartesianTest",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, XLSX

makedocs(
    format = :html,
    sitename = "XLSX.jl",
    modules = [ XLSX ],
    pages = [ "index.md", "api.md" ]
)

deploydocs(
    repo = "github.com/felipenoris/XLSX.jl.git",
    target = "build",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing
)


text:
using AdvancedTopics
using Test

@testset "AdvancedTopics.jl" begin
    # Write your own tests here.
end


text:
"""
    activationoperations(layer)::Int

Outputs the approximate mumber of operations for a Flux activation function.

```julia
layer = Conv(weight = weight,
    σ = sigmoid)

activationoperations(Conv(weight = weight,
    σ = sigmoid))
```
"""
function activationoperations(layer::Recur)
    layer = recurrentunpack(layer)
    activation = layer.σ
    if activation == sigmoid || activation == rrelu || activation == elu || activation == celu ||
        activation == softsign || activation == softplus || activation == tanh
        noofoper = 3
    elseif activation == gelu
        noofoper = 8
    elseif activation == hardsigmoid || activation == logsigmoid || activation == swish ||
            activation == selu || activation == softmax || activation == logsoftmax
        noofoper = 4
    elseif activation == hardtanh || activation == leakyrelu || activation == relu6 ||
            activation == lisht || activation == tanhshrink
        noofoper = 2
    elseif activation == logcosh || activation == mish
        noofoper = 5
    elseif activation == relu || activation == trelu || activation == softshrink
        noofoper = 1
    elseif activation == identity
        noofoper = 0
    else
        @info "Unkown activation type defaulting to identity"
        return 0
    end
    return noofoper
end

function activationoperations(layer)
    activation = layer.σ
    if activation == sigmoid || activation == rrelu || activation == elu || activation == celu ||
        activation == softsign || activation == softplus || activation == tanh
        noofoper = 4
    elseif activation == gelu
        noofoper = 9
    elseif activation == hardsigmoid || activation == logsigmoid || activation == swish ||
            activation == selu || activation == softmax || activation == logsoftmax
        noofoper = 5
    elseif activation == hardtanh || activation == leakyrelu || activation == relu6 ||
            activation == lisht || activation == tanhshrink
        noofoper = 3
    elseif activation == logcosh || activation == mish
        noofoper = 6
    elseif activation == relu || activation == trelu || activation == softshrink
        noofoper = 2
    elseif activation == identity
        noofoper = 1
    else
        @info "Unkown activation type defaulting to identity"
        return 1
    end
    return noofoper
end

"""
    gradientoperations(layer)::Int

Outputs the approximate mumber of operations for the gradient of a Flux activation
function.

```julia
layer = Conv(weight = weight,
    σ = sigmoid)

gradientoperations(Conv(weight = weight,
    σ = sigmoid))
```
"""
function gradientoperations(layer::Recur)
    layer = recurrentunpack(layer)
    activation = layer.σ
    if activation == sigmoid || activation == rrelu || activation == elu || activation == celu ||
        activation == softsign || activation == softplus || activation == tanh
        noofoper = 4
    elseif activation == gelu
        noofoper = 8
    elseif activation == hardsigmoid || activation == logsigmoid || activation == swish ||
            activation == selu || activation == softmax || activation == logsoftmax
        noofoper = 4
    elseif activation == hardtanh || activation == leakyrelu || activation == relu6 ||
            activation == lisht || activation == tanhshrink
        noofoper = 2
    elseif activation == logcosh || activation == mish
        noofoper = 5
    elseif activation == relu || activation == trelu || activation == softshrink
        noofoper = 1
    elseif activation == identity
        noofoper = 0
    else
        @info "Unkown activation type defaulting to identity"
        return 1
    end
    return noofoper
end

function gradientoperations(layer)
    activation = layer.σ
    if activation == sigmoid || activation == rrelu || activation == elu || activation == celu ||
        activation == softsign || activation == softplus || activation == tanh
        noofoper = 4
    elseif activation == gelu
        noofoper = 8
    elseif activation == hardsigmoid || activation == logsigmoid || activation == swish ||
            activation == selu || activation == softmax || activation == logsoftmax
        noofoper = 4
    elseif activation == hardtanh || activation == leakyrelu || activation == relu6 ||
            activation == lisht || activation == tanhshrink
        noofoper = 2
    elseif activation == logcosh || activation == mish
        noofoper = 5
    elseif activation == relu || activation == trelu || activation == softshrink
        noofoper = 1
    elseif activation == identity
        noofoper = 0
    else
        @info "Unkown activation type defaulting to identity"
        return 1
    end
    return noofoper
end

function recurrentunpack(layer::Recur)
    return layer.cell
end

islayer(::Any) = false
islayer(::Recur) = true
islayer(::Dense) = true
islayer(::Conv) = true
islayer(::MaxPool) = true
islayer(::MeanPool) = true
islayer(::DepthwiseConv) = true
islayer(::ConvTranspose) = true
islayer(::CrossCor) = true
islayer(::Maxout) = true


text:
# printing
using Printf, Statistics
# Oceananigans
using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.Utils
using Oceananigans.AbstractOperations
using Oceananigans.Advection

const hydrostatic = false

arch = CPU()
FT   = Float64

## units
const kilometer = 1000 # meters
const day = 86400      # seconds
const meter = 1
const hour = 60

## Domain
const Lx = 1000.0kilometer
const Ly = 2000.0kilometer
const Lz = 2985

# Discretization
Δt = 300.0
maxΔt = 300.0

end_time = 60day
advection   = WENO5()
timestepper = :RungeKutta3
# Rough target resolution
Δx = Δy = 5kilometer # 5km
Δz = 100meter
# Multiple of 16 gridpoints
const Nx = Int(192 * 0.5)
const Ny = Int(400 * 0.5)
const Nz = Int(32 * 0.5)
# Create Grid
topology = (Periodic, Bounded, Bounded)
grid = RectilinearGrid(arch, topology=topology,
                            size=(Nx, Ny, Nz),
                            x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

# Parameters
const f = -1e-4
const β = 1 * 10^(-11)
coriolis = FPlane(FT, f=f)
coriolis = BetaPlane(FT, f₀ = f, β = β)

α  = 2e-4     # [K⁻¹] Thermal expansion coefficient
g  = 9.8061   # [m/s²] gravitational constant
cᵖ = 3994.0   # [J/K]  heat capacity
ρ  = 999.8    # [kg/m³] density
const h = 1000.0     # [m] e-folding length scale for northern sponge
const ΔB = 8 * α * g # [m/s²] total change in buoyancy from surface to bottom
eos = LinearEquationOfState(FT, α=α, β=0)
buoyancy = BuoyancyTracer()

κh = 0.5e-5 # [m²/s] horizontal diffusivity
νh = 12.0   # [m²/s] horizontal viscocity
κv = 0.5e-5 # [m²/s] vertical diffusivity
νv = 3e-4   # [m²/s] vertical viscocity

vertical_closure = VerticalScalarDiffusivity(ν = νv, κ = κv)

horizontal_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh)

parameters = (
    Ly = Ly,                   # y-domain length
    τ = 0.2,                   # [N m⁻²] Zonal stress
    ρ = ρ,                     # [kg / m³]
    μ = 1.1e-3,                # [m/s]  linear drag
    H = Lz,                    # [m]
    h = h,                     # [m]    relexaction profile scale
    ΔB = ΔB,                   # [m/s²] buoyancy jump
    Lz = Lz,                   # [m]
    Lsponge = 900kilometer,   # [m]
    λᵗ = 7.0day,               # [s]
    Qᵇ = 10/(ρ * cᵖ) * α * g,  # [m² / s³]
    Qᵇ_cutoff = Ly * 5/6.      # [m]
)

# Momentum Boundary Conditions
@inline windshape(y, p) = sin( π * y / p.Ly)
@inline windstress(x, y, t, p) = - p.τ / p.ρ * windshape(y, p)
@inline ulineardrag(i, j, grid, clock, state, p) = @inbounds - p.μ * state.u[i, j, 1]
@inline vlineardrag(i, j, grid, clock, state, p) = @inbounds - p.μ * state.v[i, j, 1]

# Zonal Velocity
top_u_bc = BoundaryCondition(Flux, windstress, parameters = parameters)
bottom_u_bc =  BoundaryCondition(Flux, ulineardrag, discrete_form = true, parameters = parameters)
u_bcs = FieldBoundaryConditions(top = top_u_bc, bottom = bottom_u_bc)

# Meridional Velocity
bottom_v_bc =  BoundaryCondition(Flux, vlineardrag, discrete_form = true, parameters = parameters)
v_bcs = FieldBoundaryConditions(bottom = bottom_v_bc)

# Buoyancy Boundary Conditions Forcing. Note: Flux convention opposite of Abernathy
@inline cutoff(j, grid, p ) = grid.yᵃᶜᵃ[j] > p.Qᵇ_cutoff ? -0.0 : 1.0
@inline surface_flux(j, grid, p) = p.Qᵇ * cos(3π * grid.yᵃᶜᵃ[j] / p.Ly) * cutoff(j, grid, p)
@inline relaxation(i, j, grid, clock, state, p) = @inbounds surface_flux(j, grid, p)
top_b_bc = BoundaryCondition(Flux, relaxation, discrete_form = true, parameters = parameters)
b_bcs = FieldBoundaryConditions(top = top_b_bc)

# Save boundary conditions as named tuple
bcs = (b = b_bcs,  u = u_bcs, v = v_bcs,)

# Forcing Functions
# Sponge layers
relu(y) = (abs(y) + y) * 0.5
# Northern Wall Relaxation
@inline relaxation_profile_north(k, grid, p) = p.ΔB * ( exp(grid.zᵃᵃᶜ[k]/p.h) - exp(-p.Lz/p.h) ) / (1 - exp(-p.Lz/p.h))
function Fb_function(i, j, k, grid, clock, state, p)
    return @inbounds - (1/p.λᵗ)  * (state.b[i,j,k]
        - relaxation_profile_north(k, grid, p)) * relu( (grid.yᵃᶜᵃ[j]-p.Lsponge) / (p.Ly - p.Lsponge))
end
Fb = Forcing(Fb_function, parameters = parameters, discrete_form = true)

# Record forcings
forcings = (b = Fb, )

# Convective Parameterization
convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0,
                                                                background_κz = 5e-6,
                                                                background_νz = 3e-4)

# Model Setup
if hydrostatic
    model = HydrostaticFreeSurfaceModel(
            grid = grid,
            free_surface = ImplicitFreeSurface(),
            momentum_advection = WENO5(),
            tracer_advection = WENO5(),
            buoyancy = BuoyancyTracer(),
            coriolis = coriolis,
            closure = (horizontal_closure, vertical_closure, convective_adjustment),
            tracers = (:b,),
            boundary_conditions = bcs,
            forcing = forcings,
            )
else
    model = NonhydrostaticModel(
                    grid = grid,
                coriolis = coriolis,
                buoyancy = buoyancy,
                    closure = (horizontal_closure, vertical_closure),
                    tracers = (:b,),
        boundary_conditions = bcs,
                    forcing = forcings,
                advection = advection,
                timestepper = timestepper,
    )
end

# Timestep
Δt_wizard = TimeStepWizard(cfl = 1.0, Δt = Δt, max_change = 1.05, max_Δt = maxΔt)
cfl = AdvectiveCFL(Δt_wizard)

# Create Simulation
Ni = 1000
function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (model.clock.time / end_time)

    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    @printf("[%05.2f%%] i: %d, t: %.2e days, umax: (%6.3e, %6.3e, %6.3e) m/s, CFL: %6.4e, next Δt: %.1e s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), Δt_wizard.Δt)
    println(" ")
end

simulation = Simulation(model, Δt=Δt_wizard,
                        stop_time=end_time,
                        iteration_interval=Ni,
                        progress=print_progress)
## Run
run!(simulation)


text:
@testset "993.cousins-in-binary-tree.jl" begin
    t1 = TreeNode{Int}([1,2,3,4])
    @test is_cousins(t1, 4, 3) == false
    t1 = TreeNode{Int}([1,2,3,nothing,4,nothing,5])
    @test is_cousins(t1, 5, 4) == true
    t1 = TreeNode{Int}([1,2,3,nothing,4])
    @test is_cousins(t1, 2, 3) == false
end


text:
using KhepriRhino
using Test

@testset "KhepriRhino.jl" begin
    # Write your tests here.
end


text:
module SocialDynamics
export socdyn, socdyncoupad, socdyncoupexp, treefuncad, treefuncexp;

"""

    socdyn(nc, ns, a, p, steps, adfreq, seed, model)

Basic consistency-conformity social dynamics simulation function.

Arguments:

- `nc` : Number of non-expert citizens
- `ns` : number of experts
- `a` : number of attributes (belief, usage)
- `steps` : number of iterations each consisting of 1 interaction/person
- `adfreq` : advertising frequency
- `model` : expert influence for `model = "scientist"`, one-time media influence for `model = "ad"`, or regular media ads for `model = "regad"`

Returns:

- `n-tuse`: number of users

"""

function socdyn(nc, ns, a, p, steps, adfreq, seed, model)

srand(seed);
af = floor(adfreq*365.25/15);
tuse = Int64[];

if(model == "scientist")
	m = 1;
	n = ns + nc;	# Total number of the population

	# Set Scientist belief to 1
	soc = zeros(Bool, n,a);
	for j=1:ns
		soc[j,1]=true;
	end

elseif(model == "ad" || model == "regad")
	m = 2;
	if(model== "regad")
		m = 3;
	end
	n = nc;
	soc = zeros(Bool, n,a);
	# Informative Event
	for j=1:n
  		soc[j,1] = randbool();
	end
println("**************  End of Informative Event ***************")
end

for j = 1:steps
# Loop over all the population and update belief/usage

for i=1:n

	k=floor(rand()*n+1);

	   rcc = rand();          # sample from uniform distribution to choose CC rule
		if(m==1)
		      if(k<=ns)           # if scientist, only update second bit (usage bit)
			atr = 2;
		      else
			atr = floor(rand()*2)+1;
		      end
		elseif(m==2 || m==3)
			atr = floor(rand()*2)+1;
		else
			println("Error: unspecified model");
		end

	      if(rcc<p)           # Consistency Rule
		 soc[k,atr]= soc[k,!convert(Bool, atr-1)+1];
	      else                # Conformity Rule
		  rot = floor(rand()*n)+1;
		  soc[k,atr] = soc[rot,atr];
	      end
	end

tusage = sum(soc,1)[2];
append!(tuse , [tusage]);

if (m==3) # IF ADs are introduced regularly every certain number of steps
	if(j%af==0)
	  println("**************  NEW Informative Event ***************");
	print("steps number: ");
	println(j);
		  for j=1:n
		    soc[j,1] = randbool();
		  end
	end

# Demand << Supply
#	frac = 0.1;
#	poor = floor(frac*n);
#	if (tusage >n-1000)
#	  for j=1:poor
#	    soc[floor(rand()*n+1),2] = 0;
#	  end
#	end
end

if( tusage ==n )
  break;
end

end

return n-tuse

end

########################## COUPLED EXPERT MODEL #############
###############################################################
###############################################################

#
# Introduces the flexibility to toggle scientist belief
# Introcuces the flexibility to have non-zero initial usage and belief
# Input: expert belief (based on contamination level), social matrix belief and usage
# Belief: bool. True if belief in contamination.

"""
    socdyncoupexp(n, ns, belief, soc, p, steps)

Expert social dynamics.

Arguments:

- `n`: population size
- `ns` : number of experts
- `belief` : binary belief of water contamination
- `soc` : social record `[belief, usage, cont-accumulation, age]`.
- `p` : consistency rate; which gives the conformity rate `(1-p)`.
- `steps`: number of iterations each consisting of 1 interaction/person

Returns:

- updated `soc`

"""

function socdyncoupexp(n, ns, belief, soc, p, steps)

nc = n - ns;	# Total number of the non expert population
for j = 1:ns
	soc[j,1]=belief;
end

for j = 1:steps
# Loop over all the population and update belief/usage

	for i=1:n
		k = floor(Int, rand()*n+1);
		rcc = rand();          		# sample from uniform distribution to choose CC rule
		if(k <= ns)           		# if scientist, only update second bit (usage bit)
			atr = 2;
		else
			atr = floor(Int, rand()*2)+1;
		end

		if(rcc < p)          		# Consistency Rule
			soc[k,atr]= soc[k,!convert(Bool, atr-1)+1];
		else       			# Conformity Rule
			rot = floor(Int, rand()*n)+1;
			soc[k,atr] = soc[rot,atr];
	      end
	end
end

return soc

end

########################### COUPLED AD MODEL ##################
###############################################################
###############################################################

#
# Introduces the flexibility to toggle scientist belief
# Introcuces the flexibility to have non-zero initial usage and belief
# Input: expert belief (based on contamination level), social matrix belief and usage
# Belief: bool. True if belief in contamination.

"""
    socdyncoupad(n, adflag, belief, popfrac, soc, p, steps)

Media social dynamics.

Arguments:

- `n`: population size
- `adflag` : media influence if `adflag = 1`, do nothing otherwise.
- `belief` : binary belief of water contamination
- `popfrac` : fraction of population targetted with each ad campain.
- `soc` : social record `[belief, usage, cont-accumulation, age]`.
- `p` : consistency rate; which gives the conformity rate `(1-p)`.
- `steps`: number of iterations each consisting of 1 interaction/person

Returns:

- updated `soc`

"""

function socdyncoupad(n, adflag, belief, popfrac, soc, p, steps);

if(adflag == 1)
	for j = 1:floor(Int, popfrac*n)
		k = floor(Int, rand()*n+1)
		soc[k,1] = belief;
	end
end

for j = 1:steps
# Loop over all the population and update belief/usage

	for i = 1:n

		k = floor(Int, rand()*n+1);

		rcc = rand();  # sample from uniform distribution to choose CC rule
		atr = floor(Int, rand()*2)+1;

	      if(rcc<p)           # Consistency Rule
		 soc[k,atr] = soc[k,!convert(Bool, atr-1)+1];
	      else                # Conformity Rule
		  rot = floor(Int, rand()*n)+1;
		  soc[k,atr] = soc[rot,atr];
	      end
	end

end

return soc

end

###############################################
###############################################

#;# Inputs: simtime, pconsist, strat struct (popul), initial believers,
#;# sim = 9000;
#;# ss = [2 2 2 2];     #[indiv family community party]...
#;# iurate = 0.8;
#;# pconsist = 0.5;

"""

    treefuncexp(ss, expnum, soc, pconsist, simtime, iurate)

Expert hierarchical social dynamics according to interaction rate

Arguments:

- `ss` : vector for branching structure. For example `[2 2 2 2]` is a binary tree with 4 levels
- `expnum` : number of experts
- `soc` : social record `[belief, usage, cont-accumulation, age]`.
- `pconsist` : consistency rate
- `simtime` : simulation time (per interaction interval)
- `iurate` : initial usage rate

Returns:

- updated `soc`.

"""

function treefuncexp(ss, expnum, soc, pconsist, simtime, iurate)

strat = length(ss);  # number of stratification levels
popul = 1;
for i = 1:strat
	popul = popul * ss[i]; # Population count
end
peop = soc2tree(soc);
iupeop = iurate*popul;
for j = 1:iupeop
	k = ceil(Int, (rand()*popul)/2);
	peop[2*k] = 1;        # Choose iupopul to make use initially.
end

pconform = 1-pconsist;      # Conformity rate assumed to be subdivided all the rest
pmat = Array(Float32, strat); # probability matrix
pmat[1] = pconsist;         # first element in matrix is the lowest node

# Setting probability distribution
pp = pconform;
for i = 2:strat-1
	pp = pp/2;
	pmat[i] = pmat[i-1] + pp;
end
pmat[strat] = 1;

expindex = expnum*ss[1]; # Number of experts based on the number of attributes
	#;# Set the expert belief to 1 (for now).
for i = 1:expnum
	peop[(i-1)*ss[1]+1]=1;
end

for tt = 1:simtime
	for q = 1:popul

		j = floor(Int, popul*rand()+1);

		# If belief attribute of expert is chosen, reshuffle.
		while(mod(j,2)==1 && j<=expindex)
			j = floor(Int, popul*rand()+1);
		end
		branch = findbranch(pmat, rand()) + 1;        # Pick random branch according to pmat
		rglob = collect(findrange(branch, j, ss));        # Finds the range of under which branch falls

		rloc = findrange(branch-1, j, ss) - rglob[1] + 1;# Finds range that should be removed from br
		splice!(rglob, rloc);                        # Removes local branch from global
		otherind = rglob[floor(Int, rand()*length(rglob) + 1)];# pick a random attribute from the subbranch

		if(abs(otherind - j) >= ss[1])
			coro = 0;
			corj = 0;
			if(otherind%ss[1] == 0)        # Corr corrects for modulo when it's divisable
				coro = ss[1];
			end

			if(j%ss[1] == 0)
				corj = ss[1];
			end
			otherind = otherind + ( j%ss[1]+corj - (otherind%ss[1]+coro));
		end
		peop[j] = peop[otherind];                   # Match attributes
	end

	totcons = 0;
	for j = 1 : div(popul,2)
		totcons = totcons + peop[j*2];
	end

end

socresult = tree2soc(peop);
return socresult
end

##############################################################################
##############################################################################

"""

    treefuncad(ss, adflag, soc, pconsist, simtime, iurate, spreadrate)

Media hierarchical social dynamics according to interaction rate

Arguments:

- `ss` : vector for branching structure. For example `[2 2 2 2]` is a binary tree with 4 levels
- `adflag` : media influence if `adflag = 1`, do nothing otherwise.
- `soc` : social record `[belief, usage, cont-accumulation, age]`.
- `pconsist` : consistency rate
- `simtime` : simulation time (per interaction interval)
- `iurate` : initial usage rate
- `spreadrate` : portion of population affected by ads.

Returns:

- updated `soc`.

"""

function treefuncad(ss, adflag, soc, pconsist, simtime, iurate, spreadrate)

strat = length(ss);  # number of stratification levels
popul = 1;
for i = 1:strat
	popul = popul * ss[i];  # Population count
end
peop = soc2tree(soc);
iupeop = iurate * popul;
for j = 1:iupeop
	k = ceil(Int, (rand()*popul)/2);
	peop[2*k] = 1;        # Choose iupopul to make use initially.
end

pconform = 1-pconsist;      # Conformity rate assumed to be subdivided all the rest
pmat = Array(Float32, strat); # probability matrix
pmat[1] = pconsist;         # first element in matrix is the lowest node

# Setting probability distribution
pp = pconform;
for i = 2:strat-1
	pp = pp/2;
	pmat[i] = pmat[i-1] + pp;
end
pmat[strat] = 1;

if(adflag==1)
	for i=1:round(Int, popul*spreadrate)
		k = floor(Int, rand()*div(popul,2)) *2 + 1;
		peop[k] = 1;
	end
end

for tt = 1:simtime
	for q = 1:popul

		j = floor(Int, popul*rand() + 1);
		branch = findbranch(pmat, rand()) + 1;        	# Pick random branch according to pmat
		rglob = collect(findrange(branch, j, ss));        		# Finds the range of under which branch falls

		rloc = findrange(branch - 1, j, ss) - rglob[1] + 1;	# Finds range that should be removed from br
		splice!(rglob, rloc);                        		# Removes local branch from global
		otherind = rglob[floor(Int, rand()*length(rglob) + 1)];	# pick a random attribute from the subbranch
			if(abs(otherind - j) >= ss[1])
				coro = 0;
				corj = 0;

				if(otherind % ss[1] == 0)        # Corr corrects for modulo when it's divisable
					coro = ss[1];
        				end

        				if(j % ss[1] == 0)
          					corj = ss[1];
        				end
      				otherind = otherind + ( j%ss[1]+corj - (otherind%ss[1]+coro));
    				end

    				peop[j] = peop[otherind];                   # Match attributes
  			end

			totcons = 0;

  			for j = 1: div(popul, 2)
   				totcons = totcons + peop[j*2];
  			end

	end

socresult = tree2soc(peop);
return socresult

end

###############################################
########## TREE MODEL FUNCTIONS #############

function findbranch(A, num)
	An = A .- num;
	for k = 1:length(An)
		if(An[k] >= 0)
			return k
			break;
		end
	end
end

function findrange(branch, j, ss)
	span = 1;
	if(branch>0)
		for k = 1:branch-1
			span = span * ss[k];
		end
	end
	localind = ceil(Int, j/span);
	i0 = round(Int, (localind-1) * span + 1);
	range = i0 : (i0 + span - 1);
	return range
end

function findrange2(branch, j)
	span = 2^(branch-1);
	localind = ceil(Int, j/span);
	i0 = round(Int, (localind-1)*span + 1);
	range = i0 : (i0+span-1);
	return range
end

function soc2tree(soc)
	treevec = zeros(size(soc,1)*size(soc,2));
	for i = 1:size(soc,1)
		for j = 1:size(soc,2)
			treevec[(i-1)*size(soc,2) + j] = soc[i,j];
		end
	end
	return treevec
end

#########################################

function tree2soc(tree)
	soc = zeros(round(Int, length(tree)/2),2);
	for i = 1:size(soc,1)
		for j = 1:size(soc,2)
			soc[i,j] = tree[(i-1)*size(soc,2) + j];
		end
	end
	return soc
end

function cenergy(ss, peop, pmat)
	function baseenergy(por)
		delta = 2*sum(por)-length(por);
		return length(por)*exp(-delta);
	end

	ce = zeros(length(peop)/ss[1]);
	for i = 1:(length(peop)/ss[1])
		portion = peop[ss[1]*(i-1)+1:ss[1]*i];
		ce[i] =  baseenergy(portion);
	end

	for j = 2:(length(ss))
		for k = 1:(length(peop)/ss[j])

			tote += pmat[j]*sum(ce);
		end
	end

end

# # # #

end

max_stars_repo_name: josephbakarji/socio-hydro-abm

max_issues_repo_name: josephbakarji/socio-hydro-abm

max_forks_repo_name: josephbakarji/socio-hydro-abm


text:
abstract type Params end;

"""
    CompressionParams

Compression params

    CompressionParams(;
        compressor ::Symbol = default_compressor_name()
        level ::UInt8 = 5
        typesize ::Int32 = 8
        nthreads ::Int16 = 1
        blocksize ::Int32 = 0
        splitmode ::Bool = false
        filter_pipeline ::FilterPipeline = filter_pipeline(:shuffle)
    )

Create compression parameters

# Arguments

- `compressor::Symbol` - The name of compressor to use
- `level` - The compression level from 0 (no compression) to 9 (maximum compression)
- `typesize` - The size of type being compressed
- `nthreads` - The number of threads to use internally
- `blocksize` - The requested size of the compressed blocks (0 means auto)
- `splitmode` - Whether the blocks should be split or not.
- `filter_pipeline` - Filters pipeline
"""
struct CompressionParams <: Params
    compressor ::Compressor
    level ::UInt8
    typesize ::Int32
    nthreads ::Int16
    blocksize ::Int32
    splitmode ::Bool
    filter_pipeline ::FilterPipeline
    function CompressionParams(;
        compressor ::Symbol = default_compressor_name(),
        level = 5,
        typesize = 8,
        nthreads = 1,
        blocksize = 0,
        splitmode = false,
        filter_pipeline = filter_pipeline(:shuffle)
    )
        !(level in 0:9) && throw(ArgumentError("level must be in 0:9 range"))
        return new(
            compressor_by_name(compressor),
            level,
            typesize,
            nthreads,
            blocksize,
            splitmode,
            filter_pipeline
            )
    end
end

"""
    CompressionParams(::Type{T}; kwargs...)

Create `CompressionParams` for compressing vectors with element type `T`.
This is equivalent to `CompressionParams(;typesize = sizeof(T), kwargs...)`
"""
function CompressionParams(::Type{T}; kwargs...) where {T}
    return CompressionParams(;
        typesize = sizeof(T),
        kwargs...
    )
end

function make_cparams(p::CompressionParams)
    return Lib.blosc2_cparams(
        p.compressor.code,
        0,
        p.level,
        0,
        p.typesize,
        p.nthreads,
        p.blocksize,
        p.splitmode ? 1 : 0,
        C_NULL,
        getproperty.(p.filter_pipeline.filters, :id),
        getproperty.(p.filter_pipeline.filters, :meta),
        C_NULL,
        C_NULL,
        C_NULL
    )
end

"""
    DecompressionParams

Decompression params

    DecomplessionParams(;nthreads = 1)

Create decompression parameters
"""
struct DecompressionParams <: Params
    nthreads ::Cint
    function DecompressionParams(;nthreads = 1)
        new(nthreads)
    end
end

function make_dparams(p::DecompressionParams)
    return Lib.blosc2_dparams(
        p.nthreads, C_NULL, C_NULL, C_NULL
    )
end


text:
add_assert_scope(:padic_poly)

@doc Markdown.doc"""
    lift(a::T, K::PadicField) where T <: Union{Nemo.nmod, Generic.Res{fmpz}, gfp_elem} -> padic

Computes a lift of the element from the residue ring.
"""
function lift(a::T, K::PadicField) where T <: Union{Nemo.nmod, Nemo.fmpz_mod, Generic.Res{fmpz}, gfp_elem}
  n = modulus(parent(a))
  p = prime(K)
  v, fl = remove(n, p)
  @assert isone(fl)
  return Hecke.lift(a) + O(K, p^v)
end

@doc Markdown.doc"""
    lift(f::T, Kt) where T <: Union{nmod_poly, fmpz_mod_poly, gfp_poly} -> Generic.Poly{padic}

Computes a lift of the polynomial lifting every coefficient of the residue ring.
"""
function lift(f::T, Kt) where T <: Union{nmod_poly, fmpz_mod_poly, gfp_poly}
  K = base_ring(Kt)
  coeffs = Vector{padic}(undef, degree(f)+1)
  for i = 1:degree(f)+1
    coeffs[i] = lift(coeff(f, i-1), K)
  end
  return Kt(coeffs)
end

@doc Markdown.doc"""
    lift(x::fq_nmod, Q::QadicField) -> qadic

Computes a lift of the element from the residue ring.
"""
function lift(x::fq_nmod, Q::QadicField)
  z = Q()
  for i=0:degree(Q)-1
    setcoeff!(z, i, coeff(x, i))
  end
  return setprecision(z, 1)
end

@doc Markdown.doc"""
    lift(x::fq_nmod_poly, Kt) -> Generic.Poly{qadic}

Computes a lift of the polynomial lifting every coefficient of the residue ring.
"""
function lift(x::fq_nmod_poly, Kt)
  K = base_ring(Kt)
  coeffs = Vector{qadic}(undef, degree(x)+1)
  for i = 1:degree(x)+1
    coeffs[i] = lift(coeff(x, i-1), K)
  end
  return Kt(coeffs)
end

function _content(f::Generic.Poly{T}) where T <: Union{padic, qadic}
  K = base_ring(f)
  p = uniformizer(K)
  v = valuation(coeff(f, 0))
  for i = 1:degree(f)
    v = min(v, valuation(coeff(f, i)))
  end
  return p^v
end

function rem!(x::AbstractAlgebra.Generic.Poly{T}, y::AbstractAlgebra.Generic.Poly{T}, z::AbstractAlgebra.Generic.Poly{T}) where T <:Union{padic, qadic}
  x = rem(y, z)
  return x
end

################################################################################
#
#  Fun factor
#
################################################################################

function fun_factor(g::Generic.Poly{padic})
  K = base_ring(g)
  Kt = parent(g)
  v = precision(g)
  pv = prime(K)^v
  R = ResidueRing(FlintZZ, pv, cached = false)
  Rt = PolynomialRing(R, "t", cached = false)[1]
  fR = Rt([R(Hecke.lift(coeff(g, i))) for i = 0:degree(g)])
  u, g1 = Hecke.fun_factor(fR)
  return lift(u, Kt), lift(g1, Kt)
end

function fun_factor(f::Generic.Poly{qadic})
  K = base_ring(f)
  Kt = parent(f)
  v = precision(f)
  @assert isone(_content(f))
  if iszero(valuation(leading_coefficient(f)))
    return one(Kt), g
  end
  ind = degree(f) -1
  while !iszero(valuation(coeff(f, ind)))
    ind -= 1
  end
  g = Kt([coeff(f, i) for i = ind:degree(f)])
  h = Kt([divexact(coeff(f, i), coeff(f, ind)) for i = 0:ind])
  s = Kt(inv(coeff(g, 0)))
  t = zero(Kt)
  k = Int(clog(fmpz(v), 2))+1
  for i = 1:k
    e = f - g*h
    q, r = divrem(s*e, h)
    gn = g+t*e+q*g
    hn = h+r
    b = s*gn+t*hn-1
    c, d = divrem(s*b, hn)
    sn = s - d
    tn = t- t*b- c*gn
    g = gn
    h = hn
    s = sn
    t = tn
  end
  return g, h
end

################################################################################
#
#  Gcd
#
################################################################################

function Nemo.precision(g::Generic.Poly{T}) where T <: Union{padic, qadic}
  N = precision(coeff(g, 0))
  for i = 1:degree(g)
    N = min(N, precision(coeff(g, i)))
  end
  return N
end

function Base.gcd(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic}
  if degree(f) < degree(g)
    f, g = g, f
  end
  while true
    cf = _content(f)
    if !isone(cf)
      f = divexact(f, cf)
    end
    cg = _content(g)
    if !isone(cg)
      g = divexact(g, cg)
    end
    if !iszero(valuation(leading_coefficient(g)))
      u, g1 = fun_factor(g)
      if iszero(valuation(leading_coefficient(f)))
        g = g1#*reverse(gcd(reverse(f), reverse(u)))
      else
        v, f1 = fun_factor(f)
        return reverse(gcd(reverse(u), reverse(v)))*gcd(f1, g1)
      end
    end
    f = mod(f, g)
    if degree(f) < 1
      if iszero(f)
        return g
      else
        return divexact(f, leading_coefficient(f))
      end
    else
      f, g = g, f
    end
  end
end

################################################################################
#
#  Invmod
#
################################################################################

function invmod(u::Generic.Poly{padic}, f::Generic.Poly{padic})
  if !iszero(valuation(leading_coefficient(f)))
    error("Not yet implemented")
  end
  if !iszero(valuation(coeff(u, 0))) || !all(x -> x > 0, [valuation(coeff(u, i)) for i = 1:degree(u)])
    g, s, t = gcdx(u, f)
    return s
  end
  K = base_ring(f)
  Kt = parent(f)
  v = min(precision(f), precision(u))
  pv = prime(K)^v
  R = ResidueRing(FlintZZ, pv, cached = false)
  Rt = PolynomialRing(R, "t", cached = false)[1]
  fR = Rt([R(Hecke.lift(coeff(f, i))) for i = 0:degree(f)])
  uR = Rt([R(Hecke.lift(coeff(u, i))) for i = 0:degree(u)])
  iuR = invmod(uR, fR)
  return lift(iuR, Kt)
end

function invmod(f::Generic.Poly{qadic}, M::Generic.Poly{qadic})
  if !iszero(valuation(leading_coefficient(M)))
    error("Not yet implemented")
  end
  f = rem(f, M)
  if !iszero(valuation(coeff(f, 0))) || !all(x -> x > 0, [valuation(coeff(f, i)) for i = 1:degree(f)])
    g, s, t = gcdx(f, M)
    return s
  end
  K = base_ring(f)
  Kt = parent(f)
  v = min(precision(f), precision(M))
  g = parent(f)(inv(constant_coefficient(f)))
  c = f*g
  c = rem!(c, c, M)
  while !isone(c)
    g = mul!(g, g, 2-c)
    g = rem!(g, g, M)
    c = mul!(c, f, g)
    c = rem!(c, c, M)
  end
  return g
end

################################################################################
#
#  xgcd
#
################################################################################

#TODO: The implementation is recursive. Change it to an iterative implementation.
function gcdx(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic}
  if degree(f) < degree(g)
    r1, r2, r3 = gcdx(g, f)
    return r1, r3, r2
  end
  Kx = parent(f)
  if length(g) <= 1
    if iszero(g)
      return f, one(Kx), zero(Kx)
    else
      s = Kx(inv(coeff(g, 0)))
      @hassert one(Kx) == s*g
      return one(Kx), zero(Kx), s
    end
  end
  cf = _content(f)
  if !isone(cf)
    f1 = divexact(f, cf)
    d, u, v = gcdx(f1, g)
    @hassert :padic_poly 1 f*divexact(u, cf) + v*g == d
    return d, divexact(u, cf), v
  end
  cg = _content(g)
  if !isone(cg)
    g1 = divexact(g, cg)
    d, u, v = gcdx(f, g1)
    @hassert :padic_poly 1  f*u+divexact(v, cg)*g == d
    return d, u, divexact(v, cg)
  end
  if iszero(valuation(leading_coefficient(g)))
    q, f1 = divrem(f, g)
    d, u, v = gcdx(g, f1)
    @hassert :padic_poly 1  d == f*v+(u-v*q)*g
    return d, v, u-v*q
  end
  ug, gg = fun_factor(g)
  if iszero(valuation(leading_coefficient(f)))
    s = invmod(ug, f)
    t = divexact(1-s*ug, f)
    @hassert :padic_poly 1  t*f == 1-s*ug
    d, u, v = gcdx(f, gg)
    @hassert :padic_poly 1  d == u*f + v*gg
    @hassert :padic_poly 1  d == (u+v*t*gg)*f + v*s*g
    return d, u+v*t*gg, v*s
  end
  uf, ff = fun_factor(f)
  d, u, v = gcdx(ff, gg)
  if degree(gg) >= 1
    s = invmod(uf, gg)
    t = divexact(1-s*uf, gg)
    @hassert :padic_poly 1  t*gg == 1-s*uf
  else
    #gg = 1. Easy to compute Bezout coefficients...
    s = zero(Kx)
    t = one(Kx)
  end
  U = u*s
  V = u*ff*t + v*t*gg+s*uf*v
  d1, u1, v1 = gcdx(reverse(uf), reverse(ug))
  d1 = reverse(d1)
  u1 = reverse(u1)
  v1 = reverse(v1)
  @hassert :padic_poly 1  d1 == u1*uf+v1*ug
  if degree(ff) >= 1
    t1 = invmod(ug, ff)
    s1 = divexact(1-t1*ug, ff)
    @hassert :padic_poly 1  t1*ug + s1*ff == one(Kx)
  else
    t1 = zero(Kx)
    s1 = one(Kx)
    @hassert :padic_poly 1  s1*ff + t1*ug == one(Kx)
  end
  U1 = u1*s1
  V1 = s1*ff*v1+t1*u1*uf+t1*v1*ug

  DD = d*d1
  UU = U*U1*f + U1*V*gg+U*V1*ug
  VV = V*V1
  @hassert :padic_poly 1  DD == UU*f + VV*g
  return DD, UU, VV
end

function divexact(f::AbstractAlgebra.PolyElem{T}, g::AbstractAlgebra.PolyElem{T}) where T<: Union{padic, qadic}
   check_parent(f, g)
   f1 = deepcopy(f)
   g1 = deepcopy(g)
   iszero(g) && throw(DivideError())
   if iszero(f)
      return zero(parent(f))
   end
   lenq = length(f) - length(g) + 1
   d = Array{T}(undef, lenq)
   for i = 1:lenq
      d[i] = zero(base_ring(f))
   end
   x = gen(parent(f))
   leng = length(g)
   while length(f) >= leng
      lenf = length(f)
      q1 = d[lenf - leng + 1] = divexact(coeff(f, lenf - 1), coeff(g, leng - 1))
      f = f - shift_left(q1*g, lenf - leng)
      if length(f) == lenf # inexact case
         set_length!(f, normalise(f, lenf - 1))
      end
   end
   q = parent(f)(d)
   set_length!(q, lenq)
   K = base_ring(f)
   Kt = parent(f)
   p = prime(K)
   while !iszero(q*g1 - f1)
     q = Kt(T[coeff(q, i) + O(K, p^(precision(q)-1)) for i = 0:degree(q)])
   end
   return q
end

################################################################################
#
#  Resultant and characteristic polynomial
#
################################################################################

function rres(f::Generic.Poly{padic}, g::Generic.Poly{padic})
  Kt = parent(f)
  K = base_ring(Kt)
  p = prime(K)
  v = min(precision(f), precision(g))
  R = ResidueRing(FlintZZ, p^v, cached = false)
  cf = Vector{elem_type(R)}(undef, degree(f)+1)
  for i = 1:length(cf)
    cf[i] = R(Hecke.lift(coeff(f, i-1)))
  end
  cg = Vector{elem_type(R)}(undef, degree(g)+1)
  for i = 1:length(cg)
    cg[i] = R(Hecke.lift(coeff(g, i-1)))
  end
  Rt = PolynomialRing(R, "t", cached = false)[1]
  r = Hecke.rres_sircana_pp(Rt(cf), Rt(cg))
  return lift(r, K)
end

function resultant(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic}
  Nemo.check_parent(f, g)
  #First, we need to make the polynomials integral
  Rt = parent(f)
  R = base_ring(Rt)
  res = one(R)
  c1 = _content(f)
  if valuation(c1) < 0
    res *= c1^degree(g)
    f = divexact(f, c1)
  end
  c2 = _content(g)
  if valuation(c2) < 0
    res *= c2^degree(f)
    g = divexact(g, c2)
  end
  return res * _resultant(f, g)
end

function _resultant(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic}
  Rt = parent(f)
  R = base_ring(Rt)
  res = one(R)
  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        return zero(R)
      end
      return res
    elseif degree(f) < 1
      res *= leading_coefficient(f)^degree(g)
      return res
    elseif degree(g) < 1
      res *= leading_coefficient(g)^degree(f)
      return res
    end

    cf = _content(f)
    if !isone(cf)
      f = divexact(f, cf)
      res *= cf^degree(g)
    end

    cg = _content(g)
    if !isone(cg)
      g = divexact(g, cg)
      res *= cg^degree(f)
    end

    if degree(f) < degree(g)
      if !iszero(mod(degree(f)*degree(g), 2))
        res = -res
      end
      f, g = g, f
    end

    if valuation(leading_coefficient(g)) == 0
      r = rem(f, g)
      res *= leading_coefficient(g)^(degree(f) - degree(r))
      if !iszero(mod(degree(g)*(degree(f) - degree(r)), 2))
        res = -res
      end
      f = r
    else
      break
    end
  end
  g1, g2 = fun_factor(g)
  res1 = _resultant(f, g2)
  g1r = reverse(g1)
  fr = reverse(f)
  res2 = (constant_coefficient(g1)^(degree(f) - degree(fr)))*_resultant(fr, g1r)
  return res*res1*res2
end

degree(::FlintPadicField) = 1

@doc Markdown.doc"""
    characteristic_polynomial(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic} -> Generic.Poly{T}

Computes $\mathrm{Res}_x(f(x), t- g(x))$.
"""
function characteristic_polynomial(f::Generic.Poly{padic}, g::Generic.Poly{padic})
  Kt = parent(f)
  Ktx, x = PolynomialRing(Kt, "x")
  fcoeffs = Generic.Poly{padic}[Kt(coeff(f, i)) for i = 0:degree(f)]
  gcoeffs = Generic.Poly{padic}[Kt(-coeff(g, i)) for i = 0:degree(g)]
  f1 = Ktx(fcoeffs)
  g1 = Ktx(gcoeffs) + gen(Kt)
  return resultant(f1, g1)
end

#=
function characteristic_polynomial(f::Generic.Poly{T}, g::Generic.Poly{T}) where T <: Union{padic, qadic}
  K = base_ring(f)
  Kt = parent(f)
  p = prime(K)
  d = degree(K)
  if p^d <= degree(f)
    if d > 1
      error("Not yet implemented")
    end
    d1 = clog(fmpz(degree(f)+1), p)
    L = QadicField(p, d1, min(precision(f), precision(g)))
    Lt = PolynomialRing(L, "t")[1]
    fL = change_base_ring(f, L, Lt)
    gL = change_base_ring(g, L, Lt)
    cp = characteristic_polynomial(fL, gL)
    #cp will have coefficients over K, so I need to change the base ring.
    cf = [coeff(coeff(cp, i), 0) for i = 0:degree(cp)]
    return Kt([setprecision(K(lift(cf[i])), precision(cf[i])) for i = 1:length(cf)])
  end
  #The resultant will be a polynomial of degree degree(f1). So I need degree(f1)+1 interpolation points.
  ev_points = interpolation_points(K, degree(f)+1)
  res = Vector{typeof(f)}(undef, degree(f)+1)
  for i = 1:length(res)
    res[i] = Kt(resultant(f, ev_points[i] - g))
  end
  int_ideals = Vector{typeof(f)}(undef, length(ev_points))
  for i = 1:length(int_ideals)
    int_ideals[i] = Kt(T[-ev_points[i], K(1)])
  end
  crtctx = crt_env(int_ideals)
  resu = crt(res, crtctx)
  return resu
end
=#

function interpolation_points(K, n::Int)
  p = prime(K)
  if n < p
    return [K(i) for i = 1:n]
  end
  int_points = Vector{elem_type(K)}(undef, n)
  k, mk = ResidueField(K)
  ind = 1
  for el in k
    int_points[ind] = mk\el
    ind += 1
    if ind > n
      break
    end
  end
  @assert ind > n
  return int_points

end

################################################################################
#
#  Hensel factorization
#
################################################################################
@doc Markdown.doc"""
    Hensel_factorization(f::Generic.Poly{T}) where T <: Union{padic, qadic} -> Dict{Generic.Poly{T}, Generic.Poly{T}}

Computes a factorization of $f$ such that every factor has a unique irreducible factor over the residue field.
The output is a dictionary whose keys are lifts of the irreducible factors over the residue field and values the corresponding factors of $f$.
"""
function Hensel_factorization(f::Generic.Poly{T}) where T <: Union{padic, qadic}
  D = Dict{Generic.Poly{T}, Generic.Poly{T}}()
  Kt = parent(f)
  K = base_ring(Kt)
  k, mk = ResidueField(K)
  kt = PolynomialRing(k, "t", cached = false)[1]
  fp = kt([mk(coeff(f, i)) for i = 0:degree(f)])
  lfp = factor(fp).fac
  if length(lfp) == 1
    #The Hensel factorization is trivial...
    phi = setprecision(lift(first(keys(lfp)), Kt), precision(f))
    D[phi] = f
    return D
  end
  vlfp = Vector{dense_poly_type(elem_type(k))}(undef, length(lfp))
  ind = 1
  ks = Vector{Generic.Poly{T}}(undef, length(vlfp))
  for (k1, v) in lfp
    vlfp[ind] = k1^v
    ks[ind] = setprecision(lift(k1, Kt), precision(f))
    ind += 1
  end
  H = HenselCtxdr{T}(f, vlfp)
  lift(H, precision(f))
  for i = 1:H.n
    D[ks[i]] = H.lf[i]
  end
  return D
end

mutable struct HenselCtxdr{S}
  f::PolyElem{S}
  lf::Array{PolyElem{S}, 1}
  la::Array{PolyElem{S}, 1}
  p::S
  n::Int

  function HenselCtxdr{qadic}(f::Generic.Poly{qadic}, lfp::Vector{Generic.Poly{qadic}}, la::Vector{Generic.Poly{qadic}}, p::qadic, n::Int)
    return new(f, lfp, la, p, n)
  end

  function HenselCtxdr{S}(f::PolyElem{S}, lfp::Vector{T}) where {S, T}
    @assert sum(map(degree, lfp)) == degree(f)
    Q = base_ring(f)
    Qx = parent(f)
    i = 1
    la = Array{PolyElem{S}, 1}()
    n = length(lfp)
    while i < length(lfp)
      f1 = lfp[i]
      f2 = lfp[i+1]
      g, a, b = gcdx(f1, f2)
      @assert isone(g)
      push!(la, map_coefficients(x -> setprecision(lift(x, Q), 1), a, parent = Qx))
      push!(la, map_coefficients(x -> setprecision(lift(x, Q), 1), b, parent = Qx))
      push!(lfp, f1*f2)
      i += 2
    end
    return new(f, map(x -> map_coefficients(y -> setprecision(lift(y, Q), 1), x, parent = Qx), lfp), la, uniformizer(Q), n)
  end

  function HenselCtxdr{S}(f::PolyElem{S}) where S
    Q = base_ring(f)
    K, mK = ResidueField(Q)
    fp = change_base_ring(f, mK)
    fac = factor(fp).fac
    lfp = Vector{keytype(fac)}(undef, length(fac))
    ind = 1
    for (k, v) in fac
      lfp[ind] = k^v
      ind += 1
    end
    return HenselCtxQadic{S}(f, lfp)
  end
end

function lift(C::HenselCtxdr, mx::Int)
  p = C.p
  N = valuation(p)
#  @show map(precision, coefficients(C.f)), N, precision(parent(p))
  #have: N need mx
  ch = Int[mx]
  while ch[end] > N
    push!(ch, div(ch[end]+1, 2))
  end
  @vprint :PolyFactor 1 "using lifting chain ", ch
  for k=length(ch)-1:-1:1
    N2 = ch[k]
    i = length(C.lf)
    j = i-1
    p = setprecision(p, N2)
    while j > 0
      if i == length(C.lf)
        f = setprecision(C.f, N2)
      else
        f = setprecision(C.lf[i], N2)
      end
      #formulae and names from the Flint doc
      h = C.lf[j]
      g = C.lf[j-1]
      b = C.la[j]
      a = C.la[j-1]
      h = setprecision(h, N2)
      g = setprecision(g, N2)
      a = setprecision(a, N2)
      b = setprecision(b, N2)
      fgh = (f-g*h)*inv(p)
      G = rem(fgh*b, g)*p+g
      H = rem(fgh*a, h)*p+h
      t = (1-a*G-b*H)*inv(p)
      B = rem(t*b, g)*p+b
      A = rem(t*a, h)*p+a
      if i < length(C.lf)
        C.lf[i] = G*H
      end
      C.lf[j-1] = G
      C.lf[j] = H
      C.la[j-1] = A
      C.la[j] = B
      i -= 1
      j -= 2
    end
  end
  return nothing
end

################################################################################
#
#  Slope factorization
#
################################################################################

@doc Markdown.doc"""
    slope_factorization(f::Generic.Poly{T}) where T <: Union{padic, qadic} -> Dict{Generic.Poly{T}, Int}

Computes a factorization of $f$ such that every factor has a one-sided generalized Newton polygon.
The output is a dictionary whose keys are the factors of $f$ and the corresponding value is the multiplicity.
"""
function slope_factorization(f::Generic.Poly{T}) where T <: Union{padic, qadic}

  K = base_ring(f)
  Kt = parent(f)
  fact = Dict{Generic.Poly{T}, Int}()
  cf = _content(f)
  f = divexact(f, cf)
  if !iszero(valuation(leading_coefficient(f)))
    u, f = fun_factor(f)
    u1 = reverse(u)
    sf = slope_factorization(u1)
    for (ff, eff) in sf
      fact[reverse(ff)] = eff
    end
  end
  sqf = factor_squarefree(f)
  for (g, v) in sqf
    hg = Hensel_factorization(g)
    for (phi, fphi) in hg
      if degree(phi) == degree(fphi)
        fact[fphi] = v
        continue
      end
      NP = NewtonPolygon(fphi, phi)
      L = lines(NP)
      L1 = sort(L, rev = true, by = x -> slope(x))
      for l in L1
        if l == L1[end]
          fact[fphi] = v
          break
        end
        s = slope(l)
        mu = divexact(phi^Int(denominator(s)), uniformizer(K)^(-(Int(numerator(s)))))
        chi = characteristic_polynomial(fphi, mu)
        hchi = Hensel_factorization(chi)
        for (ppp, fff) in hchi
          if ppp == gen(Kt)
            continue
          end
          com = fff(mu)
          gc = gcd(com, fphi)
          fact[gc] = v
          fphi1 = divexact(fphi, gc)
          if gc*fphi1 != fphi
            error("problem!")
          end
          fphi = fphi1
        end
      end
    end
  end
  return fact
end


text:
const FIXED_TIME_ZONE_REGEX = r"""
    ^(?|
        Z
    |
        UTC
        (?:
            (?<sign>[+-])
            (?<hour>\d{1,2})
        )?
    |
        (?<sign>[+-])
        (?<hour>\d{2})
    |
        (?:UTC(?=[+-]))?
        (?<sign>[+-])?
        (?<hour>\d{2})
        (?|
            (?(hour)\:(?<minute>\d{2}))
            (?(minute)\:(?<second>\d{2}))?
        |
            (?(hour)(?<minute>\d{2}))
        )
    )$
    """x

"""
    FixedTimeZone

A `TimeZone` with a constant offset for all of time.
"""
struct FixedTimeZone <: TimeZone
    name::String
    offset::UTCOffset
end

"""
    FixedTimeZone(name, utc_offset, dst_offset=0) -> FixedTimeZone

Constructs a `FixedTimeZone` with the given `name`, UTC offset (in seconds), and DST offset
(in seconds).
"""
function FixedTimeZone(name::AbstractString, utc_offset::Integer, dst_offset::Integer=0)
    FixedTimeZone(name, UTCOffset(utc_offset, dst_offset))
end

function FixedTimeZone(name::AbstractString, utc_offset::Second, dst_offset::Second=Second(0))
    FixedTimeZone(name, UTCOffset(utc_offset, dst_offset))
end

# https://en.wikipedia.org/wiki/ISO_8601#Coordinated_Universal_Time_(UTC)
const UTC_ZERO = FixedTimeZone("Z", 0)

"""
    FixedTimeZone(::AbstractString) -> FixedTimeZone

Constructs a `FixedTimeZone` with a UTC offset from a string. Resulting `FixedTimeZone` will
be named like \"UTC±HH:MM[:SS]\".

# Examples
```julia
julia> FixedTimeZone(\"UTC+6\")
UTC+06:00

julia> FixedTimeZone(\"-1330\")
UTC-13:30

julia> FixedTimeZone(\"15:45:21\")
UTC+15:45:21
```
"""
function FixedTimeZone(s::AbstractString)
    s == "Z" && return UTC_ZERO

    m = match(FIXED_TIME_ZONE_REGEX, s)
    m === nothing && throw(ArgumentError("Unrecognized time zone: $s"))

    coefficient = m[:sign] == "-" ? -1 : 1
    sig = coefficient < 0 ? '-' : '+'
    hour = m[:hour] === nothing ? 0 : parse(Int, m[:hour])
    minute = m[:minute] === nothing ? 0 : parse(Int, m[:minute])
    second = m[:second] === nothing ? 0 : parse(Int, m[:second])

    if hour == 0 && minute == 0 && second == 0
        name = "UTC"
    elseif second == 0
        name = @sprintf("UTC%c%02d:%02d", sig, hour, minute)
    else
        name = @sprintf("UTC%c%02d:%02d:%02d", sig, hour, minute, second)
    end

    offset = coefficient * (hour * 3600 + minute * 60 + second)
    return FixedTimeZone(name, offset)
end

name(tz::FixedTimeZone) = tz.name
rename(tz::FixedTimeZone, name::AbstractString) = FixedTimeZone(name, tz.offset)


text:
module DigilentWaveForms

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("DigilentWaveForms was not built properly. Please run Pkg.build(\"DigilentWaveForms\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    # check_deps()
end

include(joinpath(@__DIR__, "..", "gen", "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libdwf_fixes.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_consts.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_api.jl"))

include("utils.jl")
include("enumeration.jl")
include("system.jl")

include("powersupplies.jl")
include("analogin.jl")
include("analogout.jl")

# export everything
#foreach(names(@__MODULE__, all=true)) do s
#    if startswith(string(s), "SOME_PREFIX")
#        @eval export $s
#    end
#end

end # module


text:
module SBPFD

end # module


text:
using Gumbo
using Requests
import Requests: post, get

const uri="https://utslogin.nlm.nih.gov"
const auth_endpoint = "/cas/v1/tickets/"
const service="http://umlsks.nlm.nih.gov"
const rest_uri = "https://uts-ws.nlm.nih.gov"

immutable BadResponseException <: Exception
    code::Int64
end

function Base.showerror(io::IO, e::BadResponseException)
    print(io, "BadResponseException with code: ",STATUS_CODES[e.code])
end

"""
    Credentials(user, psswd)
"""
type Credentials
    "username"
    username::ASCIIString
    "password"
    password::ASCIIString
end

function time_to_last_save(file)
    #unix time is GMT
    time_diff = Int(Dates.unix2datetime(time()) - Dates.unix2datetime(mtime(file)))/ (1000 * 60 * 60)
    return time_diff
end

"""
    get_tgt(c::Credentials)

Retrieve a ticket granting ticket (TGT) using UTS login credentials
A tgt is valid for 8 hours. Therefore, look for UTS_TGT.txt in the local
directory to see if it has been recently stored.
"""
function get_tgt(c::Credentials)
    #Commented out due to old TGT not working
    #Check if there is a valid ticket on disk
    # TGT_file = "UTS_TGT.txt"
    # if isfile(TGT_file)
    #     #check time
    #     time_elapsed = time_to_last_save(TGT_file)
    #     # Expiration time should be 8 hours - but I tend to expirience bad TGT after few hours
    #     if time_elapsed > 2.5
    #         println("UTS TGT Expired")
    #         rm(TGT_file)
    #     else
    #         # println("Using TGT from disk - saved ", time_elapsed, " hours ago")
    #         fin = open(TGT_file)
    #         lines = readlines(fin)
    #         ticket = lines[1]
    #         return ticket
    #     end
    #
    # end

    println("Requesting new TGT")
    params = Dict("username" => c.username, "password"=> c.password)
    headers = Dict("Content-type"=> "application/x-www-form-urlencoded",
    "Accept"=> "text/plain", "User-Agent"=>"julia" )
    r = post(uri*auth_endpoint,data=params,headers=headers)
    ascii_r = ASCIIString(r.data)

    doc = parsehtml(ascii_r)
    #for now - harcoded
    #TO DO:: parse and check
    try
        ticket = getattr(doc.root.children[2].children[2], "action")
        return ticket
    catch
        error("Could not get TGT: UTS response structure is wrong")
    end
    # fout = open(TGT_file, "w")
    # write(fout, ticket)
end

"""
    get_ticket(tgt)

Retrieve a single-use Service Ticket using TGT
"""
function get_ticket(TGT)
    params = Dict("service"=> service)
    h = Dict("Content-type"=> "application/x-www-form-urlencoded",
    "Accept"=> "text/plain", "User-Agent"=>"julia" )
    r = post(TGT; data=params, headers=h)
    return ASCIIString(r.data)
end

"""
    search_umls(c::Credentials, query)

Search UMLS Rest API. For more info see
[UMLS_API](https://documentation.uts.nlm.nih.gov/rest/search/)

###Arguments

- `c::Credentials`: UMLS username and password
- `query`: UMLS query containing the search term
- `version:` Optional - defaults to current

###Output

- `result_pages`: Array, where each entry is a dictionary containing a page of
results. e.g
 ```Dict{AbstractString,Any} with 3 entries:
  "pageSize"   => 25
   "pageNumber" => 1
  "result"     => Dict{AbstractString,Any}("classType"=>"searchResults","result…
  ```

###Examples

```julia
credentials = Credentials(user, psswd)
tgt = get_tgt(credentials)
term = "obesity"
query = Dict("string"=>term, "searchType"=>"exact" )
all_results= search_umls(tgt, query)
```
"""
function search_umls(tgt, query; version::ASCIIString="current", timeout=1)

    # Ticket granting ticket
    if tgt== nothing
        tgt = get_tgt(c)
    end
    page=0

    content_endpoint = "/rest/search/current"

    #each page of results is appended to the output list
    #where each entry is a dictionary containing that pages's results
    # e.g
    # Dict{AbstractString,Any} with 3 entries:
    #   "pageSize"   => 25
    #   "pageNumber" => 1
    #   "result"     => Dict{AbstractString,Any}("classType"=>"searchResults","result…
    result_pages = Array{Any,1}()

    while true

        #get a new ticket per page if necessary
        ticket = get_ticket(tgt)
        page +=1
        #append ticket to query
        query["ticket"]= ticket
        query["pageNumber"]= @sprintf("%d", page)
        r = get(rest_uri*content_endpoint, query=query, timeout=timeout)

        if r.status != 200
            throw(BadResponseException(r.status))
        end

        json_response = Requests.json(r)
        # println("No Results ", length(json_response["result"]["results"]))
        if json_response["result"]["results"][1]["ui"]=="NONE"
            break
        end
        push!(result_pages,json_response)
    end

    return result_pages

end

"""
    best_match_cui(result_pages)

Retrive the best match from array of all result pages

###Example

```julia
cui = BioMedQuery.UMLS.best_match_cui(all_results)
```
"""
function best_match_cui(result_pages)
    return result_pages[1]["result"]["results"][1]["ui"]
end

"""
    get_semantic_type(c::Credentials, cui)

Return an array of the semantic types associated with a cui

###Example

```julia
credentials = Credentials(user, psswd)
tgt = get_tgt(credentials)
cui = "C0028754"
sm = BioMedQuery.UMLS.get_semantic_type(tgt, cui)
```
"""
function get_semantic_type(tgt, cui)

    content_endpoint = "/rest/content/current/CUI/"*cui
    #get a new ticket
    ticket = get_ticket(tgt)
    r = get( rest_uri*content_endpoint,query=Dict("ticket"=> ticket))

    json_response = Requests.json(r)
    st = json_response["result"]["semanticTypes"]
    concepts = Array{ASCIIString}(length(st))
    for (ci, concept) in enumerate(st)
        concepts[ci] =concept["name"]
    end

    return concepts
end


text:
function write_ine(str::String, A::Matrix{Int}, b::Vector{Int})
    fp = open(str, "w")
    println(fp, "H-representation")
    println(fp, "begin")
    m, n = size(A)
    @assert m == length(b)
    println(fp, "$m $(n+1) integer")
    for i in 1:m
        println(fp, join([b[i] A[i,:]], " "))
    end
    println(fp, "end")
    close(fp)
    nothing
end

function write_ine(str::String, A::Matrix{Float64}, b::Vector{Float64})
    fp = open(str, "w")
    println(fp, "H-representation")
    println(fp, "begin")
    m, n = size(A)
    @assert m == length(b)
    println(fp, "$m $(n+1) real")
    for i in 1:m
        @printf(fp, "%E ", b[i])
        for j in 1:n
            @printf(fp, "%E ", A[i,j])
        end
        println(fp)
    end
    println(fp, "end")
    close(fp)
    nothing
end

function write_ext(str::String, ext::Vector)
    fp = open(str, "w")
    println(fp, "V-representation")
    println(fp, "begin")
    m = length(ext)
    n = length(ext[1])
    println(fp, "$m $n real")
    for i in 1:m
        vv = ext[i]
        @assert length(vv) == n
        println(fp, join(vv, " "))
    end
    println(fp, "end")
    println(fp, "hull")
    close(fp)
    nothing
end

function read_ext(fname::String)
    fp = open(fname, "r")
    vertices = Vector[]
    rays = Vector[]

    while true
        line = strip(readline(fp), '\n')
        line == "begin" && break
    end
    line = strip(readline(fp), '\n')
    (m,n,typ) = filter(x->!isempty(x), split(line, " "))
    line = strip(readline(fp), '\n')
    elem = []
    while line != "end"
        res  = filter(x->!isempty(x), split(line, " "))
        nums = map(float, res)
        try
            elem = convert(Vector{Int}, nums)
        catch
            elem = nums
        end
        if elem[1] == 1 # vertex
            push!(vertices, elem[2:end])
        else
            push!(rays, elem[2:end])
        end
        line = strip(readline(fp), '\n')
    end
    close(fp)
    return vertices, rays
end

function read_ine(fname::String)
    fp = open(fname, "r")

    while true
        line = strip(readline(fp), '\n')
        line == "begin" && break
    end
    line = strip(readline(fp), '\n')
    (mm,nn,typ) = filter(x->!isempty(x), split(line, " "))
    m,n = int(mm), int(nn)
    line = strip(readline(fp), '\n')
    A = Array(Float64,0,n-1)
    b = Float64[]
    while line != "end"
        res  = filter(x->!isempty(x), split(line, " "))
        elem = map(float, res)
        push!(b, elem[1])
        A = [A; -elem[2:end]']
        line = strip(readline(fp), '\n')
    end
    @assert size(A) == (m,n-1)
    close(fp)
    return b, A
end


text:
# singleton for indicating if no default arguments are present
struct NoDefault end
const NO_DEFAULT = NoDefault()

"""
    @addlogprob!(ex)

Add the result of the evaluation of `ex` to the joint log probability.
"""
macro addlogprob!(ex)
    return quote
        acclogp!($(esc(:(__varinfo__))), $(esc(ex)))
    end
end

"""
    getargs_dottilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L .~ R` or
`(~).(L, R)`, or `nothing` otherwise.
"""
getargs_dottilde(x) = nothing
function getargs_dottilde(expr::Expr)
    return MacroTools.@match expr begin
        (.~)(L_, R_) => (L, R)
        (~).(L_, R_) => (L, R)
        # Julia 1.6: see https://github.com/TuringLang/Turing.jl/issues/1525
        (L_ .~ R_) => (L, R)
        x_ => nothing
    end
end

"""
    getargs_tilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L ~ R`, or `nothing`
otherwise.
"""
getargs_tilde(x) = nothing
function getargs_tilde(expr::Expr)
    return MacroTools.@match expr begin
        (~)(L_, R_) => (L, R)
        x_ => nothing
    end
end

function to_namedtuple_expr(syms, vals=syms)
    length(syms) == 0 && return :(NamedTuple())

    names_expr = Expr(:tuple, QuoteNode.(syms)...)
    vals_expr = Expr(:tuple, vals...)
    return :(NamedTuple{$names_expr}($vals_expr))
end

#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = [r]
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
reconstruct(d::UnivariateDistribution, val::AbstractVector) = val[1]
reconstruct(d::MultivariateDistribution, val::AbstractVector) = copy(val)
function reconstruct(d::MatrixDistribution, val::AbstractVector)
    return reshape(copy(val), size(d))
end
function reconstruct!(r, d::Distribution, val::AbstractVector)
    return reconstruct!(r, d, val)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector)
    r .= val
    return r
end
function reconstruct(d::Distribution, val::AbstractVector, n::Int)
    return reconstruct(d, val, n)
end
function reconstruct(d::UnivariateDistribution, val::AbstractVector, n::Int)
    return copy(val)
end
function reconstruct(d::MultivariateDistribution, val::AbstractVector, n::Int)
    return copy(reshape(val, size(d)[1], n))
end
function reconstruct(d::MatrixDistribution, val::AbstractVector, n::Int)
    tmp = reshape(val, size(d)[1], size(d)[2], n)
    orig = [tmp[:, :, i] for i in 1:size(tmp, 3)]
    return orig
end
function reconstruct!(r, d::Distribution, val::AbstractVector, n::Int)
    return reconstruct!(r, d, val, n)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector, n::Int)
    r .= val
    return r
end

# Uniform random numbers with range 4 for robust initializations
# Reference: https://mc-stan.org/docs/2_19/reference-manual/initialization.html
randrealuni(rng::Random.AbstractRNG) = 4 * rand(rng) - 2
randrealuni(rng::Random.AbstractRNG, args...) = 4 .* rand(rng, args...) .- 2

const Transformable = Union{
    PositiveDistribution,
    UnitDistribution,
    TransformDistribution,
    SimplexDistribution,
    PDMatDistribution,
}
istransformable(dist) = false
istransformable(::Transformable) = true

#################################
# Single-sample initialisations #
#################################

inittrans(rng, dist::UnivariateDistribution) = invlink(dist, randrealuni(rng))
function inittrans(rng, dist::MultivariateDistribution)
    return invlink(dist, randrealuni(rng, size(dist)[1]))
end
inittrans(rng, dist::MatrixDistribution) = invlink(dist, randrealuni(rng, size(dist)...))

################################
# Multi-sample initialisations #
################################

inittrans(rng, dist::UnivariateDistribution, n::Int) = invlink(dist, randrealuni(rng, n))
function inittrans(rng, dist::MultivariateDistribution, n::Int)
    return invlink(dist, randrealuni(rng, size(dist)[1], n))
end
function inittrans(rng, dist::MatrixDistribution, n::Int)
    return invlink(dist, [randrealuni(rng, size(dist)...) for _ in 1:n])
end

#######################
# Convenience methods #
#######################
collectmaybe(x) = x
collectmaybe(x::Base.AbstractSet) = collect(x)


text:
@testset "Utils" begin

end


text:
export ruleSPPointMassConstraintOut

#-------------
# Update rules
#-------------

function ruleSPPointMassConstraintOut(msg_out::Message)
    m = unsafeMode(msg_out.dist)
    V = variateType(m)

    return Message(V, PointMass, m=m)
end

#--------------------------
# Custom inbounds collector
#--------------------------

function collectSumProductNodeInbounds(node::PointMassConstraint, entry::ScheduleEntry)
    inbounds = Any[]

    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if isa(inbound_interface.node, Clamp)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    return inbounds
end


text:
# Ideally we would want to use the Z3.jl interface,
# but it is currently broken.
# Ideally it looks something like this
#
# import Z3
# toZ3Type(::Type{Int64}) = Integer
# toZ3(s::Sym) = Z3.Var(toZ3Type(s._type), name = string(s.name))
# toZ3(c::T) where {T<:Real} = Z3.NumeralAst(toZ3Type(T), c)

# function symbolic(tCallsite)
#     stream = filter(t)
#     @assert isempty(t.stack)
#     @assert length(t.current) == 1
#     entry = t.current[1]

#     s = Z3.Solver()
#     vars = Dict{Sym, Any}()
#     map(entry.args) do sym
#         var = toZ3(sym)
#         vars[sym] = var
#     end

#     getZ3(s::Sym) = vars[s]
#     getZ3(c) = toZ3(c)

#     for (f, ret, args) in stream
#         z3args = map(getZ3, args)
#         if f == assert
#             @assert length(args) == 1
#             @show ret
#             Z3.solver_assert(s, z3args[1])
#         else
#             z3new = f(z3args...)
#             @assert ret isa Sym
#             vars[ret] = z3new
#         end
#     end
#     m = Z3.solver_get_model(s)
#     return s, m
# end

nbits(::Type{Int128}) = 128
nbits(::Type{Int64}) = 64
nbits(::Type{Int32}) = 32
nbits(::Type{Int16}) = 16
nbits(::Type{Int8}) = 8
nbits(::Type{UInt128}) = 128
nbits(::Type{UInt64}) = 64
nbits(::Type{UInt32}) = 32
nbits(::Type{UInt16}) = 16
nbits(::Type{UInt8}) = 8

bitsToT(::Val{128}) = UInt128
bitsToT(::Val{64}) = UInt64
bitsToT(::Val{32}) = UInt32
bitsToT(::Val{16}) = UInt16
bitsToT(::Val{8}) = UInt8
bitsToT(::Val{N}) where N = error("Can construct type with $N bits")

toZ3(::Type{T}) where T<:Integer = "(_ BitVec $(nbits(T)))"
toZ3(::Type{Bool})    = "Bool"

toZ3(x::Integer) = "(_ bv$x $(nbits(typeof(x))))"
toZ3(x::Bool) = string(x)
toZ3(x) = error("toZ3 for $x is not a thing yet")

FtoZ3(f::Function) where T = error("Can't handle $f yet")
FtoZ3(::typeof(Base.ifelse)) = "ite"
FtoZ3(::typeof(Base.:(===))) = "="
#FtoZ3(::typeof(Base.:(==)), ::Type{<:Integer}) = "bveq"
#FtoZ3(::typeof(Base.:(==)), ::Type{Bool}) = "="

function FtoZ3(f::Core.IntrinsicFunction)
    if f === Core.Intrinsics.sub_int
        return "bvsub"
    elseif f === Core.Intrinsics.add_int
        return "bvadd"
    elseif f === Core.Intrinsics.mul_int
        return "bvmul"
    elseif f === Core.Intrinsics.sdiv_int ||
           f === Core.Intrinsics.checked_sdiv_int
        return "bvudiv"
    elseif f === Core.Intrinsics.sdiv_int
        return "bvudiv"
    elseif f === Core.Intrinsics.udiv_int
        return "bvsdiv"
    elseif f === Core.Intrinsics.slt_int
        return "bvslt"
    elseif f === Core.Intrinsics.ult_int
        return "bvult"
    elseif f === Core.Intrinsics.sle_int
        return "bvsle"
    elseif f === Core.Intrinsics.ule_int
        return "bvule"
    elseif f === Core.Intrinsics.not_int
        return "bvnot"
    elseif f === Core.Intrinsics.and_int
        return "bvand"
    else
        error("Can't handle IntrinsicFunction $f yet")
    end
end

validName(name) = "|$(name)|"
declaration(s::Sym) = "(declare-const $(validName(s.name)) $(toZ3(s._type)))"

function symbolic(t::Callsite)
    return symbolic(filter(t))
end

Tunbox(x) = typeof(x)
Tunbox(s::Sym) = s._type

function symbolic(stream)
    declarations = IOBuffer()
    model = IOBuffer()

    function declare(s::Sym)
        decl = declaration(s)
        println(declarations, decl)
        return validName(s.name)
    end

    vars = Dict{Sym, String}()

    function getZ3(s::Sym)
        if !haskey(vars, s)
            vars[s] = declare(s)
        end
        vars[s]
    end
    getZ3(c) = toZ3(c)

    for (f, ret, args) in stream
        ret = anything(ret)
        z3args = join(map(getZ3, args), " ")
        if f == assert || f == prove
            @assert length(args) == 1
            @assert ret isa Bool
            z3ret = toZ3(ret)
            stmt = "(= $z3args $(z3ret))"
            if f == prove
                stmt = "(not $stmt)"
            end
            stmt = "(assert $stmt)"
        else
            @assert ret isa Sym
            z3f = FtoZ3(f)
            z3ret = getZ3(ret)
            stmt = "(assert (= $z3ret ($z3f $z3args)))"
        end
        println(model, stmt)
    end
    seekstart(model)
    write(declarations, model)
    seekstart(declarations)
    declarations
end

function runZ3(model)
    tinput = Pipe()
    output = Pipe()
    run(pipeline(`z3 -smt2 -in`, stdin=tinput, stdout=output), wait=false)
    input = IOBuffer()
    println(input, "(set-option :pp.bv-literals false)")
    write(input, model)
    write(input, "(check-sat)\n")
    write(input, "(get-model)\n")
    seekstart(input)
    @debug begin
        z3input = read(input, String)
        seekstart(input)
        """
        Z3 input:
        $z3input
        """
    end
    write(tinput, input)
    close(tinput)
    close(output.in)
    sat = readline(output)
    model = read(output, String)
    @debug """Z3 output
    $sat
    $model
    """
    return (sat == "sat", model)
end

function fromz3type(typ)
    if typ == "Int"
        return Int
    elseif typ == "Bool"
        return Bool
    else
        r_typ = r"\(_ BitVec (\d{1,3})\)"
        m = match(r_typ, typ)
        if m === nothing
            error("What even is $typ")
        else
            return bitsToT(Val(parse(Int, m.captures[1])))
        end
    end
end

function fromz3val(T, val)
    if T == Int
        return parse(Int, val)
    elseif T == Bool
        return parse(Bool, val)
    else
        return parse(T, split(val)[2][3:end])
    end
end

function stringToType(_type)
    if _type == "Int128"
        return Int128
    elseif _type == "Int64"
        return Int64
    elseif _type == "Int32"
        return Int32
    elseif _type == "Int16"
        return Int16
    elseif _type == "Int8"
        return Int8
    elseif _type == "UInt128"
        return UInt128
    elseif _type == "UInt64"
        return UInt64
    elseif _type == "UInt32"
        return UInt32
    elseif _type == "UInt16"
        return UInt16
    elseif _type == "UInt8"
        return UInt8
    elseif _type == "Bool"
        return Bool
    end
end

# This is really hacky, in reality we would want to properly parse
# s-expr.
function parseZ3(model)
    lines = split(chomp(model), '\n')[2:end-1]
    @assert length(lines) % 2 == 0
    inputs = Any[]

    r_def = r"define-fun\s+(\|.*\|)\s+\(\)\s+(\(.+?\)|Bool|Int)\s+(\(.+?\)|\S+)"
    for i in 1:2:length(lines)
        def = strip(join(lines[i:i+1]))[2:end-1]
        m = match(r_def, def)
        if m === nothing
            error("Regex didn't match: $def")
        end
        name = m.captures[1]
        T = fromz3type(m.captures[2])
        val = fromz3val(T, m.captures[3])
        push!(inputs, (name, val))
    end

    args = Any[]
    rands = Any[]
    others = Any[]

    r_arg = r"\|##arg_(\d+)#(\w+)#\d+\|"
    r_rand = r"\|##rand#(\w+)#(\d+)\|"
    r_others = r"\|##(\w+)#\d+\|"
    for (name, val) in inputs
        m = match(r_arg, name)
        if m !== nothing
            id = parse(Int, m.captures[1])
            T = stringToType(m.captures[2])
            if T <: Integer
                val = val % T
            end
            push!(args, (id, val))
            continue
        end
        m = match(r_rand, name)
        if m !== nothing
            id = parse(Int, m.captures[2])
            T = stringToType(m.captures[1])
            if T <: Integer
                val = val % T
            end
            push!(rands, (id, val))
            continue
        end
        m = match(r_others, name)
        if m !== nothing
            T = stringToType(m.captures[1])
            if T <: Integer
                val = val % T
            end
            push!(others, (name, val))
            continue
        end
        @error "Can't parse $name for $val"
    end
    args = Tuple(map(x -> x[2], sort(args))) # sort by id
    rands = map(x -> x[2], sort(rands)) # sort by id
    args, rands, others
end


text:
import ...timeevolution: dmaster_h!, dmaster_nh!, dmaster_h_dynamic!, check_master

"""
    stochastic.master(tspan, rho0, H, J, C; <keyword arguments>)

Time-evolution according to a stochastic master equation.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Deterministic part of the Hamiltonian.
* `J`: Vector containing all deterministic
        jump operators which can be of any arbitrary operator type.
* `C`: Vector containing the stochastic operators for a superoperator
        of the form `C[i]*rho + rho*Cdagger[i]`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::T, H::AbstractOperator{B,B},
                J, C;
                rates=nothing,
                Jdagger=dagger.(J), Cdagger=dagger.(C),
                fout=nothing,
                kwargs...) where {B,T<:Operator{B,B}}

    tmp = copy(rho0)

    n = length(C)

    dmaster_stoch(dx, t, rho, drho, n) = dmaster_stochastic(dx, rho, C, Cdagger, drho, n)

    isreducible = check_master(rho0, H, J, Jdagger, rates) && check_master_stoch(rho0, C, Cdagger)
    if !isreducible
        dmaster_h_determ(t, rho, drho) =
            dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp)
        integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    else
        Hnh = copy(H)
        if isa(rates, AbstractMatrix)
            for i=1:length(J), j=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im*rates[i,j])*Jdagger[i]*J[j]
            end
        elseif isa(rates, AbstractVector)
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im)*Jdagger[i]*J[i]
            end
        end
        Hnhdagger = dagger(Hnh)

        dmaster_nh_determ(t, rho, drho) =
            dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp)
        integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    end
end
master(tspan, psi0::Ket, args...; kwargs...) = master(tspan, dm(psi0), args...; kwargs...)

"""
    stochastic.master_dynamic(tspan, rho0, fdeterm, fstoch; <keyword arguments>)

Time-evolution according to a stochastic master equation with a
dynamic Hamiltonian and J.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `fdeterm`: Function `f(t, rho) -> (H, J, Jdagger)` or
        `f(t, rho) -> (H, J, Jdagger, rates)` giving the deterministic
        part of the master equation.
* `fstoch`: Function `f(t, rho) -> (C, Cdagger)` giving the stochastic
        superoperator of the form `C[i]*rho + rho*Cdagger[i]`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch` and all optional functions. If unset, the number
        is calculated automatically from the function outputs. NOTE: Set this
        number if you want to avoid an initial calculation of function outputs!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, rho0::T, fdeterm, fstoch;
                rates=nothing,
                fout=nothing,
                noise_processes::Int=0,
                kwargs...) where {B,T<:Operator{B,B}}

    tmp = copy(rho0)

    if noise_processes == 0
        fs_out = fstoch(0, rho0)
        n = length(fs_out[1])
    else
        n = noise_processes
    end

    dmaster_determ(t, rho, drho) = dmaster_h_dynamic!(drho, fdeterm, rates, rho, tmp, t)
    dmaster_stoch(dx, t, rho, drho, n) = dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n)
    integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; kwargs...)
end
master_dynamic(tspan, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...)

# Derivative functions
function dmaster_stochastic(dx::AbstractVector, rho, C, Cdagger, drho, n)
    recast!(drho,dx)
    QuantumOpticsBase.mul!(drho,C[1],rho)
    QuantumOpticsBase.mul!(drho,rho,Cdagger[1],true,true)
    drho.data .-= tr(drho)*rho.data
end
function dmaster_stochastic(dx::AbstractMatrix, rho, C, Cdagger, drho, n)
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(drho,dx_i)
        QuantumOpticsBase.mul!(drho,C[i],rho)
        QuantumOpticsBase.mul!(drho,rho,Cdagger[i],true,true)
        drho.data .-= tr(drho)*rho.data
        recast!(dx_i,drho)
    end
end

function dmaster_stoch_dynamic(dx, t, rho, f, drho, n)
    result = f(t, rho)
    QO_CHECKS[] && @assert 2 == length(result)
    C, Cdagger = result
    QO_CHECKS[] && check_master_stoch(rho, C, Cdagger)
    dmaster_stochastic(dx, rho, C, Cdagger, drho, n)
end

function integrate_master_stoch(tspan, df, dg,
                        rho0, fout,
                        n;
                        kwargs...)
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = as_vector(rho0)
    state = copy(rho0)
    dstate = copy(rho0)
    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function check_master_stoch(rho0::Operator{B,B}, C, Cdagger) where B
    # TODO: replace type checks by dispatch; make types of C known
    @assert length(C) == length(Cdagger)
    isreducible = true
    for c=C
        @assert isa(c, AbstractOperator{B,B})
        if !isa(c, DataOperator)
            isreducible = false
        end
    end
    for c=Cdagger
        @assert isa(c, AbstractOperator{B,B})
        if !isa(c, DataOperator)
            isreducible = false
        end
    end
    isreducible
end

as_vector(rho::Operator) = reshape(rho.data, length(rho))

# TODO: Speed up by recasting to n-d arrays, remove vector methods
function recast!(rho::Operator{B,B,T},x::Union{Vector, SubArray}) where {B,T}
    rho.data = reshape(x, size(rho.data))
end
recast!(x::SubArray,state::Operator{B,B}) where B = (x[:] = state.data)
recast!(x::Vector,state::Operator{B,B}) where B = nothing


text:
include("equality_test.jl")
using .Equality

include("optimization_tests.jl")
using .Optimization


text:
Random.seed!(12345)

#model = OSQP.Model()
P = [1.0 -1; -1 1]
q = [1.0, 1.0]
d = [0.0, 0]
bQP = BoxConstrainedQP(P, q, d)

sol = solve!(bQP)

@test sol == [0.0, 0.0]

P = [1.0 -1; -1 1]
q = [1.0-2.0,  1.0+2]
d = [0.0, 0]
bQP = BoxConstrainedQP(P, q, d)

sol = solve!(bQP)

@test sol == [1.0, 0.0]

P = [1.0 0; 0 0]
q = [-1.0,  0]
d = [0.0, 0]
bQP = BoxConstrainedQP(P, q, d)

sol = solve!(bQP)

@test sol[1] == 1.0
@test sol[2] >= 0

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
using BenchmarkTools, Compat, EmpiricalModeDecomposition

x = rand(1000)
t1000   = @benchmark fs = EmpiricalModeDecomposition.IMF($x)
@show t1000
println("Performance on rand(1000) $t1000")

x = rand(10000)
t10000  = @benchmark fs = EmpiricalModeDecomposition.IMF($x)
@show t10000
println("Performance on rand(10000) $t10000")

x = rand(100000)
t100000 = @benchmark fs = EmpiricalModeDecomposition.IMF($x)
println("Performance on rand(100000) $t100000")

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition