text:
"""
`Manopt.jl` – Optimization on Manifolds in Julia.
"""
module Manopt
using Colors
using ColorSchemes
using ColorTypes
using Markdown
using LinearAlgebra: I, Diagonal, eigvals, eigen, tril
using Dates: Period, Nanosecond, value
using Requires
using Random: shuffle!
using DataStructures: CircularBuffer, capacity, length, size, push!
using StaticArrays
import Base: copy, identity, &, |
import ManifoldsBase:
    ℝ,
    ℂ,
    ×,
    ^,
    _read,
    _write,
    AbstractBasis,
    AbstractPowerManifold,
    AbstractVectorTransportMethod,
    AbstractRetractionMethod,
    AbstractInverseRetractionMethod,
    CachedBasis,
    DefaultOrthonormalBasis,
    ExponentialRetraction,
    LogarithmicInverseRetraction,
    ParallelTransport,
    PowerManifold,
    AbstractManifold,
    allocate,
    allocate_result,
    allocate_result_type,
    copyto!,
    default_inverse_retraction_method,
    default_retraction_method,
    default_vector_transport_method,
    distance,
    exp,
    exp!,
    log,
    log!,
    injectivity_radius,
    inner,
    geodesic,
    get_basis,
    get_component,
    get_coordinates,
    get_vector,
    get_vectors,
    get_iterator,
    getindex,
    manifold_dimension,
    mid_point,
    mid_point!,
    NestedPowerRepresentation,
    norm,
    power_dimensions,
    project,
    project!,
    retract,
    retract!,
    inverse_retract,
    inverse_retract!,
    shortest_geodesic,
    vector_transport_to,
    vector_transport_to!,
    zero_vector,
    zero_vector!,
    DiagonalizingOrthonormalBasis,
    representation_size,
    setindex!,
    set_component!

include("plans/plan.jl")
# Functions
include("functions/bezier_curves.jl")
include("functions/adjoint_differentials.jl")
include("functions/costs.jl")
include("functions/differentials.jl")
include("functions/gradients.jl")
include("functions/Jacobi_fields.jl")
include("functions/proximal_maps.jl")
# solvers general framework
include("solvers/solver.jl")
# specific solvers
include("solvers/ChambollePock.jl")
include("solvers/conjugate_gradient_descent.jl")
include("solvers/cyclic_proximal_point.jl")
include("solvers/DouglasRachford.jl")
include("solvers/NelderMead.jl")
include("solvers/gradient_descent.jl")
include("solvers/particle_swarm.jl")
include("solvers/quasi_Newton.jl")
include("solvers/truncated_conjugate_gradient_descent.jl")
include("solvers/trust_regions.jl")
include("solvers/stochastic_gradient_descent.jl")
include("solvers/subgradient.jl")
include("solvers/debug_solver.jl")
include("solvers/record_solver.jl")
include("helpers/errorMeasures.jl")
include("helpers/exports/Asymptote.jl")
include("data/artificialDataFunctions.jl")

function __init__()
    @require Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" begin
        using .Manifolds:
            AbstractGroupManifold,
            Circle,
            Euclidean,
            Grassmann,
            Hyperbolic,
            PositiveNumbers,
            ProductManifold,
            Rotations,
            SymmetricPositiveDefinite,
            Stiefel,
            Sphere,
            TangentBundle,
            FixedRankMatrices,
            SVDMPoint,
            UMVTVector,
            ArrayPowerRepresentation,
            ProductRepr,
            submanifold_components,
            sym_rem,
            mean
        import Random: rand, randperm
        using LinearAlgebra: cholesky, det, diag, dot, Hermitian, qr, Symmetric, triu
        include("helpers/random.jl")
        # adaptions for Nonmutating manifolds
        const NONMUTATINGMANIFOLDS = Union{Circle,PositiveNumbers,Euclidean{Tuple{}}}
        include("functions/manifold_functions.jl")
        include("functions/nonmutating_manifolds_functions.jl")
        include("plans/nonmutating_manifolds_plans.jl")
        include("plans/alternating_gradient_plan.jl")
        include("solvers/alternating_gradient_descent.jl")
        export random_point, random_tangent, mid_point, mid_point!, reflect, reflect!
        export AlternatingGradientDescentOptions, AlternatingGradientProblem
        export AlternatingGradient
        export alternating_gradient_descent, alternating_gradient_descent!
    end
    return nothing
end
#
# General
export ℝ, ℂ, &, |
#
# Problems
export Problem,
    ProximalProblem,
    CostProblem,
    SubGradientProblem,
    GradientProblem,
    HessianProblem,
    PrimalDualProblem,
    StochasticGradientProblem,
    AbstractEvaluationType,
    AllocatingEvaluation,
    MutatingEvaluation
#
# Options
export Options,
    AbstractGradientOptions,
    ChambollePockOptions,
    ConjugateGradientDescentOptions,
    CyclicProximalPointOptions,
    DouglasRachfordOptions,
    GradientDescentOptions,
    AbstractHessianOptions,
    NelderMeadOptions,
    ParticleSwarmOptions,
    PrimalDualOptions,
    RecordOptions,
    StochasticGradientDescentOptions,
    SubGradientMethodOptions,
    TruncatedConjugateGradientOptions,
    TrustRegionsOptions
#
# Accessors and helpers for Options
export linesearch_backtrack
export get_cost,
    get_gradient,
    get_gradient!,
    get_subgradient,
    get_subgradient!,
    get_proximal_map,
    get_proximal_map!,
    get_options,
    get_initial_stepsize,
    get_gradients,
    get_gradients!,
    get_primal_prox,
    get_primal_prox!,
    get_dual_prox,
    get_dual_prox!,
    linearized_forward_operator,
    linearized_forward_operator!,
    adjoint_linearized_operator,
    adjoint_linearized_operator!,
    forward_operator,
    forward_operator!
export get_hessian, get_hessian!, ApproxHessianFiniteDifference
export is_options_decorator, dispatch_options_decorator
export primal_residual, dual_residual

export QuasiNewtonOptions, QuasiNewtonLimitedMemoryDirectionUpdate
export QuasiNewtonCautiousDirectionUpdate,
    BFGS, InverseBFGS, DFP, InverseDFP, SR1, InverseSR1
export InverseBroyden, Broyden
export AbstractQuasiNewtonDirectionUpdate, AbstractQuasiNewtonUpdateRule
export WolfePowellLineseach,
    StrongWolfePowellLineseach,
    operator_to_matrix,
    square_matrix_vector_product,
    WolfePowellBinaryLinesearch

export ConjugateGradientDescentOptions,
    GradientDescentOptions,
    AbstractHessianOptions,
    SubGradientMethodOptions,
    NelderMeadOptions,
    TruncatedConjugateGradientOptions,
    TrustRegionsOptions,
    ParticleSwarmOptions
export AbstractOptionsAction, StoreOptionsAction
export has_storage, get_storage, update_storage!

#
# Direction Update Rules
export DirectionUpdateRule,
    IdentityUpdateRule, StochasticGradient, AverageGradient, MomentumGradient, Nesterov
export DirectionUpdateRule,
    SteepestDirectionUpdateRule,
    HeestenesStiefelCoefficient,
    FletcherReevesCoefficient,
    PolakRibiereCoefficient,
    ConjugateDescentCoefficient,
    LiuStoreyCoefficient,
    DaiYuanCoefficient,
    HagerZhangCoefficient
#
# Solvers
export ChambollePock,
    ChambollePock!,
    conjugate_gradient_descent,
    conjugate_gradient_descent!,
    cyclic_proximal_point,
    cyclic_proximal_point!,
    DouglasRachford,
    DouglasRachford!,
    gradient_descent,
    gradient_descent!,
    NelderMead,
    NelderMead!,
    particle_swarm,
    particle_swarm!,
    quasi_Newton,
    quasi_Newton!,
    stochastic_gradient_descent,
    stochastic_gradient_descent!,
    subgradient_method,
    subgradient_method!,
    truncated_conjugate_gradient_descent,
    truncated_conjugate_gradient_descent!,
    trust_regions,
    trust_regions!
# Solver helpers
export decorate_options
export initialize_solver!, step_solver!, get_solver_result, stop_solver!
export solve
export ApproxHessianFiniteDifference
#
# Stepsize
export ConstantStepsize, DecreasingStepsize
export Linesearch, ArmijoLinesearch, NonmonotoneLinesearch
export get_stepsize, get_initial_stepsize, get_last_stepsize
#
# Stopping Criteria
export StopIfResidualIsReducedByFactor,
    StopIfResidualIsReducedByPower,
    StopWhenCurvatureIsNegative,
    StopWhenTrustRegionIsExceeded,
    StopWhenModelIncreased
export StopAfterIteration, StopWhenChangeLess, StopWhenGradientNormLess, StopWhenCostLess
export StopAfter, StopWhenAll, StopWhenAny
export get_active_stopping_criteria, get_stopping_criteria, get_reason
export are_these_stopping_critera_active
export StoppingCriterion, StoppingCriterionSet, Stepsize
#
# Data functions
export artificial_S1_signal, artificial_S1_slope_signal, artificialIn_SAR_image
export artificial_SPD_image, artificial_SPD_image2
export artificial_S2_whirl_image, artificial_S2_whirl_patch
export artificial_S2_rotation_image
export artificial_S2_whirl_patch, artificial_S2_lemniscate
export artificial_S2_composite_bezier_curve
#
# Exports
export asymptote_export_S2_signals, asymptote_export_S2_data, asymptote_export_SPD
export render_asymptote
#
# Coeffs & Helpers for differentials
export βdifferential_geodesic_startpoint, βdifferential_exp_basepoint
export βdifferential_exp_argument, βdifferential_log_basepoint, βdifferential_log_argument
export jacobi_field, adjoint_Jacobi_field
#
# Adjoint differentials
export adjoint_differential_geodesic_startpoint, adjoint_differential_geodesic_startpoint!
export adjoint_differential_geodesic_endpoint, adjoint_differential_geodesic_endpoint!
export adjoint_differential_exp_basepoint, adjoint_differential_exp_basepoint!
export adjoint_differential_exp_argument, adjoint_differential_exp_argument!
export adjoint_differential_log_basepoint, adjoint_differential_log_basepoint!
export adjoint_differential_log_argument, adjoint_differential_log_argument!
export adjoint_differential_forward_logs, adjoint_differential_forward_logs!
export adjoint_differential_bezier_control, adjoint_differential_bezier_control!
#
# Differentials
export differential_geodesic_startpoint, differential_geodesic_startpoint!
export differential_geodesic_endpoint, differential_geodesic_endpoint!
export differential_exp_basepoint, differential_exp_basepoint!
export differential_exp_argument, differential_exp_argument!
export differential_log_basepoint, differential_log_basepoint!
export differential_log_argument, differential_log_argument!
export differential_forward_logs, differential_forward_logs!
export differential_bezier_control, differential_bezier_control!
#
# Functions
export costL2TV, costL2TVTV2, costL2TV2, costTV, costTV2, costIntrICTV12
export cost_L2_acceleration_bezier, cost_acceleration_bezier
# Gradients
export grad_TV,
    grad_TV!,
    grad_TV2,
    grad_TV2!,
    grad_intrinsic_infimal_convolution_TV12,
    forward_logs,
    forward_logs!,
    grad_distance,
    grad_distance!,
    grad_acceleration_bezier,
    grad_L2_acceleration_bezier
# Proximal maps
export prox_distance, prox_distance!
export prox_TV, prox_TV!
export prox_parallel_TV, prox_parallel_TV!
export prox_TV2, prox_TV2!
export project_collaborative_TV, project_collaborative_TV!
# Error measures
export meanSquaredError, meanAverageError
#
# Bézier
export BezierSegment,
    de_casteljau,
    get_bezier_degrees,
    get_bezier_degree,
    get_bezier_inner_points,
    get_bezier_junction_tangent_vectors,
    get_bezier_junctions,
    get_bezier_points,
    get_bezier_segments
#
# Debugs
export DebugOptions, DebugAction, DebugGroup, DebugEntry, DebugEntryChange, DebugEvery
export DebugChange, DebugIterate, DebugIteration, DebugDivider
export DebugCost, DebugStoppingCriterion, DebugFactory, DebugActionFactory
export DebugGradient, DebugGradientNorm, DebugStepsize
export DebugPrimalBaseChange, DebugPrimalBaseIterate, DebugPrimalChange, DebugPrimalIterate
export DebugDualBaseChange, DebugDualBaseIterate, DebugDualChange, DebugDualIterate
export DebugDualResidual, DebugPrimalDualResidual, DebugPrimalResidual
export DebugProximalParameter
export DebugGradient, DebugGradientNorm, DebugStepsize
#
# Records - and access functions
export get_record, get_record_options, get_record_action, has_record
export RecordAction
export RecordActionFactory, RecordFactory
export RecordGroup, RecordEvery
export RecordChange, RecordCost, RecordIterate, RecordIteration
export RecordEntry, RecordEntryChange
export RecordGradient, RecordGradientNorm, RecordStepsize
export RecordPrimalBaseChange,
    RecordPrimalBaseIterate, RecordPrimalChange, RecordPrimalIterate
export RecordDualBaseChange, RecordDualBaseIterate, RecordDualChange, RecordDualIterate
export RecordProximalParameter
end


text:
module Data
using DataFrames
using Query
using Dates
using StatsPlots

export reformat, select_country, select_mainland, select_province, joindates, parsedates!, caseplot

function reformat(df::DataFrame, variable_name=:Date, value_name=:Count, cols=5)
    df = stack(df,names(df)[cols:end],
               variable_name=variable_name,
               value_name=value_name)
end

"""    select_province(df, country, province)

selects the rows of a dataframe based on their country and province
"""
function select_province(df, country, province)
    df = @from i in df begin
        @where i.Country == country && i.Province == province
        @select i
        @collect DataFrame
    end
end

"""    select_country(df, country, province)

selects the rows of a dataframe based on their country
"""
function select_country(df, country)
    df = @from i in df begin
        @where i.Country == country
        @select i
        @collect DataFrame
    end
end

"""    select_mainland(df, country, province)

selects the rows of a dataframe based on their country and excludes overseas territories
"""
function select_mainland(df, country)
    df = @from i in df begin
        @where i.Country == country && isna(i.Province)
        @select i
        @collect DataFrame
    end
end

"""     joindates(casesf, deaths, recovf)

join the cases, deaths, and recovery frames into one dataframe based on the dates and country+province
"""
function joindates(casesf, deathf, recovf)
    cdf = casesf |> @join(deathf, {_.Date, _.Province, _.Country},
    {_.Date, _.Province, _.Country},
    {_.Date, _.Province, _.Country, _.Cases, __.Deaths}) |> DataFrame

    ff = cdf |> @join(recovf, {_.Date, _.Province, _.Country},
    {_.Date, _.Province, _.Country},
    {_.Date, _.Province, _.Country, _.Cases, _.Deaths, __.Recoveries}) |> DataFrame
    return ff
end

"""     parsedates!(df::DataFrame)

replace the Dates column of a dataframe as string or symbol with the parsed Dates.Date version
"""
function parsedates!(df::DataFrame)
    df.Date = Date.(string.(df.Date), dateformat"m/d/y")
    return df
end

"""    caseplot(df)

plot the number of cases, deaths, recoveries over time
"""
function caseplot(df)
    @df df plot(:Date, [:Cases, :Deaths, :Recoveries], label=["cases" "deaths" "recoveries"])
end
end

ext: jl


text:
#
#  CollagenPolymer.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 30/03/2020.
#
#

module CollagenPolymer

# Import Julia packages
using Distributions
using LinearAlgebra
using Random
using StaticArrays

# Import local program modules
using InternalForces
using BoundaryForces
using InterMonomerForces
using CalculateNoise
using UpdateSystem
using OutputData
using Initialise
using CreateRunDirectory
using AdaptTimestep
using CellListFunctions
using Visualise

# Define function for bringing together modules to run simulation
@views function collagenPolymer(nMonomers,L,σ,ϵLJ_in,k_in,Ebend_in,boxSize,tMax,outputFlag,renderFlag)

    # Run parameters
    # nMonomers      (eg. 2    )   Number of collagen monomers
    # L              (eg. 0.5  )   Length of one monomer
    # σ              (eg. 0.05 )   Diameter of monomer = Diameter of particle within monomer/External LJ length scale (separation at which V=0) = 2*particle radius
    # ϵLJ_in         (eg. 00.0 )   External Lennard-Jones energy
    # k_in           (eg. 100.0)   Internal spring stiffness for Fraenkel spring forces between adjacent particles within a monomer
    # Ebend_in       (eg. 0.0  )   Internal bending modulus of monomer
    # boxSize        (eg. 1.0  )   Dimensions of cube in which particles are initialised
    # tMax           (eg. 0.1  )   Total simulation time
    # outputFlag     (eg. 1    )   Controls whether or not data is printed to file
    # renderFlag     (eg. 1    )   Controls whether or not system is visualised with povRay automatically

    # Thermodynamic parameters
    μ              = 1.0                                # Fluid viscosity
    kT             = 1.0                                # Boltzmann constant*Temperature

    # Derived parameters
    ϵLJ            = ϵLJ_in*kT                          # External Lennard-Jones energy
    k              = k_in*kT                            # Internal spring stiffness for forces between adjacent particles within a monomer
    Ebend          = Ebend_in*kT                        # Internal bending modulus of monomer
    outputInterval = tMax/100.0                         # Time interval for writing position data to file
    intrctnThrshld = 2.0*σ                              # Threshold separation for van der Waals interactions
    nDomains       = 1+ceil(Int64,(5.0*L)/(4.0*σ)+1.0)  # Number of particles per monomer, ensuring rₑ<=0.4σ and thus preventing corrugation issues
    rₑ             = (L-σ)/(nDomains-1)                 # Equilibrium separation of springs connecting adjacent particles within a monomer
    nParticles     = nMonomers*nDomains                 # Total number of particles in system
    rₘ             = σ*2.0^(1.0/6.0)                    # Separation at which LJ potential is minimised
    WCAthreshSq    = (rₘ)^2                             # Cutoff threshold for WCA potential - separation at which force is zero - squared
    D              = kT/(6.0*π*μ*σ)                     # Diffusion constant
    nGrid          = ceil(Int64,boxSize/intrctnThrshld) # Dimensions of cellLists array, for boundary interaction.

    # Initialise system time
    t  = 1E-10

    # Allocate system arrays and initialise monomers locations and orientations
    pos,F,W,magsF,magsW,pairsList,boundaryList,dxMatrix,neighbourCells,threadRNG,AA,BB,CC,DD,EE = initialise(nMonomers,nDomains,nParticles,rₑ,boxSize)

    # Setup data folder and output files; output initial state
    if outputFlag == 1
        foldername = createRunDirectory(nMonomers,L,nDomains,μ,kT,ϵLJ,σ,k,rₑ,Ebend,D,tMax,outputInterval,nGrid,boxSize,intrctnThrshld)
        outfile = open("output/$(foldername)/output.txt","w")
        outputData(pos,outfile)#,t,tMax)#,nMonomers,nDomains,σ)
    end

    println("|                   |                   |                   |                   |                  |")
    #Iterate over time until max system time is reached
    while t<tMax

        # Create list of particle interaction pairs based on cell lists algorithm
        pairsList,boundaryList = find_pairs(nParticles,pos,intrctnThrshld,nGrid,neighbourCells)

        # Calculate tension and bending forces within each monomer
        internalForces!(pos,F,nMonomers,nDomains,nParticles,k,rₑ,Ebend,AA,BB,CC,DD,EE)

        # Calculate van der Waals/electrostatic interactions between nearby monomer domains
        interMonomerForces!(pairsList,pos,F,nDomains,ϵLJ,σ,AA,WCAthreshSq,intrctnThrshld)

        # Calculate forces on particles from system boundary
        boundaryForces!(boundaryList,pos,F,ϵLJ,nGrid,boxSize,dxMatrix,rₘ)

        # Find stochastic term (Wiener process) for all monomers
        calculateNoise!(W,nParticles,threadRNG)

        # Adapt timestep to maximum force value
        Δt = adaptTimestep(F,W,magsF,magsW,σ)

        # Integrate system with forward euler
        t = updateSystem!(pos,F,W,t,Δt,D,kT,nParticles)

        if (t%outputInterval)<Δt
            print(">")
            outputFlag == 1 ? outputData(pos,outfile) : nothing
        end

    end
    print("\n")

    if outputFlag == 1
        close(outfile)
        renderFlag == 1 ? run(`python3 src/visualise.py output/$foldername`) : nothing
    end

end

export collagenPolymer

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
type XPCallbacks
    start_cdata::Function
    end_cdata::Function
    comment::Function
    character_data::Function
    default::Function
    default_expand::Function
    start_element::Function
    end_element::Function
    start_namespace::Function
    end_namespace::Function

    # Create an XPCallbacks instance filled with nop callbacks
    function XPCallbacks()
        start_cdata = (handler::XPStreamHandler) -> nothing
        end_cdata = (handler::XPStreamHandler) -> nothing
        comment = (handler::XPStreamHandler, txt::String) -> nothing
        character_data = (handler::XPStreamHandler, txt::String) -> nothing
        default = (handler::XPStreamHandler, txt::String) -> nothing
        default_expand = (handler::XPStreamHandler, txt::String) -> nothing
        start_element = (handler::XPStreamHandler, name::String, attrs_in::Dict{String,String}) -> nothing
        end_element = (handler::XPStreamHandler, name::String) -> nothing
        start_namespace = (handler::XPStreamHandler, prefix::String, uri::String) -> nothing
        end_namespace = (handler::XPStreamHandler, prefix::String) -> nothing

        new(start_cdata, end_cdata, comment, character_data, default,
            default_expand, start_element, end_element, start_namespace,
            end_namespace)
    end
end

type XPStreamHandler{D}
    cbs::XPCallbacks
    parser::XML_Parser
    data::D
end

function streaming_start_cdata (p_cbs::Ptr{Void})
    @DBG_PRINT ("Found StartCdata")
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    h.cbs.start_cdata(h)
    return
end
cb_streaming_start_cdata = cfunction(streaming_start_cdata, Void, (Ptr{Void},))

function streaming_end_cdata (p_cbs::Ptr{Void})
    @DBG_PRINT ("Found EndCdata")
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    h.cbs.end_cdata(h)
    return;
end
cb_streaming_end_cdata = cfunction(streaming_end_cdata, Void, (Ptr{Void},))

function streaming_cdata (p_cbs::Ptr{Void}, s::Ptr{Uint8}, len::Cint)
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    txt = bytestring(s, int(len))

    @DBG_PRINT ("Found CData : " * txt)
    h.cbs.character_data(h, txt)

    return;
end
cb_streaming_cdata = cfunction(streaming_cdata, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function streaming_start_element(p_cbs::Ptr{Void}, name::Ptr{Uint8}, attrs_in::Ptr{Ptr{Uint8}})
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler
    txt::String = bytestring(name)
    attrs::Dict{String,String} = attrs_in_to_dict(attrs_in)

    h.cbs.start_element(h, txt, attrs)

    return
end
cb_streaming_start_element = cfunction(streaming_start_element, Void, (Ptr{Void},Ptr{Uint8}, Ptr{Ptr{Uint8}}))

function streaming_end_element(p_h::Ptr{Void}, name::Ptr{Uint8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt::String = bytestring(name)
    @DBG_PRINT ("End element: $txt, current element: $(xph.pdata.name) ")

    h.cbs.end_element(h, txt)

    return
end
cb_streaming_end_element = cfunction(streaming_end_element, Void, (Ptr{Void},Ptr{Uint8}))

function streaming_comment (p_h::Ptr{Void}, data::Ptr{Uint8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = bytestring(data)
    @DBG_PRINT ("Found comment : " * txt)

    h.cbs.comment(h, txt)

    return
end
cb_streaming_comment = cfunction(streaming_comment, Void, (Ptr{Void},Ptr{Uint8}))

function streaming_default (p_h::Ptr{Void}, data::Ptr{Uint8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = bytestring(data)
    @DBG_PRINT ("Default : " * txt)

    h.cbs.default(h, txt)

    return;
end
cb_streaming_default = cfunction(streaming_default, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function streaming_default_expand (p_h::Ptr{Void}, data::Ptr{Uint8}, len::Cint)
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = bytestring(data)
    @DBG_PRINT ("Default Expand : " * txt)

    h.cbs.default_expand(h, txt)

    return;
end
cb_streaming_default_expand = cfunction(streaming_default_expand, Void, (Ptr{Void},Ptr{Uint8}, Cint))

function streaming_start_namespace (p_h::Ptr{Void}, prefix::Ptr{Uint8}, uri::Ptr{Uint8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    prefix = bytestring(prefix)
    uri = bytestring(uri)
    @DBG_PRINT ("start namespace prefix : $prefix, uri: $uri")

    h.cbs.start_namespace(h, prefix, uri)

    return;
end
cb_streaming_start_namespace = cfunction(streaming_start_namespace, Void, (Ptr{Void},Ptr{Uint8}, Ptr{Uint8}))

function streaming_end_namespace (p_h::Ptr{Void}, prefix::Ptr{Uint8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    prefix = bytestring(prefix)
    @DBG_PRINT ("end namespace prefix : $prefix")

    h.cbs.end_namespace(h, prefix)

    return;
end
cb_streaming_end_namespace = cfunction(streaming_end_namespace, Void, (Ptr{Void},Ptr{Uint8}))

# Unsupported callbacks: External Entity, NotationDecl, Not Stand Alone, Processing, UnparsedEntityDecl, StartDocType
# SetBase and GetBase

function make_parser(cbs::XPCallbacks,data=nothing,sep='\0')
    p::XML_Parser = (sep == '\0') ? XML_ParserCreate(C_NULL) : XML_ParserCreateNS(C_NULL, sep);
    if (p == C_NULL) error("XML_ParserCreate failed") end

    h = XPStreamHandler(cbs, p, data)
    p_h = pointer_from_objref(h)
    XML_SetUserData(p, p_h);

    XML_SetCdataSectionHandler(p, cb_streaming_start_cdata, cb_streaming_end_cdata)
    XML_SetCharacterDataHandler(p, cb_streaming_cdata)
    XML_SetCommentHandler(p, cb_streaming_comment)
    XML_SetDefaultHandler(p, cb_streaming_default)
    XML_SetDefaultHandlerExpand(p, cb_streaming_default_expand)
    XML_SetElementHandler(p, cb_streaming_start_element, cb_streaming_end_element)
#    XML_SetExternalEntityRefHandler(p, f_ExternaEntity)
    XML_SetNamespaceDeclHandler(p, cb_streaming_start_namespace, cb_streaming_end_namespace)
#    XML_SetNotationDeclHandler(p, f_NotationDecl)
#    XML_SetNotStandaloneHandler(p, f_NotStandalone)
#    XML_SetProcessingInstructionHandler(p, f_ProcessingInstruction)
#    XML_SetUnparsedEntityDeclHandler(p, f_UnparsedEntityDecl)
#    XML_SetStartDoctypeDeclHandler(p, f_StartDoctypeDecl)

    return h
end

function stop(h::XPStreamHandler)
    XML_StopParser(h.parser, XML_FALSE)
#    XML_ParserFree(h.parser)
end

function pause(h::XPStreamHandler)
    XML_StopParser(h.parser, XML_TRUE)
end

function resume(h::XPStreamHandler)
    XML_ResumeParser(h.parser)
end

function free(h::XPStreamHandler)
    XML_ParserFree(h.parser)
end

function parsefile(filename::String,callbacks::XPCallbacks; bufferlines=1024, data=nothing)
    h = make_parser(callbacks, data)
    # TODO: Support suspending for files too
    suspended = false
    file = open(filename, "r")
    try
        io = IOBuffer()
        while !eof(file)
            i::Int = 0
            truncate(io, 0)
            while i < bufferlines && !eof(file)
                write(io, readline(file))
                i += 1
            end
            txt = bytestring(io)
            rc = XML_Parse(h.parser, txt, length(txt.data), 0)
            if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
                # Do not fail if the user aborted the parsing
                error("Error parsing document : $rc")
            end
            if XML_GetErrorCode(h.parser) == XML_ERROR_ABORTED
                break
            end
        end
        rc = XML_Parse(h.parser, "", length(""), 1)
        #if (rc == XML_STATUS_SUSPENDED)
        #    suspended = true
        #    return XPStreamHandler(callbacks,  parser)
        #end
        if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
            # Do not fail if the user aborted the parsing
            error("Error parsing document : $rc")
        end
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(h.parser)
        @DBG_PRINT ("$e, $err, $line, $column, $pos")
        rethrow("$e, $err, $line, $column, $pos")
    finally
        if !suspended
            XML_ParserFree(h.parser)
        end
        close(file)
    end
end

function parse(txt::String,callbacks::XPCallbacks; data=nothing)
    h = make_parser(callbacks, data)
    suspended = false

    try
        rc = XML_Parse(h.parser, txt, length(txt.data), 1)
        if (rc == XML_STATUS_SUSPENDED)
            suspended = true
            return h
        end
        if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
            # Do not fail if the user aborted the parsing
            error("Error parsing document : $rc")
        end
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(h.parser)
        @DBG_PRINT ("$e, $err, $line, $column, $pos")
        rethrow("$e, $err, $line, $column, $pos")

    finally
        if !suspended
            XML_ParserFree(h.parser)
        end
    end
end


text:
#
#----------------------------------------------------------------------------#
export MeshNet
#--------------------------------------#
struct MeshNet
    nprt::Integer
    ndim::Integer

    pou::AbstractPOU
    bsp::AbstractBanachSpace
    C  ::AbstractArray
end
#--------------------------------------#
function MeshNet(pou::AbstractPOU
                ,bsp::AbstractBanachSpace)

    nprt = pou.nprt
    ndim = bsp.ndim

    init = Flux.glorot_uniform
    C = softmax(init(ndim,nprt);dims=1)

    return MeshNet(nprt,ndim,pou,bsp,C)
end
#--------------------------------------#
function MeshNet(nprt,ndim
                ,pouType::DataType=nnPOU
                ,bspType::DataType=Cosines)

    pou = pouType(nprt)
    bsp = bspType(nprt,ndim)

    return MeshNet(pou,bsp)
end
#--------------------------------------#
Flux.@functor MeshNet
DiffEqFlux.initial_params(net::MeshNet) = Flux.destructure(net)[1]
#--------------------------------------#
function (net::MeshNet)(x)  # [1,N]

    @unpack nprt,ndim = net
    @unpack pou,bsp,C = net

    N = size(x)[end]

    # partition of unity
    ϕ = pou(x)              # [nprt,N]
    ϕ = reshape(ϕ,1,nprt,N) # [1,nprt,N]

    # function space
    ψ = bsp(x)              # [ndim,nprt,N]

    u = @. ϕ * ψ * C        # [ndim,nprt,N]
    u = sum(u,dims=(1,2))   # [1,1,N]
    u = reshape(u,1,N)      # [1,N]

    return u
end
#----------------------------------------------------------------------------#
export MeshNet2D
#--------------------------------------#
struct MeshNet2D
    nprt::Integer
    ndim::Integer

    pou::AbstractPOU
    bsp::AbstractBanachSpace
    C  ::AbstractArray
end
#--------------------------------------#
function MeshNet2D(pou::AbstractPOU
                ,bsp::AbstractBanachSpace)

    nprt = pou.nprt
    ndim = bsp.ndim

    init = Flux.glorot_uniform
    C = softmax(init(ndim,nprt);dims=1)

    return MeshNet2D(nprt,ndim,pou,bsp,C)
end
#--------------------------------------#
function MeshNet2D(nprt,ndim
                ,pouType::DataType=nnPOU
                ,bspType::DataType=Cosines)

    pou = pouType(nprt)
    bsp = bspType(nprt,ndim)

    return MeshNet2D(pou,bsp)
end
#--------------------------------------#
Flux.@functor MeshNet2D
DiffEqFlux.initial_params(net::MeshNet2D) = Flux.destructure(net)[1]
#--------------------------------------#
function (net::MeshNet2D)(x,y)  # [1,N]

    @assert size(x) == size(y)

    @unpack nprt,ndim = net
    @unpack pou,bsp,C = net

    N = size(x)[end]

    # partition of unity
    ϕ = pou(x,y)             # [nprt,N]
    ϕ = reshape(ϕ,1,nprt,N) # [1,nprt,N]

    # function space
    ψ = bsp(x,y)            # [ndim,nprt,N]

    u = @. ϕ * ψ * C        # [ndim,nprt,N]
    u = sum(u,dims=(1,2))   # [1,1,N]
    u = reshape(u,1,N)      # [1,N]

    return u
end
#--------------------------------------#
function (net::MeshNet2D)(xy) # [2,N]
    x = xy[1,:]
    y = xy[2,:]
    return net(x,y)
end
#--------------------------------------#
export leastSq!
#--------------------------------------#
function leastSq!(net::MeshNet,x,u)

    @unpack nprt,ndim = net
    @unpack pou,bsp,C = net

    N = size(x)[end]

    # form g
    ϕ = pou(x)                  # [nprt,N]
    ϕ = reshape(ϕ,1,nprt,N)     # [1,nprt,N]
    ψ = bsp(x)                  # [ndim,nprt,N]
    g  = @. ϕ * ψ               # [ndim,nprt,N]

    # least square solve
    g  = reshape(g,ndim*nprt,N) # size [N,ndim*nprt]
    c  = lsq(u',g')
    C .= reshape(c,ndim,nprt)   # size [ndim,nprt]

    return
end
#--------------------------------------#
function leastSq!(net::MeshNet2D,x,y,u)

    @unpack nprt,ndim = net
    @unpack pou,bsp,C = net

    N = size(x)[end]

    # form g
    ϕ = pou(x,y)                # [nprt,N]
    ϕ = reshape(ϕ,1,nprt,N)     # [1,nprt,N]
    ψ = bsp(x,y)                # [ndim,nprt,N]
    g  = @. ϕ * ψ               # [ndim,nprt,N]

    # least square solve
    g  = reshape(g,ndim*nprt,N) # size [N,ndim*nprt]
    c  = lsq(u',g')
    C .= reshape(c,ndim,nprt)   # size [ndim,nprt]

    return
end
#----------------------------------------------------------------------------#


text:
using Forall
using Base.Test

let sum = 0
    res = forall(1:10) do i
        sum += i
    end
    @test sum === 55
    @test res === nothing
end

# This performs the same operations as above, except that it also
# allocates an array with the result of each iteration
let sum = 0
    res = map(1:10) do i
        sum += i
    end
    @test sum === 55
    @test typeof(res) === Vector{Int}
    @test length(res) == 10
    @test res[1] === 1
    @test res[end] === 55
end


text:
struct HydroPlan{T <: AbstractFloat}
    Q::Vector{T} # Discharge each hour
    S::Vector{T} # Spillage each hour
    M::Vector{T} # Reservoir content each hour
    H::Vector{T} # Resulting electricity production each hour

    function (::Type{HydroPlan})(Q::AbstractVector,S::AbstractVector,M::AbstractVector,H::AbstractVector)
        T = promote_type(eltype(Q),eltype(S),eltype(M),eltype(H),Float32)
        return new{T}(convert(Vector{T},Q),convert(Vector{T},S),convert(Vector{T},M),convert(Vector{T},H))
    end
end

struct HydroProductionPlan{T <: AbstractFloat}
    horizon::Horizon                        # Planning horizon
    modeldata::HydroPlantCollection{T,2}

    total_plan::HydroPlan{T}                     # Accumulation of individual plans
    individual_plans::Dict{Plant,HydroPlan{T}}   # Plans for each hydro power plant

    function (::Type{HydroProductionPlan})(horizon::Horizon,data::HydroPlantCollection{T,2},plan::HydroPlan{T},plans::Dict{Plant,HydroPlan{T}}) where T <: AbstractFloat
        return new{T}(horizon,data,plan,plans)
    end
end
Base.eltype(::HydroProductionPlan{T}) where T <: AbstractFloat = T

function HydroProductionPlan(model::AbstractHydroModel)
    optmodel = model.internalmodel
    horizon = HydroModels.horizon(model)
    data = hydrodata(model.data)
    plants = HydroModels.plants(model.indices)

    (haskey(optmodel.objDict,:Q) && haskey(optmodel.objDict,:S) && haskey(optmodel.objDict,:M) && haskey(optmodel.objDict,:H)) || error("Given model does not model hydro power production")

    Q = getvalue(optmodel.objDict[:Q])
    S = getvalue(optmodel.objDict[:S])
    M = getvalue(optmodel.objDict[:M])
    H = getvalue(optmodel.objDict[:H])

    individual_plans = Dict{Plant,HydroPlan{eltype(data)}}()

    for plant in plants
        individual_plans[plant] = HydroPlan([Q[plant,1,t] + Q[plant,2,t] for t in 1:nhours(horizon)],
                                            [S[plant,t] for t in 1:nhours(horizon)],
                                            [M[plant,t] for t in 1:nhours(horizon)],
                                            [data[plant].μ[1]*Q[plant,1,t] + data[plant].μ[2]*Q[plant,2,t] for t in 1:nhours(horizon)])
    end

    total_plan = HydroPlan(sum(p.Q for p in values(individual_plans)),
                           sum(p.S for p in values(individual_plans)),
                           sum(p.M for p in values(individual_plans)),
                           sum(p.H for p in values(individual_plans)))

    total_plan.H ≈ H.innerArray || error("Difference in total production: ", norm(total_plan.H-H.innerArray))

    return HydroProductionPlan(model.horizon,
                               data,
                               total_plan,
                               individual_plans)
end

# function HydroProductionPlan(horizon::Horizon,modeldata::HydroModelData,plants::Vector{Plant},model::JuMPModel,scenario::Int64)
#     @assert haskey(model.objDict,:Q) && haskey(model.objDict,:S) && haskey(model.objDict,:M) && haskey(model.objDict,:H) "Given JuMP model does not model hydro power production"

#     Q = getvalue(model.objDict[:Q])
#     S = getvalue(model.objDict[:S])
#     M = getvalue(model.objDict[:M])
#     H = getvalue(model.objDict[:H])
#     μ = modeldata.μ

#     individual_plans = Dict{Plant,Plan}()

#     for plant in plants
#         individual_plans[plant] = Plan([Q[scenario,plant,1,t] + Q[scenario,plant,2,t] for t in 1:hours(horizon)],
#                                        [S[scenario,plant,t] for t in 1:hours(horizon)],
#                                        [M[scenario,plant,t] for t in 1:hours(horizon)],
#                                        [μ[plant,1]*Q[scenario,plant,1,t] + μ[plant,2]*Q[scenario,plant,2,t] for t in 1:hours(horizon)])
#     end

#     total_plan = Plan(sum(p.Q for p in values(individual_plans)),
#                       sum(p.S for p in values(individual_plans)),
#                       sum(p.M for p in values(individual_plans)),
#                       sum(p.H for p in values(individual_plans)))

#     @assert total_plan.H ≈ H.innerArray[scenario,:] "Difference in total production"

#     return HydroProductionPlan(horizon,
#                           modeldata,
#                           individual_plans,
#                           total_plan)
# end

plants(plan::HydroProductionPlan) = collect(keys(plan.individual_plans))

discharge(plan::HydroProductionPlan) = plan.total_plan.Q
function discharge(plan::HydroProductionPlan,plant::Plant)
    if plant in plants(plan)
        return plan.individual_plans[plant].Q
    else
        throw(ArgumentError(string("Selected plant ",hour," not in production plan")))
    end
end

spillage(plan::HydroProductionPlan) = plan.total_plan.S
function spillage(plan::HydroProductionPlan,plant::Plant)
    if plant in plants(plan)
        return plan.individual_plans[plant].S
    else
        throw(ArgumentError(string("Selected plant ",hour," not in production plan")))
    end
end

reservoir(plan::HydroProductionPlan) = plan.total_plan.M
function reservoir(plan::HydroProductionPlan,plant::Plant)
    if plant in plants(plan)
        return plan.individual_plans[plant].M
    else
        throw(ArgumentError(string("Selected plant ",hour," not in production plan")))
    end
end

power(plan::HydroProductionPlan) = plan.total_plan.H
function power(plan::HydroProductionPlan,plant::Plant)
    if plant in plants(plan)
        return plan.individual_plans[plant].H
    else
        throw(ArgumentError(string("Selected plant ",hour," not in production plan")))
    end
end

function revenue(plan::HydroProductionPlan,λ::Vector{Float64})
    length(λ) >= nhours(plan.horizon) || error("Price vector is to short")
    return plan.H .* λ
end

function totalrevenue(plan::HydroProductionPlan,λ::Vector{Float64})
    length(λ) >= nhours(plan.horizon) || error("Price vector is to short")
    return plan.H⋅λ
end

function mean(plans::Vector{HydroProductionPlan{T}}) where T <: AbstractFloat
    individual_plans = Dict{Plant,HydroPlan{T}}()

    for plant in plants(plans[1])
        individual_plans[plant] = HydroPlan(mean([p.individual_plans[plant].Q for p in plans]),
                                            mean([p.individual_plans[plant].S for p in plans]),
                                            mean([p.individual_plans[plant].M for p in plans]),
                                            mean([p.individual_plans[plant].H for p in plans]))
    end

    total_plan = HydroPlan(sum(p.Q for p in values(individual_plans)),
                           sum(p.S for p in values(individual_plans)),
                           sum(p.M for p in values(individual_plans)),
                           sum(p.H for p in values(individual_plans)))

    return HydroProductionPlan(plans[1].horizon,
                               plans[1].modeldata,
                               individual_plans,
                               total_plan)
end

@recipe function f(plan::HydroProductionPlan; showHmax = false)
    Hmin = minimum(plan.total_plan.H)
    Hmax = showHmax ? sum([plan.modeldata.H̅[p] for p in plants(plan)]) : maximum(plan.total_plan.H)
    increment = mean(abs.(diff(plan.total_plan.H)))

    linewidth --> 4
    linecolor --> :black
    tickfontsize := 14
    tickfontfamily := "sans-serif"
    guidefontsize := 16
    guidefontfamily := "sans-serif"
    titlefontsize := 22
    titlefontfamily := "sans-serif"
    xlabel := "Hour"
    ylabel := "Energy Volume [MWh]"
    xlims := (1,nhours(plan.horizon))
    ylims := (Hmin-increment,Hmax+increment)
    xticks --> 1:1:nhours(plan.horizon)
    yticks := Hmin:increment:Hmax
    yformatter := (d) -> @sprintf("%.2f",d)

    @series begin
        label --> "H"
        1:1:nhours(plan.horizon),plan.total_plan.H
    end

    if showHmax
        @series begin
            label --> "Hmax"
            linestyle --> :dash
            linecolor --> :black
            linewidth --> 2
            0:1:nhours(plan.horizon),fill(Hmax,nhours(plan.horizon))
        end
    end
end

@recipe function f(plan::HydroProductionPlan,plant::Plant; showHmax = false)
    if !(plant in plants(plant))
        throw(ArgumentError(string("Selected plant ",hour," not in production plan")))
    end
    Hmin = minimum(plan.total_plan.H)
    Hmax = showHmax ? sum([plan.modeldata.H̅[p] for p in plants(plan)]) : maximum(plan.total_plan.H)
    increment = std(plan.total_plan.H)

    linewidth --> 4
    linecolor --> :black
    xlims := (-1,nhours(horizon))
    ylims := (Hmin-increment,Hmax+increment)
    xticks := 0:1:nhours(plan.horizon)
    yticks := Hmin:increment:Hmax

    @series begin
        plan.individual_plans[plant].H
    end

    if showHmax
        @series begin
            linestyle --> :dash
            linecolor --> :black
            linewidth --> 2
            fill(plan.modeldata.H̅[plant],nhours(plan.horizon))
        end
    end
end

function show(io::IO, ::MIME"text/plain", plan::HydroProductionPlan)
    show(io,plan)
end

function show(io::IO, plan::HydroProductionPlan)
    if get(io, :multiline, false)
        print(io,"Production Plan")
    else
        println(io,"Hydro Power Production Plan")
        println(io,"Power production:")
        Base.print_matrix(io,plan.total_plan.H)
    end
end


text:
jl_to_np_type(o::PyObject) = numpy.bool_
jl_to_np_type(o::Type{Bool}) = numpy.uint32
jl_to_np_type(o::Type{UInt32}) = numpy.uint32
jl_to_np_type(o::Type{UInt64}) = numpy.uint64
jl_to_np_type(o::Type{Int32}) = numpy.int32
jl_to_np_type(o::Type{Int64}) = numpy.int64
jl_to_np_type(o::Type{Float32}) = numpy.float32
jl_to_np_type(o::Type{Float64}) = numpy.float64
jl_to_np_type(o::Type{ComplexF32}) = numpy.complex64
jl_to_np_type(o::Type{ComplexF64}) = numpy.complex128


text:
chn = sample(mdl1, sim.data, HMC(0.05, 10), MCMCThreads(), 1000, 2)
post = posterior(mdl1, chn, sim.data)

@test expectation(post).latent.α ≈ 0.5 atol=0.1

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
include("bnf_utils.jl")
include("bnf.jl")


text:
struct DistortedView{T,P<:AbstractMatrix,E<:AbstractExtrapolation,G,D} <: AbstractArray{T,2}
    parent::P
    etp::E
    grid::G
    field::D

    function DistortedView(parent::AbstractMatrix{T}, grid::AbstractArray{Float64,3}) where T
        @assert size(grid,1) == 2
        etp = ImageTransformations.box_extrapolation(parent, Flat())
        field = ImageTransformations.box_extrapolation(grid, 0.0)
        new{T,typeof(parent),typeof(etp),typeof(grid),typeof(field)}(parent, etp, grid, field)
    end
end

Base.parent(A::DistortedView) = A.parent
Base.size(A::DistortedView) = map(length, indices(A.parent))

function ShowItLikeYouBuildIt.showarg(io::IO, A::DistortedView)
    print(io, typeof(A).name, '(')
    showarg(io, parent(A))
    print(io, ", ")
    showarg(io, A.grid)
    print(io, " as ", size(A.field,2), '×', size(A.field,3), " vector field")
    print(io, ')')
end

# showargs for SubArray{<:Colorant} is already implemented by ImageCore
function ShowItLikeYouBuildIt.showarg(io::IO, A::SubArray{<:Number,N,<:DistortedView}) where N
    print(io, "view(")
    showarg(io, parent(A))
    print(io, ", ")
    for (i, el) in enumerate(A.indexes)
        print(io, el)
        i < length(A.indexes) && print(io, ", ")
    end
    print(io, ')')
end

Base.summary(A::DistortedView) = summary_build(A)
Base.summary(A::SubArray{<:Number,N,<:DistortedView}) where {N} = summary_build(A)

# inline speeds up ~30%
@inline function Base.getindex(A::DistortedView, i::Int, j::Int)
    # unpack member variables
    parent = A.parent
    etp    = A.etp
    field  = A.field
    # size of the parent array
    indsy, indsx = indices(parent)
    leny,  lenx  = length(indsy), length(indsx)
    # grid size of the field field
    _, gh, gw = size(field)
    # map array indices to grid indices
    gi, gj = (i-1)/(leny-1)*(gh-1)+1, (j-1)/(lenx-1)*(gw-1)+1
    # compute parent indices offset
    @inbounds dy = field[1, gi, gj] * leny
    @inbounds dx = field[2, gi, gj] * lenx
    # compute parent indices and return value
    # Note: we subtract instead of add, because the vector field is
    #       specified in forward mode (i.e. in reference to the source
    #       pixel), while here we need to compute using backward mode
    #       (in reference to the target pixel)
    checkbounds(indsy, i)
    checkbounds(indsx, j)
    @inbounds y = indsy[i] - dy
    @inbounds x = indsx[j] - dx
    @inbounds res = etp[y, x]
    res
end


text:
@kwdef mutable struct SFixedGainDivHandler <: AbstractHandler
    id = 0
end

@register(SFixedGainDivHandler, ÷, begin
    [SBit, Number] => [SBit]
    [Vector{SBit}, Number] => [Vector{SBit}]
    [Matrix{SBit}, Number] => [Matrix{SBit}]
end)

function (handler::SFixedGainDivHandler)(netlist::Netlist,
                                         inputs::Vector{Variable},
                                         outputs::Vector{Variable})
    # update netlist with inputs
    setsigned!(netlist, getname(inputs[1]), true)

    # compute output size
    outsize = getsize(netlist, getname(outputs[1]))

    # add output net to netlist
    setsigned!(netlist, getname(outputs[1]), true)

    outstring = """
        $stdcomment
        // BEGIN fdiv$(handler.id)
        stoch_fixed_gain_div_mat #(
                .COUNTER_SIZE(8),
                .GAIN($(inputs[2].name)),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) fdiv$(handler.id)_p (
                .CLK(CLK),
                .nRST(nRST),
                .A($(inputs[1].name)_p),
                .Y($(outputs[1].name)_p)
            );
        stoch_fixed_gain_div_mat #(
                .COUNTER_SIZE(8),
                .GAIN($(inputs[2].name)),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) fdiv$(handler.id)_m (
                .CLK(CLK),
                .nRST(nRST),
                .A($(inputs[1].name)_m),
                .Y($(outputs[1].name)_m)
            );
        // END fdiv$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

allowconstreplace(::Type{SFixedGainDivHandler}) = false


text:
module ClosedIntervals

import Base.show, Base.isempty, Base.in, Base.length, Base.<<, Base.>>
import Base.isequal, Base.isless
import Base.*, Base.+, Base.==, Base.hash, Base.issubset
import Base.⊊, Base.⊇, Base.⊋

import AbstractLattices: ∨, ∧

export ClosedInterval, EmptyInterval, ⊇, ⊋
export show, left, right, .., ±, ∨, ∧

# Create the ClosedInterval type
struct ClosedInterval{T}
    L::T       # left end point
    R::T       # right end point
    nil::Bool  # signal if this is an empty interval
end

# Construct from two distinct end points
function ClosedInterval(l,r)
    (a,b) = promote(l,r)
    if a>b
        a,b = b,a
    end
    return ClosedInterval(a,b,false)
end

# dot-dot Construction
(..)(a,b) = ClosedInterval(a,b)

# plus-minus consruction
±(a,b) = ClosedInterval(a-b,a+b)

# Construct from a 2-tuple
function ClosedInterval(ab::Tuple{S,T}) where {S,T}
    return ClosedInterval(ab[1],ab[2]) # use 2-arg to test order
end

# Construct from one end point: assume L and R are the same
ClosedInterval(a) = ClosedInterval(a,a,false)

# Construction with no specified end points: assume [0,1]
function ClosedInterval(T::DataType = Float64)
    ClosedInterval(zero(T),one(T),false)
end

# Create an empty interval
"""
`EmptyInterval(T::DataType = Float64)` creates an empty
`ClosedInterval` of a given type.
"""
function EmptyInterval(T::DataType = Float64)
    return ClosedInterval(zero(T),zero(T), true)
end

# Fetch the left end point
"""
For a `ClosedInterval` `I`, `left(I)` returns its left end point.
"""
function left(J::ClosedInterval)
    if isempty(J)
        error("An empty interval does not have a left end point")
    end
    return J.L
end

# Fetch the right end point
"""
For a `ClosedInterval` `I`, `right(I)` returns its right end point.
"""
function right(J::ClosedInterval)
    if isempty(J)
        error("An empty interval does not have a right end point")
    end
    return J.R
end

# Is this an empty interval?
"""
For a `ClosedInterval I`, `isempty(I)` tests if `I` is an empty interval.
"""
isempty(J::ClosedInterval) = J.nil

# Print as a closed interval should be printed
function show(io::IO, J::ClosedInterval)
    if J.nil
        print(io,"[]")
    else
        print(io,"[", J.L, ",", J.R, "]")
    end
end

# The length of an interval is R-L unless it's empty. In that case, we
# throw an error.
"""
`length(I)` is the length of the `ClosedInterval` `I`.
"""
function length(J::ClosedInterval{T}) where T
    if isempty(J)
      return zero(T)
    end
    return J.R - J.L
end

# Test point for membership
"""
For a number `x` and a `ClosedInterval` `I` (of the same type)
`in(x,I)` tests if `x` is contained in the interval `I`.
"""
function in(x, J::ClosedInterval)
    return J.L <= x <= J.R
end

# The intersection of two intervals is the largest interval contained
# in both. If the intervals are disjoint, we return an empty interval.
"""
For `ClosedInterval`s `J` and `K`, `J*K` is their intersection.
Also available as `J ∧ K`.
"""
function *(J::ClosedInterval, K::ClosedInterval)

    a = max(J.L, K.L)   # left end point of result
    b = min(J.R, K.R)   # right end point of result

    a,b = promote(a,b)
    # if either interval is nil, so is their *
    if J.nil || K.nil
        return EmptyInterval(typeof(a))
    end

    if a>b  # uh oh, they're disjoint
        return EmptyInterval(typeof(a))
    end

    return ClosedInterval(a,b,false)
end

"""
`J ∧ K` is the largest `ClosedInterval` contained in both.
See also `*`.
"""
(∧)(J::ClosedInterval, K::ClosedInterval) = J*K

# The + of two intervals is the smallest interval containing them
# both. This is the same as their union if they overlap.

"""
For `ClosedInterval`s `J` and `K`, `J+K` is the smallest `ClosedInterval`
containing them both. This is also available as `J ∨ K`.
"""
function +(J::ClosedInterval, K::ClosedInterval)
    # The empty interval acts as an identity element for this
    # operation, so we check if either interval is empty first.

    if J.nil
        return K
    end
    if K.nil
        return J
    end

    a = min(J.L, K.L)   # left end point of result
    b = max(J.R, K.R)   # right end point of result

    return ClosedInterval(a,b,false)
end

"""
`J ∨ K` is the smalles `ClosedInterval` containing both. See also `+`.
"""
(∨)(J::ClosedInterval, K::ClosedInterval)= J+K

# Compare intervals for equality
function isequal(I::ClosedInterval, J::ClosedInterval)
    if I.nil && J.nil
        return true
    end
    if I.nil != J.nil
        return false
    end
    return I.L==J.L && I.R==J.R
end

==(I::ClosedInterval, J::ClosedInterval) = isequal(I,J)

# Sort intervals lexicographically, but put empty intervals at the
# bottom of the order.
"""
Lexicographic ordering of `ClosedInterval`s.
"""
function isless(I::ClosedInterval, J::ClosedInterval)
    if J.nil
        return false
    end
    if I.nil
       return true
    end

    if I.L < J.L
        return true
    end

    if I.L > J.L
        return false
    end

    return I.R < J.R
end

# We use << to mean "completely to the left of"
"""
For `ClosedInterval`s `I` and `J`, `I<<J` tests if `I` is completely
to the left of `J`.
"""
function <<(I::ClosedInterval, J::ClosedInterval)
    if I.nil || J.nil
        return false
    end
    return I.R < J.L
end

# Likewise, >> means "completely to the right of"
"""
For `ClosedInterval`s `I` and `J`, `I>>J` tests if `I` is completely
to the right of `J`.
"""
>>(I::ClosedInterval, J::ClosedInterval) = J << I

function hash(J::ClosedInterval, h::UInt64=UInt64(0))
  if isempty(J)
    return hash(false,h)
  end
  v = [J.L, J.R]
  return hash(v,h)
end

function issubset(J::ClosedInterval, K::ClosedInterval)
  if isempty(J)
    return true
  end
  if isempty(K)
    return false
  end

  return in(J.L,K) && in(J.R,K)
end

(⊇)(J::ClosedInterval,K::ClosedInterval) = issubset(K,J)

(⊊)(J::ClosedInterval,K::ClosedInterval) = issubset(J,K) && !(J==K)
(⊋)(J::ClosedInterval,K::ClosedInterval) = issubset(K,J) && !(J==K)

end # module


text:
getnames(check_function::Function) = filter(check_function, uint32(0:65534))

# gets all the names currently boundo to programs
getProgramNames()	  = getnames(glIsProgram)
getShaderNames() 	  = getnames(glIsShader)
getVertexArrayNames() = getnames(glIsVertexArray)

# display info for all active uniforms in a program
function getUniformsInfo(p::GLProgram)
	program = p.id
	# Get uniforms info (not in named blocks)
	@show activeUnif = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)

	for i=0:activeUnif-1
		@show index = glGetActiveUniformsiv(program, i, GL_UNIFORM_BLOCK_INDEX)
		if (index == -1)
			@show name 		     = glGetActiveUniformName(program, i)
			@show uniType 	   	 = glGetActiveUniformsiv(program, i, GL_UNIFORM_TYPE)

			@show uniSize 	   	 = glGetActiveUniformsiv(program, i, GL_UNIFORM_SIZE)
			@show uniArrayStride = glGetActiveUniformsiv(program, i, GL_UNIFORM_ARRAY_STRIDE)

			auxSize = 0
			if (uniArrayStride > 0)
				@show auxSize = uniArrayStride * uniSize
			else
				@show auxSize = spGLSLTypeSize[uniType]
			end
		end
	end
	# Get named blocks info
	@show count = glGetProgramiv(program, GL_ACTIVE_UNIFORM_BLOCKS)

	for i=0:count-1
		# Get blocks name
		@show name 	 		 = glGetActiveUniformBlockName(program, i)
		@show dataSize 		 = glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_DATA_SIZE)

		@show index 	 		 = glGetActiveUniformBlockiv(program, i,  GL_UNIFORM_BLOCK_BINDING)
		@show binding_point 	 = glGetIntegeri_v(GL_UNIFORM_BUFFER_BINDING, index)

		@show activeUnif   	 = glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_ACTIVE_UNIFORMS)

		indices = zeros(GLuint, activeUnif)
		glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_ACTIVE_UNIFORM_INDICES, indices)
		@show indices
		for ubindex in indices
			@show name 		   = glGetActiveUniformName(program, ubindex)
			@show uniType 	   = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_TYPE)
			@show uniOffset    = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_OFFSET)
			@show uniSize 	   = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_SIZE)
			@show uniMatStride = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_MATRIX_STRIDE)
		end
	end
end

# display the values for uniforms in the default block
function getUniformInfo(p::GLProgram, uniName::Symbol)
	# is it a program ?
	@show program 				  = p.id
	@show loc 	  				  = glGetUniformLocation(program, uniName)
	@show name, typ, uniform_size = glGetActiveUniform(program, loc)
end

# display the values for a uniform in a named block
function getUniformInBlockInfo(p::GLProgram,
				blockName,
				uniName)

	program = p.id

	@show index = glGetUniformBlockIndex(program, blockName)
	if (index == GL_INVALID_INDEX)
		println("$uniName is not a valid uniform name in block $blockName")
	end
	@show bindIndex 		= glGetActiveUniformBlockiv(program, index, GL_UNIFORM_BLOCK_BINDING)
	@show bufferIndex 		= glGetIntegeri_v(GL_UNIFORM_BUFFER_BINDING, bindIndex)
	@show uniIndex 			= glGetUniformIndices(program, uniName)

	@show uniType 			= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_TYPE)
	@show uniOffset 		= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_OFFSET)
	@show uniSize 			= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_SIZE)
	@show uniArrayStride 	= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_ARRAY_STRIDE)
	@show uniMatStride 		= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_MATRIX_STRIDE)
end

# display information for a program's attributes
function getAttributesInfo(p::GLProgram)

	program = p.id
	# how many attribs?
	@show activeAttr = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
	# get location and type for each attrib
	for i=0:activeAttr-1
		@show name, typ, siz = glGetActiveAttrib(program,	i)
		@show loc = glGetAttribLocation(program, name)
	end
end

# display program's information
function getProgramInfo(p::GLProgram)
	# check if name is really a program
	@show program = p.id
	# Get the shader's name
	@show shaders = glGetAttachedShaders(program)
	for shader in shaders
		@show info = GLENUM(convert(GLenum, glGetShaderiv(shader, GL_SHADER_TYPE))).name
	end
	# Get program info
	@show info = glGetProgramiv(program, GL_PROGRAM_SEPARABLE)
	@show info = glGetProgramiv(program, GL_PROGRAM_BINARY_RETRIEVABLE_HINT)
	@show info = glGetProgramiv(program, GL_LINK_STATUS)
	@show info = glGetProgramiv(program, GL_VALIDATE_STATUS)
	@show info = glGetProgramiv(program, GL_DELETE_STATUS)
	@show info = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
	@show info = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)
	@show info = glGetProgramiv(program, GL_ACTIVE_UNIFORM_BLOCKS)
	@show info = glGetProgramiv(program, GL_ACTIVE_ATOMIC_COUNTER_BUFFERS)
	@show info = glGetProgramiv(program, GL_TRANSFORM_FEEDBACK_BUFFER_MODE)
	@show info = glGetProgramiv(program, GL_TRANSFORM_FEEDBACK_VARYINGS)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.453305806680525 <= q <= 37.480198521751326)

                     add_NL_constraint(m, :(softplus(-0.5060529352779994 + 0.34605798298901425*softplus(-0.49627455608242954 + -0.46375151869975495*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.8504066842806357*softplus(-0.2321501407718678 + 0.5600091754803458*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.462211279142271*softplus(-0.9807289687619161 + -0.5968595287514287*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.0428411186394575*softplus(0.7899732054652109 + -0.5363355140855166*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + softplus(0.02992678336456356 + -0.8781954986928167*softplus(-0.49627455608242954 + -0.46375151869975495*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.84297522225665*softplus(-0.2321501407718678 + 0.5600091754803458*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.44937020244474857*softplus(-0.9807289687619161 + -0.5968595287514287*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.012040065223148932*softplus(0.7899732054652109 + -0.5363355140855166*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + softplus(-0.25761034404046734 + -0.19911102915318102*softplus(-0.49627455608242954 + -0.46375151869975495*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7457041629864394*softplus(-0.2321501407718678 + 0.5600091754803458*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.7638178262546895*softplus(-0.9807289687619161 + -0.5968595287514287*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.07889233901952375*softplus(0.7899732054652109 + -0.5363355140855166*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) + softplus(0.1153458350428731 + 0.9297971624965538*softplus(-0.49627455608242954 + -0.46375151869975495*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + 0.4918746175115887*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.38157591143072933*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.5439511060803288*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.7017891371439648*softplus(-0.2321501407718678 + 0.5600091754803458*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.0034071169835434567*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.6765465160367254*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + 0.8220752699067027*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + -0.8093911002605614*softplus(-0.9807289687619161 + -0.5968595287514287*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.3897253830589995*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + 0.5397807246416386*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.24558591797668194*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))) + 0.43763294749909276*softplus(0.7899732054652109 + -0.5363355140855166*softplus(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])) + -0.6727460735415476*softplus(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])) + -0.4129142282119118*softplus(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])) + -0.4150508411017726*softplus(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module FineShiftTests

using Test, FineShift, InterpolationKernels

function randu(T::Type{<:AbstractFloat}; max::Real=1, min::Real=0)
    _min = T(min)
    _max = T(max)
    return (_max - _min)*rand(T) + _min
end

function randu(T::Type{<:AbstractFloat}, dims::Dims; max::Real=1, min::Real=0)
    A = rand(T, dims)
    if max != 1 || min != 0
        a = T(max) - T(min)
        b = T(min)
        @inbounds @simd for i in eachindex(A)
            A[i] = a*A[i] + b
        end
    end
    return A
end

function FineShift.fineshift(dims::NTuple{N,Int},
                             src::AbstractArray{T,N},
                             ker::Kernel{T},
                             off::NTuple{N,Real},
                             adj::Bool = false) where {T<:AbstractFloat,N}
   return _fineshift(dims, src, ker, off, adj, Val(N))
end

function _fineshift(dims::NTuple{N,Int},
                    src::AbstractArray{T,N},
                    ker::Kernel{T},
                    off::NTuple{N,Real},
                    adj::Bool,
                    ::Val{1}) where {T<:AbstractFloat,N}
    return fineshift(dims[1], src, ker, off[1], 1, adj)
end

function _fineshift(dims::NTuple{N,Int},
                    src::AbstractArray{T,N},
                    ker::Kernel{T},
                    off::NTuple{N,Real},
                    adj::Bool,
                    ::Val{D}) where {T<:AbstractFloat,N,D}
    return fineshift(dims[D], _fineshift(dims, src, ker, off, adj, Val(D-1)),
                     ker, off[D], D, adj)
end

function vdot(x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tx<:Real,Ty<:Real,N}
    T = float(promote_type(Tx,Ty))
    local s::T = 0
    @assert axes(x) == axes(y)
    @inbounds for i in eachindex(x, y)
        s += T(x[i])*T(y[i])
    end
    return s
end

@testset "Fast correlation code" begin
    function splitindices(m::Int, n::Int, S::Int, l::Int)
        # This mimics the splitting of the 1D correlation code in 4 different
        # intervals.
        p = S - l
        i1 = 1 - p
        i2 = n - p
        i3 = l - 1 + n
        I1 = 1:min(i1,m)
        I2 = max(i1+1,1):min(i2,m)
        I3 = max(i2+1,1):min(i3-1,m)
        I4 = max(i3,1):m
        J = zeros(Int, S)
        if length(I2) > 0
            j = first(I2) - l
            for k in 1:S-1
                J[k+1] = clamp(j+k,1,n)
            end
        elseif length(I3) > 0
            j = first(I3) - l
            for k in 1:S-1
                J[k+1] = clamp(j+k,1,n)
            end
        end
        for i in I1
            checksplitindices(m, n, S, l, i, 1)
        end
        for i in I2
            J[1:S-1] = J[2:S]
            J[S] = p + i
            checksplitindices(m, n, S, l, i, J)
        end
        for i in I3
            J[1:S-1] = J[2:S]
            checksplitindices(m, n, S, l, i, J)
        end
        for i in I4
            checksplitindices(m, n, S, l, i, n)
        end
        return true
    end

    function checksplitindices(m::Int, n::Int, S::Int, l::Int,
                               i::Int, Jp::Vector{Int})
        for k = 1:S
            j = clamp(i-l+k,1,n)
            jp = Jp[k]
            jp == j || error("(i,k)=($i,$k) j=$j, not $jp")
        end
    end

    function checksplitindices(m::Int, n::Int, S::Int, l::Int,
                               i::Int, jp::Int)
        for k = 1:S
            j = clamp(i-l+k,1,n)
            jp == j || error("(i,k)=($i,$k) j=$j, not $jp")
        end
    end

    # For offset t and kernel support size S, l = floor(S+1+t), hence to
    # explore all cases we want t in -(S+1):(S+1), that is l in 0:2(S+1)
    # with S = 4 (as assumed in the tests) we consider l in -1:10
    @testset "(l,m,n) = ($l,$m,$n)" for m in (3,6), n in (2,4,8), l in -1:10
        @test splitindices(m, n, 4, l)
    end
end

@testset "Fine shift with zero offset" begin
    @testset "Dimensions: $dims" for dims in ((100,), (20,30), (10,20,30))
        T = Float64
        N = length(dims)
        ker = CatmullRomSpline{T}()
        src = randu(T, dims; min=-1, max=1)
        off = ntuple(i -> zero(T), N)
        dst = fineshift(dims, src, ker, off)
        @test dst == src
    end
end

@testset "Direct/adjoint fine shift" begin
    @testset "Dimensions: $dims" for dims in ((100,), (20,30), (10,20,30))
        T = Float64
        N = length(dims)
        ker = CatmullRomSpline{T}()
        xdims = dims
        ydims = ntuple(d -> (isodd(d) ? dims[d] - 1 : dims[d] + 2), N)
        x = randu(T, xdims; min=-1, max=1)
        y = randu(T, ydims; min=-1, max=1)
        off = ntuple(i -> randu(T; min=-3, max=3), N)
        @test (vdot(y, fineshift(ydims, x, ker, off, false)) ≈
               vdot(x, fineshift(xdims, y, ker, off, true )))
    end
end

@testset "Fine shift of a smooth function" begin
    function f(x::T) where {T<:AbstractFloat}
        q = T(0.03)
        r = T(1.2)
        return exp(-q*x*x)*cos(x - r)
    end
    s = 5.3 # offset in number of samples
    t = -40:0.05:40
    y1 = f.(t)
    z1 = fineshift(length(t), y1, CatmullRomSpline(), s)
    z2 = f.(t .- s*step(t))
    @test maximum(abs.(z1 .- z2)) < 1e-5
end

end # module


text:
"""
    This implementation is based on the port from GSL to Python by
    J. Michael Burgess for his pynchrotron package:
    https://github.com/grburgess/pynchrotron
"""

module SynchrotronKernel

    include("F.jl")
    include("G.jl")

    export synchrotron_kernel,
           synchrotron_intensity,
           synchrotron_polarisation,
           ℱ, 𝒢

    """
        synchrotron_kernel(x::Real)

    Computes the first synchrotron function and the polarisation components at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.
    Returns a tuple `(K_tot, K_ort, K_par)`.

    ```julia
    K_tot = F(x)
    K_ort = 0.5 * (F(x) + G(x))
    K_par = 0.5 * (F(x) - G(x))
    ```

    ## Synchrotron Functions
    ``F(x) = x \\int_x^\\infty K_{\\frac{5}{3}}(t) dt``

    ``G(x) = x K_{\\frac{2}{3}}(x)``

    """
    function synchrotron_kernel(x::Real)
        # first synchrotron function
        K  = F(x)
        # second synchrotron function
        pK = G(x)

        # orthogonal and parallel polarisation components
        K_ort = 0.5 * (K + pK)
        K_par = 0.5 * (K - pK)

        return K, K_ort, K_par
    end

    """
        synchrotron_intensity(x::Real)

    Computes the total synchrotron kernel, without polarisation components.
    Wrapper for [`F`](@ref).

    ``F(x) = x \\int_x^\\infty K_{\\frac{5}{3}}(t) dt``
    """
    synchrotron_intensity(x::Real) = F(x)

    """
        synchrotron_polarisation(x::Real)

    Computes the second synchrotron function at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.
    Returns a tuple of `(K_ort, K_par)`.

    ```julia
    K_ort = 0.5 * (F(x) + G(x))
    K_par = 0.5 * (F(x) - G(x))
    ```

    ## Synchrotron Functions
    ``F(x) = x \\int_x^\\infty K_{\\frac{5}{3}}(t) dt``

    ``G(x) = x K_{\\frac{2}{3}}(x)``
    """
    function synchrotron_polarisation(x::Real)

        # first synchrotron function
        K  = F(x)
        # second synchrotron function
        pK = G(x)

        K_ort = 0.5 * (K + pK)
        K_par = 0.5 * (K - pK)

        # orthogonal and parallel polarisation components
        return K_ort, K_par
    end

end # module


text:
n = 512;

f = rand(n,n);

h = filt_gen(WT.db2)

g = [0, h[length(h):-1:2]...] .* ( (-1).^(1:length(h)));

Jmax = round(Int64,log2(n))-1;
Jmin = 1;
fW = copy(f);

for j=Jmax:-1:Jmin
    A = fW[1:2^(j+1),1:2^(j+1)];
    for d=1:2
        Coarse = subsampling(cconvol(A,h,d),d);
        Detail = subsampling(cconvol(A,g,d),d);
        A = cat(Coarse, Detail,dims=d);
    end
    fW[1:2^(j+1),1:2^(j+1)] = A;
    j1 = Jmax-j;
end

@test size(fW) == size(f)


text:
module JSONPointer

import OrderedCollections

include("pointer.jl")

export @j_str

end # module


text:
using Documenter, Dispersal, Weave, IJulia

basedir = @__DIR__

example = joinpath(basedir, "src/example.jmd")

mdpath = joinpath(basedir, "src/example.md")
notebookdir = joinpath(basedir, "build/notebook")

mkpath(joinpath(basedir, "build/assets"))
mkpath(notebookdir)

# Generate examples latex and images
weave(example, out_path=mdpath, doctype="github")

# Generate examples notebook
convert_doc(example, joinpath(notebookdir, "example.ipynb"))

# Generate HTML docs
makedocs(
    modules = [Dispersal],
    sitename = "Dispersal.jl",
    pages = [
        "Home" => "index.md",
        "Examples" => "example.md"
    ],
    clean = false,
    checkdocs = :all,
    strict = true,
)

deploydocs(
    repo = "github.com/cesaraustralia/Dispersal.jl.git",
)


text:
Actions = Dict('N' => [0, 1], 'S' => [0, -1], 'E' => [1, 0], 'W' => [-1, 0])

struct Instr
    act::Char
    val::Int
end

puzzleinput = readlines("inputs/day12.txt")
puzzleinput = map(s -> Instr(s[1], parse(Int, s[2:end])), puzzleinput)

# First Half

struct Ship
    pos::Vector{Int} # Ship position on the plane (x, y)
    dir::Vector{Int} # Ship direction as unit vector (eg. ovest = (-1,0))
end
# The ship status

ship = Ship([0, 0], [1, 0]) # The ship starts by facing east.

function handle1(s::Ship, i::Instr)
    # From a ship status s, handles an instruction i returning the new ship status

    act, val = i.act, i.val
    if act == 'L' || act == 'R'
        val = (act == 'R' ? -1 : 1) * val
        #  Right means clockwise rotation = negatives angles in trigs functions

        θ = deg2rad(val)
        R = trunc.(Int, [cos(θ) -sin(θ); sin(θ) cos(θ)])
        return Ship(s.pos, R * s.dir)
        # To perform rotation, multiply by rotation matrix

    end
    if act == 'F'
        v = val * s.dir
    else
        v = val * Actions[act]
    end
    Ship(s.pos + v, s.dir)
end

wherelead = reduce(handle1, puzzleinput, init=ship)

sum(map(abs, wherelead.pos))

# Second Half

function handle2(sw::Tuple{Ship,Array{Int,1}}, i::Instr)
    # From a ship status s and a waypoint w,
    # handles an instruction i returning the new ship status and waypoint

    s, w = sw[1], sw[2]
    act, val = i.act, i.val
    if act == 'L' || act == 'R'
        val = (act == 'R' ? -1 : 1) * val
        θ = deg2rad(val)
        R = trunc.(Int, [cos(θ) -sin(θ); sin(θ) cos(θ)])
        return (s, R * w)
    end
    if act == 'F'
        v = val * w
        s = Ship(s.pos + v, s.dir)
    else
        v = val * Actions[act]
        w += v
    end
    (s, w)
end

waypoint = ship.pos + [10, 1]
# The waypoint starts 10 units east and 1 unit north relative to the ship.

wherelead = reduce(handle2, puzzleinput, init=(ship, waypoint))

sum(map(abs, wherelead[1].pos))


text:
import Knet.Ops21: gelu, geluback
import Base.Broadcast: broadcasted
import Knet
using Knet.KnetArrays: KnetArray, DevArray, Bcasted
using CUDA: CuArray, CuPtr
using Knet.LibKnet8: @knet8

for (R,P) in ((KnetArray,Ptr), (CuArray,CuPtr)), T in (Float32,Float64); S = sizeof(T) * 8
    for f in ("gelu",)
        J, F = Symbol(f), "$(f)_$S"; M = which(@__MODULE__,J)
        @eval begin
            function broadcasted(::typeof($J),x::$R{$T})
                y = similar(x)
                @knet8($F,(Cint,$P{$T},$P{$T}),length(y),x,y)
                return y
            end
            # Bcasted methods -- only needed for KnetArray
            ($M).$J(x::Bcasted{<:$R{$T}}) = broadcasted($J, x.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}) = broadcasted($J, x.value) |> Bcasted
        end
    end
    for f in ("geluback",)
        J, F = Symbol(f), "$(f)_$(S)_11"; M = which(@__MODULE__,J)
        @eval begin
            function broadcasted(::typeof($J),x::$R{$T},y::$R{$T})
                z = similar(x)
                @knet8($F,(Cint,$P{$T},$P{$T},$P{$T}),length(z),x,y,z)
                return z
            end
            # Bcasted methods -- only needed for KnetArray
            ($M).$J(x::Bcasted{<:$R{$T}}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            ($M).$J(x::$R{$T}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            ($M).$J(x::Bcasted{<:$R{$T}}, y::$R{$T}) = broadcasted($J, x.value, y) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            broadcasted(::typeof($J),x::$R{$T}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}, y::$R{$T}) = broadcasted($J, x.value, y) |> Bcasted
        end
    end
end


text:
using LinearAlgebra
using Multigraphs, Graphs
using OMEinsum

@enum SCGate begin
    Single
    FSim
    Id
    Noise
end
struct RQCLayout{VT, PT, GT}
    g::Multigraph{VT}
    edge_patterns::Dict{Edge{VT}, PT}
    gates::Dict{VT, Vector{GT}}
    pattern_loop::Vector{PT}
    qubit_locs::Dict{VT, Tuple{Float64, Float64}}
end

function RQCLayout{VT, PT, GT}(vs::Vector{VT}, edge_pattern::Dict{Edge{VT}, PT}, pattern_loop::Vector{PT},
        qubit_locs::Dict{VT, Tuple{Float64, Float64}} = Dict{VT, Tuple{Float64, Float64}}()) where {VT, PT, GT}
    g = Multigraph{VT}(Dict(v => VT[] for v in vs), maximum(vs))
    for (e, et) in edge_pattern
        add_edge!(g, e)
    end
    gates = Dict(v => GT[] for v in vertices(g))
    return RQCLayout{VT, PT, GT}(g, edge_pattern, gates, pattern_loop, qubit_locs)
end
RQCLayout{VT, PT, GT}(nbits::VT, edge_pattern::Dict{Edge{VT}, PT}, pattern_loop::Vector{PT},
        qubit_locs::Dict{VT, Tuple{Float64, Float64}} = Dict{VT, Tuple{Float64, Float64}}()) where {VT, PT, GT} =
    RQCLayout{VT, PT, GT}(VT[i for i = 1:nbits], edge_pattern, pattern_loop, qubit_locs)

Base.copy(layout::RQCLayout) = deepcopy(layout)
Base.deepcopy(layout::RQCLayout{VT, PT, GT}) where {VT, PT, GT} =
    RQCLayout{VT, PT, GT}(deepcopy(layout.g), deepcopy(layout.edge_patterns),
        deepcopy(layout.gates), deepcopy(layout.pattern_loop), deepcopy(layout.qubit_locs))
Graphs.vertices(layout::RQCLayout) = vertices(layout.g)
Graphs.has_vertex(layout::RQCLayout, v) = has_vertex(layout.g, v)
Graphs.neighbors(layout::RQCLayout, v) = neighbors(layout.g, v)
function Graphs.connected_components(g::RQCLayout{VT, PT, SCGate}) where {VT, PT}
    vs_remain = Set(collect(vertices(g)))
    comps = []
    while !isempty(vs_remain)
        v0 = pop!(vs_remain)
        comp_v0 = [v0]
        comp_remain = [v0]
        while !isempty(comp_remain)
            v = pop!(comp_remain)
            gates_v = gates(g, v)
            for i = 1:length(gates_v)
                gates_v[i] === FSim || continue
                u = partner(g, v, depth_to_cycle(i))
                !(u in comp_v0) && push!(comp_remain, u)
                !(u in comp_v0) && push!(comp_v0, u)
                delete!(vs_remain, u)
            end
        end
        push!(comps, sort!(comp_v0))
    end
    return sort!(comps; by = x -> first(x))
end

gates(layout::RQCLayout, v) = layout.gates[v]
cycle_pattern(layout::RQCLayout, ncycles) = getindex(layout.pattern_loop, rem(ncycles-1, length(layout.pattern_loop))+1)
edge_pattern(layout::RQCLayout, u, v) = layout.edge_patterns[Edge(min(u, v), max(u, v))]
function partner(layout::RQCLayout, v, ncycles)
    pt = cycle_pattern(layout, ncycles)
    nbs = neighbors(layout, v)
    idx = findfirst(u -> edge_pattern(layout, u, v) === pt, nbs)
    idx !== nothing && return nbs[idx]
    return nothing
end
has_partner(g, v, ncycles) = (partner(g, v, ncycles) !== nothing)
depth_to_cycle(d) = d÷2

function depth(g::RQCLayout)
    v0 = first(keys(g.gates))
    return length(gates(g, v0))
end

function flux(g::RQCLayout{VT, PT, SCGate}, vs = collect(vertices(g))) where {VT, PT}
    M = [gates(g, v)[i] !== Noise for v in vs, i in 1:depth(g)]
    return minimum(sum(M, dims = 1))
end