text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function add_elec_load_balance_constraints(m, p; _n="")

	##Constraint (8a): Electrical Load Balancing with Grid
    if isempty(p.s.electric_tariff.export_bins)
        conrefs = @constraint(m, [ts in p.time_steps_with_grid],
            sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.techs.elec) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.elec ) +
            sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) ==
            sum( sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec)
                + m[Symbol("dvCurtail"*_n)][t, ts] for t in p.techs.elec)
            + sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec)
            + sum(m[Symbol("dvThermalProduction"*_n)][t, ts] / p.cop[t] for t in p.techs.cooling)
            + p.s.electric_load.loads_kw[ts]
            - p.s.cooling_load.loads_kw_thermal[ts] / p.cop["ExistingChiller"]
        )
    else
        conrefs = @constraint(m, [ts in p.time_steps_with_grid],
            sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.techs.elec) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.elec ) +
            sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) ==
            sum(  sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec)
                + sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t])
                + m[Symbol("dvCurtail"*_n)][t, ts]
            for t in p.techs.elec)
            + sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec)
            + sum(m[Symbol("dvThermalProduction"*_n)][t, ts] / p.cop[t] for t in p.techs.cooling)
            + p.s.electric_load.loads_kw[ts]
            - p.s.cooling_load.loads_kw_thermal[ts] / p.cop["ExistingChiller"]
        )
    end

	for (i, cr) in enumerate(conrefs)
		JuMP.set_name(cr, "con_load_balance"*_n*string("_t", i))
	end

	##Constraint (8b): Electrical Load Balancing without Grid
	if !p.s.settings.off_grid_flag # load balancing constraint for grid-connected runs
        @constraint(m, [ts in p.time_steps_without_grid],
            sum(p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.techs.elec) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.elec )  ==
            sum( sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) +
            m[Symbol("dvCurtail"*_n)][t, ts] for t in p.techs.elec) +
            (p.s.electric_load.critical_loads_kw[ts] )
        )
    else # load balancing constraint for off-grid runs
        @constraint(m, [ts in p.time_steps_without_grid],
            sum(p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.techs.elec) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.elec )  ==
            sum( sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) +
            m[Symbol("dvCurtail"*_n)][t, ts] for t in p.techs.elec) +
            (p.s.electric_load.critical_loads_kw[ts] * m[Symbol("dvOffgridLoadServedFraction"*_n)][ts] )
        )
        ##Constraint : For off-grid scenarios, annual load served must be >= minimum percent specified
        @constraint(m, sum(m[Symbol("dvOffgridLoadServedFraction"*_n)][ts] * p.s.electric_load.critical_loads_kw[ts] for ts in p.time_steps_without_grid) >=
			sum(p.s.electric_load.critical_loads_kw) * p.s.electric_load.min_load_met_annual_pct
		)
    end

end

function add_production_constraints(m, p; _n="")
	# Constraint (4d): Electrical production sent to storage or export must be less than technology's rated production
    if isempty(p.s.electric_tariff.export_bins)
        @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec)
          + m[Symbol("dvCurtail"*_n)][t, ts]
         <= p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
        )
    else
        @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec)
          + m[Symbol("dvCurtail"*_n)][t, ts]
          + sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t])
         <= p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
        )
    end

	# Constraint (4e): Electrical production sent to storage or grid must be less than technology's rated production - no grid
	@constraint(m, [t in p.techs.elec, ts in p.time_steps_without_grid],
		sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec)  <=
		p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
	)

end

function add_thermal_load_constraints(m, p; _n="")

    m[Symbol("binFlexHVAC"*_n)] = 0
    m[Symbol("dvTemperature"*_n)] = 0
    m[Symbol("dvComfortLimitViolationCost"*_n)] = 0

    if !isnothing(p.s.flexible_hvac)
        #= FlexibleHVAC does not require equality constraints for thermal loads. The thermal loads
        are instead a function of the energy required to keep the space temperature within the
        comfort limits.
        =#
        add_flexible_hvac_constraints(m, p, _n=_n)

	##Constraint (5b): Hot thermal loads
    else
        if !isempty(p.techs.heating)

            # if !isempty(p.SteamTurbineTechs)
            #     @constraint(m, HotThermalLoadCon[ts in p.time_steps],
            #             # sum(m[Symbol("dvThermalProduction"*_n)][t,ts] - m[:dvThermalToSteamTurbine][t,ts] for t in p.CHPTechs) +
            #             # sum(m[Symbol("dvThermalProduction"*_n)][t,ts] for t in p.SteamTurbineTechs) +
            #             sum(p.production_factor[t,ts] * (m[Symbol("dvThermalProduction"*_n)][t,ts] - m[:dvThermalToSteamTurbine][t,ts]) for t in p.techs.boiler)
            #             # + sum(p.GHPHeatingThermalServed[g,ts] * m[:binGHP][g] for g in p.GHPOptions)
            #             # + sum(m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.hot)
            #             ==
            #             p.HeatingLoad[ts] * p.s.existing_boiler.efficiency
            #             # + sum(m[:dvProductionToWaste][t,ts] for t in p.CHPTechs) +
            #             # sum(m[:dvProductionToStorage][b,t,ts] for b in p.s.storage.types.hot, t in p.techs.heating) +
            #             # sum(m[Symbol("dvThermalProduction"*_n)][t,ts] / p.thermal_cop[t] for t in p.techs.absorption_chiller)
            #     )
            # else
                @constraint(m, [ts in p.time_steps],
                        sum(m[Symbol("dvThermalProduction"*_n)][t,ts] for t in union(p.techs.boiler, p.techs.chp)) +
                        # TODO do all thermal techs have production_factor = 1 ? get rid of it if so
                        + sum(m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.hot)
                        # + sum(p.GHPHeatingThermalServed[g,ts] * m[:binGHP][g] for g in p.GHPOptions)
                        ==
                        (p.s.dhw_load.loads_kw[ts] + p.s.space_heating_load.loads_kw[ts])
                        + sum(m[Symbol("dvProductionToWaste"*_n)][t,ts] for t in p.techs.chp) +
                        sum(m[:dvProductionToStorage][b,t,ts] for b in p.s.storage.types.hot, t in union(p.techs.boiler, p.techs.chp))  +
                        sum(m[Symbol("dvThermalProduction"*_n)][t,ts] / p.thermal_cop[t] for t in p.techs.absorption_chiller)
                )
            # end

        end

        if !isempty(p.techs.cooling)

            ##Constraint (5a): Cold thermal loads
            @constraint(m, [ts in p.time_steps_with_grid],
                    sum(m[Symbol("dvThermalProduction"*_n)][t,ts] for t in p.techs.cooling) +
                    sum(m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.s.storage.types.cold) ==
                    p.s.cooling_load.loads_kw_thermal[ts] +
                    sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for b in p.s.storage.types.cold, t in p.techs.cooling) #-
                    # sum(p.GHPCoolingThermalServed[g,ts] * m[:binGHP][g] for g in p.GHPOptions) +
            )

            # @constraint(m, [ts in p.time_steps_with_grid],
            #     sum(m[Symbol("dvThermalProduction"*_n)][t, ts] for t in p.techs.cooling) ==
            #     p.s.cooling_load.loads_kw_thermal[ts]
            # )
        end
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
using Hive
using Base.Test
using DataFrames

function open_database(f::Function)
    session = HiveSession()
    try
        result = execute(session, "use default")
        @test round(Int, result) == 0
        f(session)
    finally
        close(session)
    end
end

function fetch_server_metadata(session)
    server_name = get_info(session, InfoType.CLI_SERVER_NAME)
    @test length(server_name) > 0
    println("Server name: $server_name")

    cats = catalogs(session)
    @test isa(cats, DataFrame)
    println("Catalogs:")
    println(cats)

    sch = schemas(session)
    # at least the default schema should be present
    @test prod(size(sch)) > 1
    println("Schemas:")
    println(sch)

    tbls = tables(session)
    @test isa(tbls, DataFrame)
    println("Tables:")
    println(tbls)

    ttypes = tabletypes(session)
    @test isa(ttypes, DataFrame)
    println("Table types:")
    println(ttypes)

    cols = columns(session)
    @test isa(cols, DataFrame)
    println("Columns:")
    println(cols)

    fns = functions(session, "%")
    # at least the default functions should be present
    @test prod(size(fns)) > 1
    println("Functions:")
    println(fns)
    nothing
end

function create_table(session)
    rs = execute(session, "show tables like 'twitter_small'")
    table_exists = prod(size(dataframe(rs))) > 0

    if table_exists
        println("Use existing table: twitter_small")
    else
        println("Create table: twitter_small")
        result = execute(session, "create table twitter_small (fromid int, toid int) row format delimited fields terminated by ',' lines terminated by '\n' stored as textfile")
        @test round(Int, result) == 0.0
    end

    rs = execute(session, "select count(*) from twitter_small")
    rowcount = dataframe(rs)[1,1]
    if rowcount > 0
        println("Use existing data: $rowcount rows")
    else
        rowcount = 10^4
        println("Insert data: $rowcount rows")
        filename = "/tmp/twitter_small"
        open(filename, "w") do io
            writedlm(io, convert(Array{Int}, rand(UInt16, rowcount, 2)), ',')
        end
        result = execute(session, "load data local inpath '$filename' into table twitter_small")
        @test round(Int, result) == 0.0
    end
end

function fetch_records(session)
    rs = execute(session, "select min(fromid), max(fromid) from twitter_small")
    minmax = dataframe(rs)
    maxval = ceil(Int, mean([minmax[1,1], minmax[1,2]]))
    lim = 10000

    println("Execute, record iterator:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim")
    cnt = 0
    for rec in records(rs)
       println(rec)
        cnt += 1
    end
    close(rs)
    @test cnt <= lim

    println("Execute, dataframe iterator:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim")
    cnt = 0
    for frames in dataframes(rs)
       println(frames)
        cnt += size(frames, 1)
    end
    close(rs)
    @test cnt <= lim

    println("Execute, async:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim"; async=true)
    while !isready(rs)
        println("waiting...")
        sleep(10)
    end
    rs = result(rs)
    df = dataframe(rs)
    @test size(df, 1) <= lim
    println(df)
    nothing
end

open_database() do session
    create_table(session)
    fetch_server_metadata(session)
    fetch_records(session)
end

ext: jl


text:
module PETSc

using MPI, LinearAlgebra, SparseArrays

using PETSc_jll

include("const.jl")
include("lib.jl")
include("init.jl")
include("ref.jl")
include("viewer.jl")
include("options.jl")
include("vec.jl")
include("mat.jl")
include("matshell.jl")
include("ksp.jl")
include("pc.jl")
include("snes.jl")

end


text:
#=
Author : Ng Yeesian
Modified : Guy Farmer / Zander Tedjo / Michael Hilborn
defines various deployment and dispatch models and functions that will be run on them
=#
# using Gurobi
abstract type DeploymentModel end

abstract type DispatchModel end

#=function respond_to!(redeploy::RedeployModel, i::Int, t::Int)
    @assert length(redeploy.ambulances[i]) > 0 "$i: $(redeploy.ambulances[i])"
    amb = popfirst!(redeploy.ambulances[i])
    @assert amb != 0
    @assert redeploy.status[amb] == :available "$amb: $(redeploy.status[amb])"
    redeploy.status[amb] = :responding
    redeploy.fromtime[amb] = t
    amb
end

function arriveatscene!(redeploy::RedeployModel, amb::Int, t::Int)
    @assert redeploy.status[amb] == :responding "$amb: $(redeploy.status[amb])"
    @assert redeploy.hospital[amb] == 0
    redeploy.status[amb] = :atscene
    redeploy.fromtime[amb] = t
end

function conveying!(redeploy::RedeployModel, amb::Int, hosp::Int, t::Int)
    @assert redeploy.status[amb] == :atscene "$amb: $(redeploy.status[amb])"
    @assert redeploy.hospital[amb] != 0
    redeploy.status[amb] = :conveying
    redeploy.fromtime[amb] = t
    redeploy.hospital[amb] = hosp
end

function returning_to!(redeploy::RedeployModel, amb::Int, t::Int)
    @assert redeploy.status[amb] == :conveying "$amb: $(redeploy.status[amb])"
    @assert redeploy.hospital[amb] != 0
    redeploy.status[amb] = :returning
    redeploy.fromtime[amb] = t
    redeploy.assignment[amb]
end

function redeploying_to!(redeploy::RedeployModel, amb::Int, i::Int, t::Int)
    # DEBUG:
    println("redeploying amb $amb from $(redeploy.assignment[amb]) to $i")
    ambulances = redeploy.ambulances[redeploy.assignment[amb]]
    @assert !in(amb, redeploy.ambulances[i])
    @assert in(amb, ambulances)
    deleteat!(ambulances, findfirst(ambulances, amb))
    @assert !in(amb, ambulances)
    redeploy.assignment[amb] = i
    redeploy.status[amb] = :redeploying
    redeploy.fromtime[amb] = t
end

function returned_to!(redeploy::RedeployModel, amb::Int, t::Int)
    @assert in(redeploy.status[amb], (:returning, :redeploying)) "$amb: $(redeploy.status[amb])"
    redeploy.hospital[amb] = 0
    redeploy.status[amb] = :available
    redeploy.fromtime[amb] = t
    @assert !in(amb, redeploy.ambulances[redeploy.assignment[amb]])
    push!(redeploy.ambulances[redeploy.assignment[amb]], amb)
end

function redirected!(redeploy::RedeployModel, amb::Int, t::Int)
    @assert in(redeploy.status[amb], (:returning, :redeploying)) "$amb: $(redeploy.status[amb])"
    redeploy.hospital[amb] = 0
    redeploy.status[amb] = :responding
    redeploy.fromtime[amb] = t
end
=#

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: michaelhilborn/AmbulanceDeployement


text:
using SumOfSquares

const polyhedra_test = joinpath(Pkg.dir("SumOfSquares"), "test")

sostest(s) = include(joinpath(polyhedra_test, "$s.jl"))

sostest("solvers")

sostest("certificate")

sostest("motzkin")

# SOSTools demos
sostest("sospoly")
sostest("sosdemo2")
sostest("sosdemo3")
sostest("sosdemo4")
sostest("sosdemo5")
sostest("sosdemo6")
sostest("domain")
sostest("sosmatrix")

function sosprogram(vartable, decvartable)

end

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
using Documenter, BitConverter

makedocs(sitename="BitConverter",
         doctest=true)


text:
module Tableaus

    import GenericLinearAlgebra

    using CompactBasisFunctions
    using QuadratureRules
    using Reexport

    @reexport using RungeKutta

    using GeometricBase
    using GeometricBase.Config
    using GeometricBase.Utils

    using ..Integrators
    using ..Integrators.SPARK
    using ..Integrators.VPRK

    using GeometricIntegrators.Integrators.SPARK: get_ã_vspark_primary,
                                                  get_α_vspark_primary,
                                                  compute_ã_vspark_primary,
                                                  compute_α_vspark_primary

    export TableauGauss

    include("tableaus/coefficients_glrk.jl")
    include("tableaus/coefficients_lob.jl")

    include("tableaus/tableaus_spark.jl")

    export TableauSPARKGLRK,
           TableauSPARKLobattoIIIAIIIB,
           TableauSPARKLobattoIIIBIIIA,
           TableauSPARKGLRKLobattoIIIAIIIB,
           TableauSPARKGLRKLobattoIIIBIIIA,
           TableauSPARKLobatto,
           TableauSPARKLobABC,
           TableauSPARKLobABD,
           TableauSPARKVPRK,
           TableauSPARKGLVPRK

    include("tableaus/tableaus_vprk.jl")

    export TableauVPGLRK,
           TableauVPLobattoIIIA,
           TableauVPLobattoIIIB,
           TableauVPLobattoIIIC,
           TableauVPLobattoIIIC̄,
           TableauVPLobattoIIID,
           TableauVPLobattoIIIE,
           TableauVPLobattoIIIF,
           TableauVPLobattoIIIG,
           TableauVPSRK3,
           TableauVPLobattoIIIAIIIA,
           TableauVPLobattoIIIBIIIB,
           TableauVPRadauIIAIIA

    export TableauSymplecticProjection,
           TableauLobattoIIIAIIIBpSymplectic,
           TableauLobattoIIIBIIIApSymplectic,
           TableauGausspSymplectic

    include("tableaus/tableaus_vpark.jl")

    export TableauVSPARKLobattoIIIAIIIBProjection,
           TableauVSPARKLobattoIIIBIIIAProjection,
           TableauVSPARKModifiedLobattoIIIAIIIBProjection,
           TableauVSPARKModifiedLobattoIIIBIIIAProjection,
           TableauVSPARKInternalProjection,
           TableauVSPARKModifiedInternalProjection,
           TableauVSPARKMidpointProjection,
           TableauVSPARKModifiedMidpointProjection,
           TableauVSPARKSymmetricProjection,
           TableauVSPARKSymplecticProjection,
           TableauVSPARKGLRKpLobattoIIIAIIIB,
           TableauVSPARKGLRKpLobattoIIIBIIIA,
           TableauVSPARKGLRKpModifiedLobattoIIIAIIIB,
           TableauVSPARKGLRKpModifiedLobattoIIIBIIIA,
           TableauVSPARKGLRKpInternal,
           TableauVSPARKGLRKpModifiedInternal,
           TableauVSPARKGLRKpMidpoint,
           TableauVSPARKGLRKpModifiedMidpoint,
           TableauVSPARKGLRKpSymmetric,
           TableauVSPARKGLRKpSymplectic,
           TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB,
           TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB,
           TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA,
           TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA,
           TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB,
           TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA,
           TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB,
           TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA,
           TableauVSPARKLobattoIIIAIIIBpMidpoint,
           TableauVSPARKLobattoIIIBIIIApMidpoint,
           TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint,
           TableauVSPARKLobattoIIIBIIIApModifiedMidpoint,
           TableauVSPARKLobattoIIIAIIIBpSymmetric,
           TableauVSPARKLobattoIIIBIIIApSymmetric,
           TableauVSPARKLobABCCD,
           TableauVSPARKLobABCCE,
           TableauVSPARKLobABDE,
           TableauVSPARKLobABED,
           TableauVSPARKLobABD,
           TableauVSPARKLobABE,
           TableauVSPARKLobDE,
           TableauVSPARKLobED

    include("tableaus/tableaus_vspark_primary.jl")

    export TableauVSPARKLobattoIIIAB,
           TableauVSPARKLobattoIIIBA,
           TableauVSPARKLobattoIIICC̄,
           TableauVSPARKLobattoIIIC̄C,
           TableauVSPARKLobattoIIID,
           TableauVSPARKLobattoIIIE,
           TableauVSPARKGLRKLobattoIIIAB,
           TableauVSPARKGLRKLobattoIIIBA,
           TableauVSPARKGLRKLobattoIIICC̄,
           TableauVSPARKGLRKLobattoIIIC̄C,
           TableauVSPARKGLRKLobattoIIID,
           TableauVSPARKGLRKLobattoIIIE

    include("tableaus/tableaus_vspark_secondary.jl")

    export getTableauHPARK,
           TableauHPARKGLRK,
           TableauHPARKLobattoIIIAIIIB,
           TableauHPARKLobattoIIIBIIIA

    include("tableaus/tableaus_hpark.jl")

    export TableauHSPARKSymmetricProjection,
           TableauHSPARKGLRKpSymmetric,
           TableauHSPARKLobattoIIIAIIIBpSymmetric,
           TableauHSPARKLobattoIIIBIIIApSymmetric

    include("tableaus/tableaus_hspark_primary.jl")

    export TableauHSPARKLobattoIIIAB,
           TableauHSPARKLobattoIIIBA,
           TableauHSPARKLobattoIIICC̄,
           TableauHSPARKLobattoIIIC̄C,
           TableauHSPARKLobattoIIID,
           TableauHSPARKLobattoIIIE,
           TableauHSPARKGLRKLobattoIIIAB,
           TableauHSPARKGLRKLobattoIIIBA,
           TableauHSPARKGLRKLobattoIIICC̄,
           TableauHSPARKGLRKLobattoIIIC̄C,
           TableauHSPARKGLRKLobattoIIID,
           TableauHSPARKGLRKLobattoIIIE

    include("tableaus/tableaus_hspark_secondary.jl")

    export TableauSLRKLobattoIIIAB,
           TableauSLRKLobattoIIIBA,
           TableauSLRKLobattoIIICC̄,
           TableauSLRKLobattoIIIC̄C,
           TableauSLRKLobattoIIID,
           TableauSLRKLobattoIIIE

    include("tableaus/tableaus_slrk.jl")

    export TableauLieA,
           TableauLieB,
           TableauMarchuk,
           TableauMcLachlan2,
           TableauMcLachlan4,
           TableauStrang,
           TableauStrangA,
           TableauStrangB,
           TableauSuzukiFractal,
           TableauTripleJump

    include("tableaus/tableaus_splitting.jl")

    function __init__()
        add_config(:tab_compensated_summation, true)
    end

end


text:
struct SemDiffProximal{A, B, C, D} <: SemDiff
    algorithm::A
    options::B
    operator_g::C
    operator_h::D
end

SemDiffProximal(;algorithm = ProximalAlgorithms.PANOC(), options = Dict{Symbol, Any}(), operator_g, operator_h = nothing, kwargs...) =
    SemDiffProximal(algorithm, options, operator_g, operator_h)

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemDiffProximal)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
#!/usr/bin/env julia

#=
-----------------
Utility Functions
-----------------
=#
"""
    trapz(y, x)

Array based trapezoidal integration of y with respect to x.

Limits of integration defined by the start and end points of the arrays.
"""
function trapz(y::Vector{RheoFloat}, x::Vector{RheoFloat})

    n = length(x)

    @assert n==length(y) "X and Y array length must match."
    n==1 && return zero(RheoFloat)

    r = zero(RheoFloat)
    @inbounds for i in 2:n
        r += (y[i-1] + y[i])*(x[i] - x[i-1])
    end

    r/2

end

"""
    derivCD(y, x)

Given two arrays of data, x and y, calculate dy/dx using central difference
method and forward and backward difference for array boundaries.
"""
function derivCD(y::Vector{RheoFloat}, x::Vector{RheoFloat})

    # get length
    N = length(x)

    # assert y and x arrays are same length
    @assert length(y)==N "X and Y Array lengths must match."

    # initialise zero array of length y
    ydot = similar(y)

    # assume 'imaginary' previous point is 0.0, and Δx is the same as the next one ahead
    # this is a physical assumption that material is at rest before first data point.
    # Could be problematic in some cases if sudden jump as we are actually missing
    # important information about how quickly that jump happened.
    @inbounds ydot[1] = y[1]/(x[2] - x[1])

    # central difference with uneven spacing for general case of constant or variable sample rate
    @inbounds for i in 2:(N-1)
        Δx₁ = x[i] - x[i-1]
        Δx₂ = x[i+1] - x[i]
        ydot[i] = (y[i+1]*Δx₁^2 + (Δx₂^2 - Δx₁^2)*y[i] - y[i-1]*Δx₂^2)/(Δx₁*Δx₂*(Δx₁ + Δx₂))
    end

    # 1st order backwards difference for last element
    ydot[N] = (y[N] - y[N-1])/(x[N] - x[N-1])

    return ydot

end

"""
    derivBD(y, x)

Given two arrays of data, x and y, calculate dy/dx using 1st order
backward difference. Assumes y==0 at a previous point, i.e.
y is 'at rest'. Captures instantaneous loading where derivCD will smooth.
"""
function derivBD(y::Vector{RheoFloat}, x::Vector{RheoFloat})

    # get length
    N = length(x)

    # assert y and x arrays are same length
    @assert length(y)==N "X and Y Array lengths must match."

    # initialise zero array of length y
    ydot = similar(y)

    # assume 'imaginary' previous point is 0.0, and Δx is the same as the next one ahead
    # this is a physical assumption that material is at rest before first data point.
    # Could be problematic in some cases if sudden jump as we are actually missing
    # important information about how quickly that jump happened.
    @inbounds ydot[1] = y[1]/(x[2] - x[1])

    # backwards difference method for rest of points
    @inbounds for i in 2:N
        ydot[i] = (y[i] - y[i-1])/(x[i] - x[i-1])
    end

    return ydot

end

"""
    doublederivCD(y, x)

Given two arrays of data, x and y, calculate d^2(y)/dx^2 using 2nd order difference. Data
should conform to the following constraints:
- No step loading (signal should start at 0)
"""
function doublederivCD(y::Vector{RheoFloat}, x::Vector{RheoFloat})
    # get length
    N = length(x)

    # assert y and x arrays are same length
    @assert length(y)==N "X and Y Array lengths must match."

    # initialise zero array of length y
    yddot = similar(y)

    # no physical assumptions made, unlike for first order BD derivative.
    # forward difference for first element, backwards difference for last element
    @inbounds yddot[1] = (y[3] - 2*y[2] + y[1])/(x[2] - x[1])^2
    @inbounds yddot[N] = (y[N] - 2*y[N-1] + y[N-2])/(x[N] - x[N-1])^2

    # central difference method for rest of points
    @inbounds for i in 2:(N-1)
        yddot[i] = (y[i+1] - 2*y[i] + y[i-1])/(x[i] - x[i-1])^2
    end

    return yddot
end

"""
    derivFbasic(α, ydot, yddot, t, dt)

Compute the Caputo fractional derivative using a simple convolution on the integrated-by-parts
form which does not have a singularity in the kernel. Sample rate must be constant, and this
attribute is not checked in the function body for computational efficiency.
"""
function derivFbasic(α, ydot, yddot, t, dt)
    I1 = ydot[1]*t.^(1-α)
    I2 = conv(t.^(1-α), yddot)[1:length(t)]*dt
    return (I1 + I2)/gamma(2 - α)
end

function constantcheck(t::Vector{RheoFloat})
    # get array of backward differences
    diff = t[2:end] - t[1:end-1]
    # check if any element is not approximately equal to 1st element
    check = all(x -> x ≈ diff[1], diff)
end

function getsampleperiod(t::Vector{RheoFloat})
    # check sample rate is constant, otherwise sample period varies
    @assert constantcheck(t) "Sample-rate must be constant"
    # return sample period
    sampleperiod = t[2] - t[1]
end

"""
    closestindex(x, val)

Find the index of the array element closest to val.
"""
function closestindex(x::Vector{T}, val::Real) where {T<:Real}

    # intialise closest match variable, assuming best match is index 1
    ibest = 1

    # diff between value and current element
    dxbest = abs(x[ibest]-val)

    # loop through all elements, looking for smallest difference
    @inbounds for I in eachindex(x)
        dx = abs(x[I]-val)
        if dx < dxbest
            dxbest = dx
            ibest = I
        end
    end

    ibest
end

"""
    closestindices(x, vals)

Uses `closestindex` iteratively to find closest index for each value in `vals` array,
returns array of indices.
"""
closestindices(x::Vector{T1}, vals::Vector{T2}) where {T1<:Real, T2<:Real} = broadcast(closestindex, (x,), vals)

#=
--------------------------------
Preprocessing base functionality
--------------------------------
=#
"""
    getsigma(τ, samplerate)

Generate sigma/std deviation for gaussian smoothing kernel.

Acts as a low pass filter. Information of time scale τ will be half power,
faster will be increasingly cut. Called by smoothgauss function.
"""
function getsigma(τ::Real, samplerate::Real)

    # get freq, reduce to half power and generate gaussian std. deviation
    smoothfreq = 1.0/τ

    sF_halfpower = smoothfreq/sqrt(2.0*log(2.0))

    σ = samplerate/(2.0*π*sF_halfpower)

    RheoFloat(σ)
end

# """
#     fixed_resample(x::Vector{T}, y::Vector{T}, boundaries::Vector{U}, elperiods::Union{Vector{U}, U}) where T<:RheoFloat where U<:Integer
#
# Resample two arrays with new sample rate(s).
# It can either reduce the sample rate by not taking every array element or it can increase the sample
# rate by inserting new elements. If the number of elements per period (elperiods) is a positive number, then it upsamples; if
# it is negative it downsamples.
#
# The function returns the resampled x and y arrays.
# """
# function fixed_resample(x::Vector{T}, y::Vector{T}, boundaries::Vector{U}, elperiods::Union{Vector{U}, U}; includelastel::Bool=true) where T<:RheoFloat where U<:Integer
#
#     # assert correct function signature
#     @assert length(x)==length(y) "X and Y arrays must have same length."
#     @assert length(elperiods)==length(boundaries)-1 "Number of different sample periods must be 1 less than boundaries provided"
#
#     # initialise resampled arrays as empty
#     xᵦ = Vector{RheoFloat}(undef, 0)
#     yᵦ = Vector{RheoFloat}(undef, 0)
#
#     # loop through boundaries
#     for i in 1:length(boundaries)-1
#
#         # upsampling, starts at each element then intepolates up N times
#         if !signbit(elperiods[i])
#             for k in boundaries[i]:(boundaries[i+1]-1)
#                     x_add  = range(x[k], x[k+1]; length=elperiods[i]+1)[1:end-1]
#                     y_add =  y[k] .+ (x_add[1:end] .- x[k]).*(y[k+1] .- y[k])./(x[k+1] .- x[k])
#                     xᵦ = append!(xᵦ, x_add)
#                     yᵦ = append!(yᵦ, y_add)
#              end
#
#         # downsampling, simply takes every N element as in downsample function
#         elseif signbit(elperiods[i])
#             append!(xᵦ,x[boundaries[i]:abs(elperiods[i]):(boundaries[i+1]-1)])
#             append!(yᵦ,y[boundaries[i]:abs(elperiods[i]):(boundaries[i+1]-1)])
#         end
#     end
#
#     # last element may be missed, so check and add if needed
#     lastel = boundaries[end]
#     if includelastel && xᵦ[end]<x[lastel]
#         append!(xᵦ,x[lastel])
#         append!(yᵦ,y[lastel])
#     end
#
#     return xᵦ, yᵦ
# end

#=
-----------------------------------------
Fitting and predicting base functionality
-----------------------------------------
=#
function singularitytest(modulus; t1::RheoFloat=zero(RheoFloat))

    startval = modulus(t1)

    if isnan(startval) || startval == Inf
        return true
    else
        return false
    end

end

"""
    boltzintegral_nonsing(modulus, time_series, prescribed_dot)

Calculate Boltzmann Superposition integral using direct integration method.

This is much slower and slightly less accurate (depending on sample resolution)
than the convolution method. However, it works for variable sample rate.

# Arguments

- `modulus`: Viscoelastic modulus function in which the parameters have been substituted
- `time_series`: The array of times
- `prescribed_dot`: Derivative of (usually prescribed) variable inside the integration kernel
"""
function boltzintegral_nonsing(modulus, time_series::Vector{Float64}, prescribed_dot::Vector{Float64})

    # need to add an additional 'previous' time point to capture any instantaneous loading
    time_previous = time_series[1] - (time_series[2] - time_series[1])
    time_mod = vcat([time_previous], time_series)
    # material is assumed at rest at this 'previous' time
    prescribed_dot_mod = vcat([0.0], prescribed_dot)

    I = zeros(length(time_mod))
    @inbounds for (i,v) in enumerate(time_mod)
        # generate integral for each time step
        τ = time_mod[1:i]
        Modulus_arg = v .- τ
        Modulusᵢ = modulus(Modulus_arg)
        df_dtᵢ = prescribed_dot_mod[1:i]
        intergrand = Modulusᵢ.*df_dtᵢ

        I[i] = trapz(intergrand, τ)
    end
    # fix initial point
    I[2] = (prescribed_dot[1]*modulus(time_series)*(time_series[2] - time_series[1]))[1]
    I[2:end]

end

"""
    obj_var_nonsing(params, grad, modulus, time_series, prescribed_dot, measured; _insight::Bool = false)

Generate the sum-of-squares of the difference between `measured` and the Boltzmann
convolution integral of the selected `modulus` and `prescribed_dot`. Used when
sample rate is variable and no singularity in model.

# Arguments

- `params`: Array of parameters sent to `modulus`
- `grad`: Gradient argument used by NLOpt
- `modulus`: Viscoelastic modulus function
- `time_series`: Array of time data
- `prescribed_dot`: Convolved with `modulus` over `time_series`, usually dϵ/dt for stress relaxation and dσ/dt for creep
- `measured`: Data for comparison against, usually σ for stress relaxation and ϵ for creep
- `_insight`: Declare whether insight info should be shown when this function is called, true or false
"""
function obj_var_nonsing(params, grad, modulus, time_series, prescribed_dot, measured; _insight::Bool = false)

    _insight && println("Current Parameters: ", params)

    mod = (t->modulus(t,params))
    convolved = boltzintegral_nonsing(mod, time_series, prescribed_dot)

    cost = sum((measured - convolved).^2)

end

"""
    boltzintegral_sing(modulus, time_series, prescribed_dot)

Calculate Boltzmann Superposition integral using direct integration method.

This is much slower and slightly less accurate (depending on sample resolution)
than the convolution method. However, it works for variable sample rate.

# Arguments

- `modulus`: Viscoelastic modulus function in which the parameters have been substituted
- `time_series`: The array of times
- `prescribed_dot`: Derivative of (usually prescribed) variable inside the integration kernel
"""
function boltzintegral_sing(modulus, time_series::Vector{Float64}, prescribed_dot::Vector{Float64})

    offset_t0 = (time_series[2] - time_series[1])/singularity_offset
    # need to add an additional 'previous' time point to capture any instantaneous loading
    time_previous = time_series[1] - (time_series[2] - time_series[1])
    time_mod = vcat([time_previous], time_series)
    # material is assumed at rest at this 'previous' time
    prescribed_dot_mod = vcat([0.0], prescribed_dot)

    offset = offset_t0
    I = zeros(length(time_mod))
    for (i,v) in enumerate(time_mod)
        if i>1
            offset = (time_mod[i] - time_mod[i-1])/singularity_offset
        end
        τ = time_mod[1:i]
        Modulus_arg = v .- τ
        Modulus_arg[end] = offset
        Modulusᵢ = modulus(Modulus_arg)
        df_dtᵢ = prescribed_dot_mod[1:i]

        intergrand = Modulusᵢ.*df_dtᵢ

        I[i] = trapz(intergrand, τ)
    end
    # fix initial point
    I[2] = (prescribed_dot[1]*modulus([offset_t0])*(time_series[2] - time_series[1]))[1]
    I[2:end]

end

"""
    obj_var_sing(params, grad, modulus, time_series, prescribed_dot, measured; _insight::Bool = false)

Generate the sum-of-squares of the difference between `measured` and the Boltzmann
convolution integral of the selected `modulus` and `prescribed_dot`. Used when
sample rate is variable and there IS singularity in model.

# Arguments

- `params`: Array of parameters sent to `modulus`
- `grad`: Gradient argument used by NLOpt
- `modulus`: Viscoelastic modulus function
- `time_series`: Array of time data
- `prescribed_dot`: Convolved with `modulus` over `time_series`, usually dϵ/dt for stress relaxation and dσ/dt for creep
- `measured`: Data for comparison against, usually σ for stress relaxation and ϵ for creep
- `_insight`: Declare whether insight info should be shown when this function is called, true or false
"""
function obj_var_sing(params, grad, modulus, time_series, prescribed_dot, measured; _insight::Bool = false)

    _insight && println("Current Parameters: ", params)

    mod = (t->modulus(t,params))
    convolved = boltzintegral_sing(mod, time_series, prescribed_dot)

    # skip first point for error computation, workaround for singularity approximation error
    cost = sum((measured[2:end] - convolved[2:end]).^2)

end

"""
    boltzconvolve(modulus, time_series, dt::RheoFloat,prescribed_dot)

Calculate Boltzmann Superposition integral using convolution method.

This is much faster and slightly more accurate (depending on sample resolution)
than the integral method. However, it only works for constant sample rate.

# Arguments

- `modulus`: Viscoelastic modulus function in which the parameters have been substituted
- `time_series`: The array of times
- `dt`: Constant time step (sample period)
- `prescribed_dot`: Derivative of (usually prescribed) variable inside the integration kernel
"""
function boltzconvolve(modulus, time_series, dt, prescribed_dot)

    Modulus = modulus(time_series)
    # fast convolution
    β = conv(Modulus, prescribed_dot)
    # pick out relevant elements (1st half) and multiply by dt
    β = β[1:length(time_series)]*dt

end

"""
    obj_const_nonsing(params, grad, modulus, time_series, dt, prescribed_dot, measured; _insight::Bool = false)

Generate the sum-of-squares of the difference between `measured` and the Boltzmann
convolution integral of the selected `modulus` and `prescribed_dot`. Used when
sample rate is constant and model does not feature singularity.

# Arguments

- `params`: Array of parameters sent to `modulus`
- `grad`: Gradient argument used by NLOpt
- `modulus`: Viscoelastic modulus function
- `time_series`: Array of time data
- `dt`: Constant time step (sample period)
- `prescribed_dot`: Convolved with `modulus` over `time_series`, usually dϵ/dt for stress relaxation and dσ/dt for creep
- `measured`: Data for comparison against, usually σ for stress relaxation and ϵ for creep
- `_insight`: Declare whether insight info should be shown when this function is called, true or false
"""
function obj_const_nonsing(params, grad, modulus, time_series, dt, prescribed_dot, measured; _insight::Bool = false)

    _insight && println("Current Parameters: ", params)

    # Replace parameters inside the viscoelastic modulus
    mod = (t->modulus(t,params))
    convolved = boltzconvolve(mod, time_series, dt, prescribed_dot)

    cost = sum((measured - convolved).^2)
    return cost

end

"""
    obj_const_sing(params, grad, modulus, time_series, dt, prescribed_dot, measured; _insight::Bool = false)

Generate the sum-of-squares of the difference between `measured` and the Boltzmann
convolution integral of the selected `modulus` and `prescribed_dot`. Used when
sample rate is constant and model does feature singularity.

# Arguments

- `params`: Array of parameters sent to `modulus`
- `grad`: Gradient argument used by NLOpt
- `modulus`: Viscoelastic modulus function
- `time_series`: Array of time data
- `dt`: Constant time step (sample period)
- `prescribed_dot`: Convolved with `modulus` over `time_series`, usually dϵ/dt for stress relaxation and dσ/dt for creep
- `measured`: Data for comparison against, usually σ for stress relaxation and ϵ for creep
- `_insight`: Declare whether insight info should be shown when this function is called, true or false
"""
function obj_const_sing(params, grad,modulus, time_series,dt, prescribed_dot, measured; _insight::Bool = false)

    _insight && println("Current Parameters: ", params)

    mod = (t->modulus(t,params))
    # convolved = boltzconvolve_sing(modulus, time_series, dt, params, prescribed_dot)
    convolved = boltzconvolve(mod, time_series, dt, prescribed_dot)

    # skip first point for error computation, workaround for singularity approximation error
    cost = sum((measured[2:end] - convolved[2:end]).^2)

end

"""
    obj_const_weighted(params, grad, modulus, time_series, dt, prescribed_dot, measured_weighted, weights; _insight::Bool = false)

Post-convolution weighted constant sample-rate cost function for emphasising regions without introducing variable sample-rates.
"""
function obj_const_weighted(params, grad, modulus, time_series, dt, prescribed_dot, measured_weighted, weights; _insight::Bool = false)

    _insight && println("Current Parameters: ", params)

    # Replace parameters inside the viscoelastic modulus
    mod = (t->modulus(t,params))
    convolved = boltzconvolve(mod, time_series, dt, prescribed_dot)

    cost = sum((measured_weighted - convolved[weights]).^2)
    return cost
end

"""
    leastsquares_init(params_init, low_bounds, hi_bounds, modulus, time_series, dt, prescribed_dot, measured; insight = false, sampling = "constant", singularity = false)

Initialise then begin a least squares fitting of the supplied data.

# Arguments

- `params_init`: Initial parameters to be used (starting guess)
- `low_bounds`: Lower bounds for parameters
- `hi_bounds`: Higher bounds for parameters
- `modulus`: Viscoelastic modulus function
- `time_series`: Array of time data
- `dt`: Constant time step (sample period)
- `prescribed_dot`: Convolved with `modulus` over `time_series`, usually dϵ/dt for stress relaxation and dσ/dt for creep
- `measured`: Data for comparison against, usually σ for stress relaxation and ϵ for creep
- `constant_sampling`: true if sample rate is constant, false otherwise
- `insight`: Declare whether insight info should be shown when this function is called, true or false
- `singularity`: Presence of singularity in model
- `indweight`: indices weighting, see `modelfit` for more discussion
"""
function leastsquares_init(params_init::Vector{RheoFloat}, low_bounds::RheovecOrNone,
                           hi_bounds::RheovecOrNone, modulus,
                           time_series::Vector{RheoFloat}, dt::RheoFloat,
                           prescribed_dot::Vector{RheoFloat}, measured::Vector{RheoFloat};
                           insight::Bool = false, constant_sampling::Bool=true,
                           singularity::Bool = false, _rel_tol = 1e-4, indweights=nothing)

    # initialise NLOpt.Opt object with :LN_SBPLX Subplex algorithm
    opt = Opt(:LN_SBPLX, length(params_init))

    # set lower bounds and upper bounds unless they take null value
    if !isnothing(low_bounds)
        low_bounds = convert(Vector{Float64},low_bounds)
        lower_bounds!(opt, low_bounds)
    end

    if !isnothing(hi_bounds)
        hi_bounds = convert(Vector{Float64}, hi_bounds)
        upper_bounds!(opt, hi_bounds)
    end

    # set relative tolerance
    xtol_rel!(opt, _rel_tol)

    # Convert to float64 to avoid conversion by NLOpt
    params_init = convert(Vector{Float64},params_init)
    time_series = convert(Vector{Float64},time_series)
    prescribed_dot = convert(Vector{Float64},prescribed_dot)
    measured = convert(Vector{Float64},measured)
    dt = convert(Float64, dt)

    # set Opt object as a minimisation objective. Use a closure for additional
    # arguments sent to object objectivefunc
    if constant_sampling && !singularity && isnothing(indweights)
        min_objective!(opt, (params, grad) -> obj_const_nonsing(params, grad, modulus,
                                                            time_series, dt,
                                                            prescribed_dot, measured;
                                                            _insight = insight))

    elseif constant_sampling && singularity && isnothing(indweights)
        # remove singularity, just go close to it, 1/10th over first sample period
        time_series[1] = 0.0 + (time_series[2] - time_series[1])/singularity_offset
        min_objective!(opt, (params, grad) -> obj_const_sing(params, grad, modulus,
                                                        time_series, dt,
                                                        prescribed_dot, measured;
                                                        _insight = insight))

    elseif constant_sampling && !singularity && !isnothing(indweights)
        measured_weighted = measured[indweights]
        min_objective!(opt, (params, grad) -> obj_const_weighted(params, grad, modulus,
                                                            time_series, dt,
                                                            prescribed_dot, measured_weighted, indweights;
                                                            _insight = insight))

    elseif constant_sampling && singularity && !isnothing(indweights)
        indsingular = indweights[indweights.>1]
        measured_weighted = measured[indsingular]
        # remove singularity, just go close to it, 1/10th over first sample period
        time_series[1] = 0.0 + (time_series[2] - time_series[1])/singularity_offset
        min_objective!(opt, (params, grad) -> obj_const_weighted(params, grad, modulus,
                                                        time_series, dt,
                                                        prescribed_dot, measured_weighted, indsingular;
                                                        _insight = insight))

    elseif !constant_sampling && !singularity
        min_objective!(opt, (params, grad) -> obj_var_nonsing(params, grad, modulus,
                                                        time_series, prescribed_dot,
                                                        measured; _insight = insight))

    elseif !constant_sampling && singularity
        @warn "Note that large changes in sample rate, particularly from low sample rate to
        much higher sample rate can introduce significant innacuracies due to the singularity
        approximation method being used. To avoid this, please resample your data to a constant
        sample rate before fitting."
        min_objective!(opt, (params, grad) -> obj_var_sing(params, grad, modulus,
                                                        time_series, prescribed_dot,
                                                        measured; _insight = insight))

    end

    # minimise objective func, minx are the parameters resulting in minimum
    (minf, minx, ret) = NLopt.optimize(opt, params_init)

    # return all
    return (convert(RheoFloat,minf), convert(Vector{RheoFloat},minx), ret)

end

#=
----------------------------------------------
Step fitting and predicting base functionality
----------------------------------------------
=#
function obj_step_nonsing(params, grad, modulus, t, prescribed::Float64, measured::Vector{Float64}; _insight=false)
    if _insight
        println("Current Parameters: ", params)
    end

    mod = (t->modulus(t,params))
    estimated = prescribed*mod(t)

    cost = sum((measured - estimated).^2)
end

function obj_step_weighted(params, grad, modulus, t, prescribed::Float64, measured::Vector{Float64}, weights; _insight=false)
    if _insight
        println("Current Parameters: ", params)
    end

    mod = (t->modulus(t,params))
    estimated = prescribed*mod(t)

    cost = sum((measured - estimated[weights]).^2)
end

function leastsquares_stepinit(params_init::Vector{RheoFloat}, low_bounds::RheovecOrNone,
                           hi_bounds::RheovecOrNone, modulus,
                           time_series::Vector{RheoFloat}, prescribed::RheoFloat,
                           measured::Vector{RheoFloat}; insight::Bool = false,
                           singularity::Bool = false, _rel_tol = 1e-4, indweights=nothing)

    # initialise NLOpt.Opt object with :LN_SBPLX Subplex algorithm
    opt = Opt(:LN_SBPLX, length(params_init))

    # set lower bounds and upper bounds unless they take null value
    if !isnothing(low_bounds)
        low_bounds = convert(Vector{Float64}, low_bounds)
        lower_bounds!(opt, low_bounds)
    end

    if !isnothing(hi_bounds)
        hi_bounds = convert(Vector{Float64}, hi_bounds)
        upper_bounds!(opt, hi_bounds)
    end

    # set relative tolerance
    xtol_rel!(opt, _rel_tol)

    # set Opt object as a minimisation objective. Use a closure for additional
    # arguments sent to object objectivefunc

    params_init = convert(Vector{Float64},params_init)
    time_series = convert(Vector{Float64},time_series)
    prescribed = convert(Float64,prescribed)
    measured = convert(Vector{Float64},measured)

    if isnothing(indweights) && !singularity
        min_objective!(opt, (params, grad) -> obj_step_nonsing(params, grad, modulus,
                                                            time_series, prescribed, measured;
                                                            _insight = insight))

    elseif isnothing(indweights) && singularity
        indices = collect(Integer, 2:length(measured))
        measured_weighted = measured[indices]
        min_objective!(opt, (params, grad) -> obj_step_weighted(params, grad, modulus,
                                                        time_series, prescribed, measured_weighted, indices;
                                                        _insight = insight))

    elseif !isnothing(indweights) && !singularity
        measured_weighted = measured[indweights]
        min_objective!(opt, (params, grad) -> obj_step_weighted(params, grad, modulus,
                                                            time_series, prescribed, measured_weighted, indweights;
                                                            _insight = insight))

    elseif !isnothing(indweights) && singularity
        indsingular = indweights[indweights.>1]
        measured_weighted = measured[indsingular]
        min_objective!(opt, (params, grad) -> obj_step_weighted(params, grad, modulus,
                                                        time_series, prescribed, measured_weighted, indsingular;
                                                        _insight = insight))

    end

    # minimise objective func, minx are the parameters resulting in minimum
    (minf, minx, ret) = NLopt.optimize(opt, params_init)

    # return all
    return (convert(RheoFloat,minf), convert(Vector{RheoFloat},minx), ret)

end


text:
using Random, Plots, Distributions, Measures; pyplot()

N = 5*10^3
alpha = 0.05
confLevel = 1 - alpha
z = quantile(Normal(),1-alpha/2)

function randCI(n,p)
    sample = rand(n) .< p
    pHat = sum(sample)/n
    serr = sqrt(pHat*(1-pHat)/n)
    (pHat - z*serr, pHat + z*serr)
end
cover(p,ci) = ci[1] <= p && p <= ci[2]

pGrid = 0.1:0.01:0.9
nGrid = 5:1:50
errs = zeros(length(nGrid),length(pGrid))

for i in 1:length(nGrid)
    for j in 1:length(pGrid)
        Random.seed!(0)
        n, p = nGrid[i], pGrid[j]
        coverageRatio = sum([cover(p,randCI(n,p)) for _ in 1:N])/N
        errs[i,j] = confLevel - coverageRatio
    end
end

default(xlabel = "p", ylabel = "n",
    xticks =([1:5:length(pGrid);], minimum(pGrid):0.05:maximum(pGrid)),
    yticks =([1:5:length(nGrid);], minimum(nGrid):5:maximum(nGrid)))

p1 = heatmap(errs, c=cgrad([:white, :black]))
p2 = heatmap(abs.(errs) .<= 0.04, legend = false, c=cgrad([:black, :white]))
plot(p1,p2, size = (1000,400), margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using AoC
using Test

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


text:
#
# Implementation of AtomsBase interface in a struct-of-arrays style.
#
export FastSystem

struct FastSystem{D, L <: Unitful.Length, M <: Unitful.Mass} <: AbstractSystem{D}
    box::SVector{D, SVector{D, L}}
    boundary_conditions::SVector{D, BoundaryCondition}
    positions::Vector{SVector{D, L}}
    atomic_symbols::Vector{Symbol}
    atomic_numbers::Vector{Int}
    atomic_masses::Vector{M}
end

# Constructor to fetch the types
function FastSystem(box, boundary_conditions, positions, atomic_symbols, atomic_numbers, atomic_masses)
    FastSystem{length(box),eltype(eltype(positions)),eltype(atomic_masses)}(
        box, boundary_conditions, positions, atomic_symbols, atomic_numbers, atomic_masses
    )
end

# Constructor to take data from another system
function FastSystem(system::AbstractSystem)
    FastSystem(bounding_box(system), boundary_conditions(system), position(system),
               atomic_symbol(system), atomic_number(system), atomic_mass(system))
end

# Convenience constructor where we don't have to preconstruct all the static stuff...
function FastSystem(particles, box, boundary_conditions)
    D = length(box)
    if !all(length.(box) .== D)
        throw(ArgumentError("Box must have D vectors of length D=$D."))
    end
    if length(boundary_conditions) != D
        throw(ArgumentError("Boundary conditions should be of length D=$D."))
    end
    if !all(n_dimensions.(particles) .== D)
        throw(ArgumentError("Particles must have positions of length D=$D."))
    end
    FastSystem(box, boundary_conditions, position.(particles), atomic_symbol.(particles),
               atomic_number.(particles), atomic_mass.(particles))
end

function Base.show(io::IO, system::FastSystem)
    print(io, "FastSystem")
    show_system(io, system)
end

bounding_box(sys::FastSystem)        = sys.box
boundary_conditions(sys::FastSystem) = sys.boundary_conditions
Base.length(sys::FastSystem)         = length(sys.positions)
Base.size(sys::FastSystem)           = size(sys.positions)

species_type(sys::FS) where {FS <: FastSystem} = AtomView{FS}
Base.getindex(sys::FastSystem, index::Int)     = AtomView(sys, index)

position(s::FastSystem)       = s.positions
atomic_symbol(s::FastSystem)  = s.atomic_symbols
atomic_number(s::FastSystem)  = s.atomic_numbers
atomic_mass(s::FastSystem)    = s.atomic_masses
velocity(s::FastSystem)       = missing

position(s::FastSystem, i)      = s.positions[i]
atomic_symbol(s::FastSystem, i) = s.atomic_symbols[i]
atomic_number(s::FastSystem, i) = s.atomic_numbers[i]
atomic_mass(s::FastSystem, i)   = s.atomic_masses[i]
velocity(s::FastSystem, i)      = missing


text:
function transpose_11!(A)
    n,m = size(A)
    B=Matrix{eltype(A)}(undef,m,n)
    for j in 1:n
        for i in 1:m
            B[i,j] = A[j,i]
        end
    end
    A=B
    return A
end

function transpose_12!(A)
    n,m = size(A)
    for j in 1:n
        for i in j:m
            temp = A[i,j]
            A[i,j] = A[j,i]
            A[j,i] = temp
        end
    end
    return A
end


text:
export
    BicycleModel

immutable BicycleModel
    def::VehicleDef
    a::Float64 # distance between cg and front axle [m]
    b::Float64 # distance between cg and rear axle [m]
end
function BicycleModel(def::VehicleDef;
    a::Float64 = 1.5,
    b::Float64 = 1.5,
    )

    return BicycleModel(def, a, b)
end


text:
module ACE

using Base: NamedTuple
using Reexport

# external imports that are useful for all submodules
include("imports.jl")

@extimports
@baseimports

# TODO
# - move to imports

using ACEbase.ObjectPools: acquire!, release!, VectorPool
using ForwardDiff: derivative
import ChainRules: rrule, ZeroTangent, NoTangent
import ACEbase: evaluate, evaluate_d
import ACEbase: gradtype, valtype
import ACEbase: acquire_B!, release_B!, acquire_dB!, release_dB!

# TODO: gradtype should have a standard fallback

abstract type AbstractACEModel end

abstract type AbstractProperty end

function coco_init end
function coco_zeros end
function coco_filter end
function coco_dot end
function coco_type end

# TODO
# * decide on rand(basis) interface

# * move these the following definitions to ACEbase
function _rrule_evaluate end
function _rrule_evaluate_d end

include("auxiliary.jl")

include("rotations3d.jl")
using ACE.Rotations3D

include("polynomials/wigner.jl")

include("states.jl")
include("symmetrygroups.jl")
include("properties.jl")

include("prototypes.jl")

# basic polynomial building blocks
include("polynomials/sphericalharmonics.jl")
include("polynomials/transforms.jl"); @reexport using ACE.Transforms
include("polynomials/orthpolys.jl"); @reexport using ACE.OrthPolys

# The One-particle basis is the first proper building block
include("oneparticlebasis.jl")

# three specific 1p-bases that are always useful
include("Ylm1pbasis.jl")
include("Rn1pbasis.jl")
include("scal1pbasis.jl")
include("discrete1pbasis.jl")

include("product_1pbasis.jl")

# basis selectors used to specify finite subsets of basis functions
include("basisselectors.jl")
# ... amongst other things used to initialize sparse basis sets
include("sparsegrids.jl")

# the permutation-invariant, and symmerized bases
include("pibasis.jl")
include("symmbasis.jl")

# some experimental stuff
include("multiplier.jl")
include("bonds.jl")

# models and model evaluators

include("linearmodel.jl")

include("evaluator.jl")
# include("grapheval.jl")

include("random.jl")
@reexport using ACE.Random

include("utils.jl")
@reexport using ACE.Utils

include("testing/testing.jl")

include("ad.jl")

# ---------------- some extra experimental dispatching

evaluate(basis::SymmetricBasis, Xs::AbstractVector) =
      evaluate(basis, ACEConfig(Xs))

evaluate_d(basis::SymmetricBasis, Xs::AbstractVector) =
      evaluate_d(basis, ACEConfig(Xs))

evaluate(model::LinearACEModel, Xs::AbstractVector) =
      evaluate(model, ACEConfig(Xs))

grad_config(model::LinearACEModel, Xs::AbstractVector) =
      grad_config(model, ACEConfig(Xs))

grad_params(model::LinearACEModel, Xs::AbstractVector) =
      grad_params(model, ACEConfig(Xs))

end # module

# TODO:
# # orthogonal basis
# include("orth.jl")

# include("utils/importv5.jl")
# include("compat/compat.jl")
# include("export/export.jl")
# - bond model
# - pure basis
# - real basis
# - regularisers
# - descriptors
# - random potentials
# TODO -> move elsewhere!!!
# include("rpi/rpi_degrees.jl")


text:
function soft_coulomb(x; Z=1.0)
    @. Z / sqrt(x^2 + 1.0)
end

function harmonic(x; k=1.0)
    @. 0.5 * k * x^2
end


text:
function CliqueExpansion(H::hypergraph,weighted::Bool=true,binary::Bool=false)
    """
    Weighted clique expansion where a hyperedge e is expanded to a
    weighted clique with each edge having weight 1/(|e| - 1)
    """
    n = length(H.D)
    I = Int64[]
    J = Int64[]
    V = Float64[]
    ks = setdiff(keys(H.E),1)
    for k in ks
        for edge in keys(H.E[k])
            weight = H.E[k][edge]
            for i = 1:k-1
                ei = edge[i]
                for j = i+1:k
                    ej = edge[j]
                    push!(I,ei)
                    push!(J,ej)
                    if weighted
                        push!(V, weight / (k - 1))
                    else
                        push!(V, weight)
                    end
                end
            end
        end
    end
    A = SparseArrays.sparse(I,J,V,n,n)
    for i = 1:n
        A[i, i] = 0.0
    end
    SparseArrays.dropzeros!(A)
    A = SparseArrays.sparse(A+A')
    if binary
        I, J, V = SparseArrays.findnz(A)
        A = SparseArrays.sparse(I, J, 1, n, n)
    end
    return A
end

function CliqueExpansionModularity(H::hypergraph,gamma::Float64=1.0,weighted::Bool=true,randflag::Bool=false,binary::Bool=false,clusterpenalty::Float64=0.0,maxits::Int64=10000)
    return CliqueExpansionModularity(H,gamma;weighted=weighted,randflag=randflag,binary=binary,clusterpenalty=clusterpenalty,maxits=maxits)
end

function CliqueExpansionModularity(H::hypergraph,gamma::Float64=1.0;weighted::Bool=true,randflag::Bool=false,binary::Bool=false,clusterpenalty::Float64=0.0,maxits::Int64=10000)
    """
    Perform a clique expansion on the hypergraph H and then run vanilla
    modularity on the resulting graph.
    """
    A = CliqueExpansion(H,weighted,binary)
    return VanillaModularity(A,gamma,randflag,clusterpenalty,maxits)
end

function StarExpansionModularity(H::hypergraph,gamma::Float64=1.0;weighted::Bool=true,randflag::Bool=false,binary::Bool=false,maxits::Int64=100)
    """
    Perform a clique expansion on the hypergraph H and then run vanilla
    modularity on the resulting graph.
    """
    He2n, w = hypergraph2incidence(H)
    m,n = size(He2n)
    A = [spzeros(n,n) He2n'; He2n spzeros(m,m) ]
    Za = VanillaModularity(A,gamma,randflag,maxits)
    return Za[1:n]
end

function VanillaModularity(A::SparseArrays.SparseMatrixCSC{Float64,Int64},gamma::Float64=1.0,randflag::Bool=false,clusterpenalty::Float64=0.0,maxits::Int64=10000)
    """
    Vanilla modularity algorithm, obtained by calling the LambdaLouvain algorithm
    implementation from:

    Parameterized Correlation Clustering in Hypergraphs and Bipartite Graphs
    https://arxiv.org/abs/2002.09460

    Code: https://github.com/nveldt/ParamCC/blob/master/src/Graph_Louvain.jl
    """

    d = vec(sum(A,dims = 2))
    n = length(d)
    vol = sum(d)
    lam = gamma/vol
    Cs = LambdaLouvain(A,d,lam,randflag,maxits,clusterpenalty)

    c = Cs[:,end]
    @assert(length(c) == n)

    return c
end

function computeDyadicResolutionParameter(H, Z; mode = "γ", weighted=true, binary=false, pseudocount = 0.0)
    """
    compute the dyadic resolution parameter associated to a partition using the formula from Newman (2016): https://arxiv.org/abs/1606.02319
    """

    G = CliqueExpansion(H, weighted, binary)
    I, J = SparseArrays.findnz(G)

    n = length(H.D)  # number of nodes
    m = sum(G)/2     # number of edges

    # form degree sequence and edge counts
    D = vec(sum(G, dims=1))

    m_in = pseudocount
    m_out = pseudocount

    for k in 1:length(I)
        if Z[I[k]] == Z[J[k]]
            m_in  += G[I[k], J[k]]/2
        else
            m_out += G[I[k], J[k]]/2
        end
    end

    # compute resolution parameter
    V = [sum(D[findall(Z .== c)]) for c in unique(Z)]
    ωᵢ = 4*m*m_in / (sum(V.^2))
    ωₒ = (2m_out)/(2m - (sum(V.^2)/(2m)))

    if mode == "γ"
        γ = (ωᵢ - ωₒ)/(log(ωᵢ) - log(ωₒ))
        return γ
    else
        return(ωᵢ, ωₒ)
    end
end

function dyadicModularity(H, Z, γ; weighted=false, binary=false)
    """
    Compute the normalized dyadic modularity a given partition in the graph obtained from projecting a hypergraph.
    The resolution parameter and basic properties of the dyadic graph projection can be specified.

    # Arguments

    - H::hypergraph, a hypergraph
    - Z::Vector{Int4}, a label vector
    - γ::Float64, the resolution parameter. Must be strictly positive.
    - weighted::bool, whether to weight the clique obtained by projecting a hypergraph according to a factor of 1/(k-1), where k is the clique size. Choosing to do so preserves the degrees of nodes in the dyadic projection. If false and if binary is false, then each edge within the clique has weight 1.
    - binary::bool, whether to convert the graph to a binary graph in which present edges have weight 1 and absent edges have weight 0. If false and if weighted is false, then the weight of an edge between two nodes in the return value is equal to the number of hyperedges containing both of those nodes.

    # Returns

    - Q::Float64, a normalized modularity score between -1 and 1.

    """
    G = CliqueExpansion(H, weighted, binary)
    d = vec(sum(G, dims=1))

    # non-degree (cut) term
    edge_obj = 0.0
    for (i, j, v) in zip(SparseArrays.findnz(G)...)
        if Z[i] == Z[j]
            edge_obj += v
        end
    end

    # volume terms
    vols = Dict{Int64, Float64}()
    for c in unique(Z)
        vols[c] = 0.0
    end
    for i = 1:length(d)
        vols[Z[i]] += d[i]
    end

    Q = edge_obj
    volG = sum(d)
    vol_term = 0.0
    for c in unique(Z)
        Q -= γ * vols[c]^2 / volG
    end

    return Q / volG
end

function dyadicLogLikelihood(H, Z, ωᵢ, ωₒ; weighted=false, binary=false, constants = false)
    G = CliqueExpansion(H, weighted, binary)
    d = vec(sum(G, dims=1))

    # Eq. (14) from https://arxiv.org/pdf/1606.02319.pdf
    γ = (ωᵢ - ωₒ) / (log(ωᵢ) - log(ωₒ))
    Q = dyadicModularity(H, Z, γ; weighted=weighted)
    m = sum(d) / 2
    B = m * log(ωᵢ / ωₒ)
    C = m * (ωₒ + log(ωₒ))

    if constants
        K = sum(SpecialFunctions.loggamma.(SparseArrays.nonzeros(G).+1)) / 2
        return B * Q + C - K
    end
    return B * Q + C
end

function dyadicLogLikelihoodDirect(H, Z, ωᵢ, ωₒ)
    """
    A naive version of the dyadic log-likelihood calculation, used only for testing the fast one.
    """
    G = CliqueExpansion(H, false, false)
    d = vec(sum(G, dims=1))
    m = sum(d)/2
    L = 0.0
    for (i, j, v) in zip(SparseArrays.findnz(G)...)
        if i < j
            a = G[i,j]
            ω = Z[i] == Z[j] ? ωᵢ : ωₒ
            L += log(poisson_pdf(Int(a), ω*d[i]*d[j]/(2m)))
        end
    end
    return L
end

# Below are functions for just graph louvain, nothing about hypergraphs

function ConstructAdj(C::SparseArrays.SparseMatrixCSC,n::Int64)
    """
    ConstructAdj: Construct Adjacency List
    This takes in a sparse adjacency matrix for a graph, and returns an adjacency
    list. While it seems inefficient to store the matrix in multiple ways, as long
    as there are no terrible memory issues, it is VERY handy to be able to quickly
    get a list of neighbors for each node.

    The function also returns the degree vector d. This is NOT the weighted degree,
        d[i] = total number of neighbors of node i
    """
    rp = C.rowval
    ci = C.colptr
    Neighbs = Vector{Vector{Int64}}()
    d = zeros(Int64,n)
    for i = 1:n
        # chop up the rp vector and put it in Neighbs
        push!(Neighbs,rp[ci[i]:ci[i+1]-1])
        d[i] = ci[i+1]-ci[i]
    end

    # d is the number of neighbors. This is the unweighted degree,
    # but note importantly that if the original graph is weighted this is
    # not the same as the degree vector d we will sometimes use
    return Neighbs, d
end

"""
Evaluates the LambdaCC objective, which is equivalent to modularity in graphs.

https://dl.acm.org/doi/pdf/10.1145/3178876.3186110

A = adjacency matrix for a graph
c = clustering vector
w = weights vector (often taken to be degrees of nodes)
lam = resolution paramter

This can be replaced with our function for checking the hypergraph modularity
objective.
"""
function LamCCobj(A::SparseArrays.SparseMatrixCSC,c,w,lam)

    w_volA = sum(w) # weighted volume
    obj = (lam*(w_volA)^2-lam*sum(w.^2))/2 # constant with respect to the clustering

    for i = 1:maximum(c)
        S = findall(x->x==i,c)
        AS = A[:,S]
        vol = sum(AS.nzval);
        SAS = AS[S,:]
        edges = sum(SAS.nzval);
        cut = vol-edges
        w_volS = sum(w[S])
        obj += 1/2*(cut - lam*w_volS*(w_volA-w_volS))
    end
    return obj
end

"""
Run the Louvain algorithm many times, taking the result with the best objective.
"""
function Many_Louvain(A::SparseArrays.SparseMatrixCSC{Float64,Int64},w::Vector{Float64},lam::Float64,numtimes::Int64,maxits::Int64=10000)

    n = size(A,1)
    BestObj = Inf
    cBest = collect(1:n)

    for k = 1:numtimes

        ######
        Cs = LambdaLouvain(A,w,lam,true,maxits)
        ######

        c = Cs[:,end]

        ######
        obj = LamCCobj(A,c,w,lam)
        ######

        if obj < BestObj
            BestObj = obj
            cBest = c
        end
    end

    return cBest, BestObj
end

"""
The full Louvain algorithm. Lambda is the resolution parameter; think of it as
lambda = gamma/(vol(G)), where gamma is the more traditional resolution parameter
associated with the modularity objective.

w = weights function, w[i] = weight for node i. This is often the degree of node
    i, but doesn't have to be.
"""
function LambdaLouvain(A::SparseArrays.SparseMatrixCSC{Float64,Int64},w::Vector{Float64},lam::Float64,randflag::Bool=false,maxits::Int64=10000,clusterpenalty::Float64=0.0)

    @assert(LinearAlgebra.issymmetric(A))
    n = size(A,1)

    # Step 1: greedy moves until no more improvement
    c, improved = LambdaLouvain_Step(A,w,lam,randflag,maxits,clusterpenalty)

    # Keeps track of the clustering that is found after each call to Step 1.
    if improved
        Cs = c
        c_old = copy(c)
    else
        Cs = c
    end

    # As long as something is still improving each time you call step 1, keep going
    while improved

        # Step 2: Collapse the clustering into supernodes

        Anew, wnew = collapse_clustering(A,w,c_old)

        # Step 1: Go back to greedy local moves, this time on the reduced graph
        cSuper, improved = LambdaLouvain_Step(Anew,wnew,lam,randflag,maxits,clusterpenalty)
        N = length(wnew)    # N = number of supernodes = number clusters from last round

        # Undo the procedure that merged nodes into super nodes, so that you get a clustering vector of length n
        # Do this only if the last call to Step 1 led to at least one greedy move.

        if improved

            # Extract what that new clustering is
            c_new = zeros(Int64,n)

            Clusters = clusters_from_cvec(c_old)

            # For each supernode, place all original node IDs that make it up
            # into a cluster of the supernodes label
            for i = 1:N

                # Get the cluster that supernode i is in
                SuperI_Cluster = cSuper[i]

                # Get individual node ID that are in supernode i.

                SuperI_nodes = Clusters[i]
                c_new[SuperI_nodes] .= SuperI_Cluster
            end
            Cs = [Cs c_new]
            c_old = copy(c_new)
        end
    end

    return Cs
end

"""
From a cluster indicator vector, extract a vector of vectors storing clusters
"""

function clusters_from_cvec(c::Vector{Int64})

    Clusters = Vector{Vector{Int64}}()
    for v = 1:maximum(c)
        push!(Clusters, Vector{Int64}())
    end
    for i = 1:length(c)
        push!(Clusters[c[i]],i)
    end
    return Clusters
end

"""
Run Step 1 of the Louvain algorithm: iterate through nodes and greedily move
nodes to adjacent clusters.
"""
function LambdaLouvain_Step(A::SparseArrays.SparseMatrixCSC{Float64,Int64},w::Vector{Float64},lam::Float64,randflag::Bool=false,maxits::Int64=Inf,clusterpenalty::Float64=0.0)
    n = size(A,1)

    # This permutes the node labels, to add randomization in the Louvain
    # algorithm so that you don't always traverse the nodes in the same order
    if randflag
        p = Random.randperm(n)
        A = A[p,p]
        undop = sortperm(p)
        w = w[p]
    end

    # All nodes start in singleton clusters
    c = collect(1:n)
    @assert(size(w,1) == n)
    improving = true
    its = 0

    # number of clusters
    K = n

    # Neighs[i] lists nodes that share an edge with node i
    # degvec[i] = number of neighbors of node i
    Neighbs, degvec = ConstructAdj(A,n)

    # Store a vector of n clusters
    Clusters = Vector{Vector{Int64}}()
    for v = 1:n
        push!(Clusters, Vector{Int64}())
        push!(Clusters[v],v)    # singleton clusters--one for each node
    end

    # As long as we can keep improving, continue
    while improving && its < maxits
        its += 1
        # println("Iteration $its")
        improving = false
        nextclus = 1

        # visit each node in turn
        for i = 1:n

            # Cluster index for node i
            Ci_ind = c[i]

            # Get the indices of nodes in i's cluster
            Ci = Clusters[Ci_ind]

            clus_size = length(Ci)

            Ni = Neighbs[i]

            ## First few lines here essentially compute the current contribution
            # this node i has on the objective, based on its current cluster placement

            # Weight of negative mistakes currently at i
            neg_inner = w[i]*(sum(w[Ci]) - w[i])

            # Weight of positive mistakes if we remove i from Ci
            pos_inner = sum(A[i,Ci])

            # Increase in mistakes if we remove i from Ci
            total_inner = pos_inner - lam*neg_inner

            BestImprove = 0
            BestC = Ci_ind

            Ai = A[:,i]

            # Get the neighboring clusters of i:
            # there are the only places we would even consider moving node i to
            NC = unique(c[Ni])

            # Now let's see if it's better to move to a nearby cluster, Cj
            for j = 1:length(NC)
                Cj_ind = NC[j]

                # Check how much it would improve to move i to to cluster j
                if Cj_ind == Ci_ind
                    change = 0
                else

                    Cj = Clusters[Cj_ind]

                    # Moving i from Ci to Cj adds negative mistakes
                    neg_outer = w[i]*(sum(w[Cj]))

                    # Moving i from Ci to Cj decreases positive mistakes

                    pos_outer = mysum(Ai,Cj)

                    total_outer = lam*neg_outer - pos_outer

                    # Change in cluster part of objective
                    if clus_size == 1
                        Δclus = clusterpenalty*(log(K-1)-log(K))
                    else
                        Δclus = 0
                    end

                    # This is the overall change in objective if we move i to Cj
                    change = total_outer + total_inner + Δclus

                end

                # Check if this is currently the best possible greedy move to make
                if change < BestImprove
                    BestImprove = change
                    BestC = Cj_ind
                    improving = true
                end
            end
            # Done checking whether there's a better place to greedily move the cluster to.

            # Move i to the best cluster to move it to, do this
            # only if it strictly improves the objective
            if BestImprove < 0
                ci_old = c[i]
                c[i] = BestC

                # Remove i from its old cluster...
                Clusters[ci_old] = setdiff(Clusters[ci_old],i)

                # ...and add it to its new cluster
                push!(Clusters[BestC],i)
                sort!(Clusters[BestC])

                improving = true # we have a reason to keep iterating!
                if clus_size == 1
                    # we moved a singleton to another cluster
                    K -= 1
                end
            end

        end
    end
    if its == 1
        improved = false
    else
        improved = true
        c, Clusters = renumber(c,Clusters)
    end

    if randflag
        c = c[undop]    # if you previously mixed up node order, re-order them now
    end

    return c, improved

end

function mysum(v::SparseArrays.SparseVector{Float64,Int64},inds::Vector{Int64})
    """
    Less numerically stable, but for the purposes of a greedy clustering
    algorithm works as well as it needs to, and is faster.
    """
    s = 0
    for i in inds
        s += v[i]
    end
    return s
end

"""
Step 2 of the Louvain algorithm:
Collapse a clustering into a new network of supernodes and weighted edges, so
that you can then run the same greedy node-moving on supernodes.
"""
function collapse_clustering(A::SparseArrays.SparseMatrixCSC{Float64,Int64},w::Vector{Float64},c::Vector{Int64})

    n = size(A,1)
    Clusters = Vector{Vector{Int64}}()
    ClusterNeighbs = Vector{Set{Int64}}()
    Neighbs, degvec = ConstructAdj(A,n)

    for v = 1:maximum(c)
        push!(Clusters, Vector{Int64}())
        push!(ClusterNeighbs, Set{Int64}())
    end

    for i = 1:n
        push!(Clusters[c[i]],i)
        # Also keep track of which clusters are adjacent clusters
        for j in Neighbs[i]
            if c[j] != c[i]
                push!(ClusterNeighbs[c[i]],c[j])
            end
        end
    end

    # Number of supernodes to form = number of clusters
    N = round(Int64,maximum(c))

    I = Vector{Int64}()
    J = Vector{Int64}()
    V = Vector{Float64}()
    wnew = zeros(N)
    # Construct a new sparse matrix with new node weights
    for i = 1:N
        Ci = sort(Clusters[i])
        wnew[i] = sum(w[Ci])
        ACi = A[:,Ci]
        for j in ClusterNeighbs[i]
            Cj = Clusters[j]
            Eij = sum(ACi[Cj,:])
            push!(I,i)
            push!(J,j)
            push!(V,Eij)
        end
    end

    Anew = SparseArrays.sparse(I,J,V,N,N)
    Anew = Anew+Anew'

    return Anew, wnew
end

# Sort sizes: Return a set of cluster sizes, arranged in order
function ClusterSizes(c)
    c = renumber(c)

    sizes = Vector{Int64}()
    for i = 1:maximum(c)
        inds = findall(x->x==i,c)
        push!(sizes,length(inds))
    end

    return sort(sizes,rev=true)
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

struct __enum_AclEntryProto_AclEntryScopeProto <: ProtoEnum
    ACCESS::Int32
    DEFAULT::Int32
    __enum_AclEntryProto_AclEntryScopeProto() = new(0,1)
end #struct __enum_AclEntryProto_AclEntryScopeProto
const AclEntryProto_AclEntryScopeProto = __enum_AclEntryProto_AclEntryScopeProto()

struct __enum_AclEntryProto_AclEntryTypeProto <: ProtoEnum
    USER::Int32
    GROUP::Int32
    MASK::Int32
    OTHER::Int32
    __enum_AclEntryProto_AclEntryTypeProto() = new(0,1,2,3)
end #struct __enum_AclEntryProto_AclEntryTypeProto
const AclEntryProto_AclEntryTypeProto = __enum_AclEntryProto_AclEntryTypeProto()

struct __enum_AclEntryProto_FsActionProto <: ProtoEnum
    NONE::Int32
    EXECUTE::Int32
    WRITE::Int32
    WRITE_EXECUTE::Int32
    READ::Int32
    READ_EXECUTE::Int32
    READ_WRITE::Int32
    PERM_ALL::Int32
    __enum_AclEntryProto_FsActionProto() = new(0,1,2,3,4,5,6,7)
end #struct __enum_AclEntryProto_FsActionProto
const AclEntryProto_FsActionProto = __enum_AclEntryProto_FsActionProto()

mutable struct AclEntryProto <: ProtoType
    _type::Int32
    scope::Int32
    permissions::Int32
    name::AbstractString
    AclEntryProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AclEntryProto
const __req_AclEntryProto = Symbol[:_type,:scope,:permissions]
meta(t::Type{AclEntryProto}) = meta(t, __req_AclEntryProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct AclStatusProto <: ProtoType
    owner::AbstractString
    group::AbstractString
    sticky::Bool
    entries::Base.Vector{AclEntryProto}
    permission::FsPermissionProto
    AclStatusProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AclStatusProto
const __req_AclStatusProto = Symbol[:owner,:group,:sticky]
meta(t::Type{AclStatusProto}) = meta(t, __req_AclStatusProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct AclEditLogProto <: ProtoType
    src::AbstractString
    entries::Base.Vector{AclEntryProto}
    AclEditLogProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AclEditLogProto
const __req_AclEditLogProto = Symbol[:src]
meta(t::Type{AclEditLogProto}) = meta(t, __req_AclEditLogProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ModifyAclEntriesRequestProto <: ProtoType
    src::AbstractString
    aclSpec::Base.Vector{AclEntryProto}
    ModifyAclEntriesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ModifyAclEntriesRequestProto
const __req_ModifyAclEntriesRequestProto = Symbol[:src]
meta(t::Type{ModifyAclEntriesRequestProto}) = meta(t, __req_ModifyAclEntriesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ModifyAclEntriesResponseProto <: ProtoType
    ModifyAclEntriesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ModifyAclEntriesResponseProto

mutable struct RemoveAclRequestProto <: ProtoType
    src::AbstractString
    RemoveAclRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveAclRequestProto
const __req_RemoveAclRequestProto = Symbol[:src]
meta(t::Type{RemoveAclRequestProto}) = meta(t, __req_RemoveAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct RemoveAclResponseProto <: ProtoType
    RemoveAclResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveAclResponseProto

mutable struct RemoveAclEntriesRequestProto <: ProtoType
    src::AbstractString
    aclSpec::Base.Vector{AclEntryProto}
    RemoveAclEntriesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveAclEntriesRequestProto
const __req_RemoveAclEntriesRequestProto = Symbol[:src]
meta(t::Type{RemoveAclEntriesRequestProto}) = meta(t, __req_RemoveAclEntriesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct RemoveAclEntriesResponseProto <: ProtoType
    RemoveAclEntriesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveAclEntriesResponseProto

mutable struct RemoveDefaultAclRequestProto <: ProtoType
    src::AbstractString
    RemoveDefaultAclRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveDefaultAclRequestProto
const __req_RemoveDefaultAclRequestProto = Symbol[:src]
meta(t::Type{RemoveDefaultAclRequestProto}) = meta(t, __req_RemoveDefaultAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct RemoveDefaultAclResponseProto <: ProtoType
    RemoveDefaultAclResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RemoveDefaultAclResponseProto

mutable struct SetAclRequestProto <: ProtoType
    src::AbstractString
    aclSpec::Base.Vector{AclEntryProto}
    SetAclRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SetAclRequestProto
const __req_SetAclRequestProto = Symbol[:src]
meta(t::Type{SetAclRequestProto}) = meta(t, __req_SetAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct SetAclResponseProto <: ProtoType
    SetAclResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SetAclResponseProto

mutable struct GetAclStatusRequestProto <: ProtoType
    src::AbstractString
    GetAclStatusRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetAclStatusRequestProto
const __req_GetAclStatusRequestProto = Symbol[:src]
meta(t::Type{GetAclStatusRequestProto}) = meta(t, __req_GetAclStatusRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GetAclStatusResponseProto <: ProtoType
    result::AclStatusProto
    GetAclStatusResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetAclStatusResponseProto
const __req_GetAclStatusResponseProto = Symbol[:result]
meta(t::Type{GetAclStatusResponseProto}) = meta(t, __req_GetAclStatusResponseProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export AclEntryProto_AclEntryScopeProto, AclEntryProto_AclEntryTypeProto, AclEntryProto_FsActionProto, AclEntryProto, AclStatusProto, AclEditLogProto, ModifyAclEntriesRequestProto, ModifyAclEntriesResponseProto, RemoveAclRequestProto, RemoveAclResponseProto, RemoveAclEntriesRequestProto, RemoveAclEntriesResponseProto, RemoveDefaultAclRequestProto, RemoveDefaultAclResponseProto, SetAclRequestProto, SetAclResponseProto, GetAclStatusRequestProto, GetAclStatusResponseProto


text:
module OpenTelemetryExporterOtlpProtoGrpc

export OtlpProtoGrpcTraceExporter

using gRPCClient

import OpenTelemetryAPI
import OpenTelemetrySDK
import OpenTelemetryProto

const API = OpenTelemetryAPI
const SDK = OpenTelemetrySDK
const Otlp = OpenTelemetryProto.OpentelemetryClients
const Trace = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.trace.v1
const Resource = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.resource.v1
const Common = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.common.v1

struct OtlpProtoGrpcTraceExporter{T} <: SDK.AbstractExporter
    client::T
end

function OtlpProtoGrpcTraceExporter(;
    url = "http://localhost:4317",
    is_blocking = true,
    kw...,
)
    if is_blocking
        client = Otlp.TraceServiceBlockingClient(url; kw...)
    else
        client = Otlp.TraceServiceClient(url; kw...)
    end
    OtlpProtoGrpcTraceExporter(client)
end

function SDK.export!(se::OtlpProtoGrpcTraceExporter, sp)
    res, status = Otlp.Export(se.client, convert(Otlp.ExportTraceServiceRequest, sp))
    if gRPCClient.gRPCCheck(status; throw_error = false)
        SDK.EXPORT_SUCCESS
    else
        SDK.EXPORT_FAILURE
    end
end

Base.convert(t::Type{Otlp.ExportTraceServiceRequest}, s::API.Span) = convert(t, [s])

function Base.convert(::Type{Otlp.ExportTraceServiceRequest}, spans::Vector)
    Otlp.ExportTraceServiceRequest(
        resource_spans = [
            Trace.ResourceSpans(
                resource = convert(Resource.Resource, spans[1].tracer.provider.resource),
                schema_url = spans[1].tracer.provider.resource.schema_url,
                # Typically only one element is contained
                # TODO: maybe group spans by resource first?
                instrumentation_library_spans = [
                    convert(Trace.InstrumentationLibrarySpans, spans),
                ],
            ),
        ],
    )
end

function Base.convert(::Type{Resource.Resource}, r::SDK.Resource)
    Resource.Resource(
        attributes = convert(Vector{Common.KeyValue}, r.attributes),
        dropped_attributes_count = API.n_dropped(r.attributes),
    )
end

function Base.convert(::Type{Trace.InstrumentationLibrarySpans}, spans::Vector)
    Trace.InstrumentationLibrarySpans(
        instrumentation_library = convert(
            Common.InstrumentationLibrary,
            spans[1].tracer.instrumentation,
        ),
        schema_url = spans[1].tracer.provider.resource.schema_url,
        spans = [convert(Trace.Span, s) for s in spans],
    )
end

function Base.convert(::Type{Common.InstrumentationLibrary}, info::API.InstrumentationInfo)
    Common.InstrumentationLibrary(name = info.name, version = string(info.version))
end

# TODO: use method call instead of accessing field directly
function Base.convert(::Type{Trace.Span}, s::API.Span)
    Trace.Span(
        trace_id = reinterpret(UInt8, [API.span_context(s).trace_id]),
        span_id = reinterpret(UInt8, [API.span_context(s).span_id]),
        trace_state = string(API.span_context(s).trace_state),
        parent_span_id = if isnothing(s.parent_span_context)
            UInt8[]
        else
            reinterpret(UInt8, [s.parent_span_context.span_id])
        end,
        name = s.name[],
        kind = Int32(s.kind),
        start_time_unix_nano = s.start_time,
        end_time_unix_nano = s.end_time[],
        attributes = convert(Vector{Common.KeyValue}, s.attributes),
        dropped_attributes_count = UInt32(API.n_dropped(s.attributes)),
        events = [convert(Trace.Span_Event, e) for e in s.events],
        dropped_events_count = UInt32(API.n_dropped(s.events)),
        links = [convert(Trace.Span_Link, e) for e in s.links],
        dropped_links_count = UInt32(API.n_dropped(s.links)),
        status = convert(Trace.Status, s.status[]),
    )
end

function Base.convert(
    ::Type{Vector{Common.KeyValue}},
    attrs::Union{API.StaticAttrs,API.DynamicAttrs},
)
    [
        Common.KeyValue(key = k, value = convert(Common.AnyValue, v)) for
        (k, v) in pairs(attrs)
    ]
end

Base.convert(::Type{Common.AnyValue}, v::String) = Common.AnyValue(string_value = v)
Base.convert(::Type{Common.AnyValue}, v::Bool) = Common.AnyValue(bool_value = v)
Base.convert(::Type{Common.AnyValue}, v::Int) = Common.AnyValue(int_value = v)
Base.convert(::Type{Common.AnyValue}, v::Float64) = Common.AnyValue(double_value = v)
Base.convert(::Type{Common.AnyValue}, v::Vector{UInt8}) = Common.AnyValue(bytes_value = v)
Base.convert(::Type{Common.AnyValue}, v::Vector) =
    Common.AnyValue(array_value = convert(Common.ArrayValue, v))
Base.convert(::Type{Common.AnyValue}, v::Union{API.StaticAttrs,API.DynamicAttrs}) =
    Common.AnyValue(kvlist_value = convert(Common.KeyValueList, v))

Base.convert(::Type{Common.ArrayValue}, v::Vector) =
    Common.ArrayValue(values = [convert(Common.AnyValue, x) for x in v])
Base.convert(::Type{Common.KeyValueList}, v::Union{API.StaticAttrs,API.DynamicAttrs}) =
    Common.KeyValueList(values = convert(Vector{Common.KeyValue}, v))

function Base.convert(::Type{Trace.Span_Event}, event::API.Event)
    Trace.Span_Event(
        time_unix_nano = event.timestamp,
        name = event.name,
        attributes = convert(Vector{Common.KeyValue}, event.attributes),
        dropped_attributes_count = API.n_dropped(event.attributes),
    )
end

function Base.convert(::Type{Trace.Span_Link}, link::API.Link)
    Trace.Span_Link(
        trace_id = reinterpret(UInt8, [link.context.trace_id]),
        span_id = reinterpret(UInt8, [link.context.span_id]),
        trace_state = string(link.context.trace_state),
        attributes = convert(Vector{Common.KeyValue}, link.attributes),
        dropped_attributes_count = API.n_dropped(link.attributes),
    )
end

function Base.convert(::Type{Trace.Status}, status::SDK.SpanStatus)
    Trace.Status(code = Int32(status.code), message = something(status.description, ""))
end

end # module


text:
#
# Computes the total potential energy of the system
#

function utotal(p,data)
  utotal = 0.
  for ip in 1:(data.N)-1
    for jp in ip+1:(data.N)
      rij = pbcseparation(p[ip],p[jp],data.side)
      if rij < data.cutoff
        utotal += upair(rij,data)
      end
    end
  end
  return utotal
end


text:
#=
Functions to implement proportional selection.
Returns a new population where the individuals are selected in proportion to their fitness.
=#

export  propsel, propsel!
@doc """ function propsel( pop::Population, fitness::Function )
Apply proportional selection to Population pop using fitness,
and return the result.
"""
function propsel( pop::Population, dfe::Function, fitness_table::Dict{Int64,Float64} )
  new_pop = deepcopy(pop)
  propsel!( new_pop, dfe, fitness_table )
  new_pop
end

@doc """ function propsel( pop::Population, fitness::Function )
Apply proportional selection to Population pop using fitness to generate
and return a new selected population of size n.
"""
function propsel( pop::Population, n::Int64, dfe::Function, fitness_table::Dict{Int64,Float64} )
  new_pop = zeros(Int64,n)
  fmax = 0.0
  for p in pop
    fitp = dfe_fitness( p, dfe, fitness_table )
    if fitp > fmax
      fmax = fitp
    end
  end
  if fmax == 0.0
    # all elements have fitness zero
    return
  end

  N = length(pop)
  k = 0
  while k < n
    i = rand(1:N)
    w = dfe_fitness(pop[i],dfe, fitness_table ) / fmax
    if rand() < w
      new_pop[k + 1] = pop[i]
      k += 1
    end
  end
  new_pop
end

#=
@doc """function propsel!(p::Population, fitness::Vector{Float64} )
Conduct proportional selection in-place.
"""
function propsel!( pop::Population, fitness::Vector{Float64}  )
  fmax = maximum(fitness)
  if fmax == 0
    # all elements have fitness zero
    return
  end

  n = length(pop)
  selected = zeros(Int64, n)
  k = 0
  while k < n
    i = rand(1:n)
    w = fitness[pop[i]] / fmax
    if rand() < w
      selected[k + 1] = i
      k += 1
    end
  end
  pop[:] = [ pop[selected[i]] for i = 1:n ]
end
=#

@doc """function propsel!(p::Population, fitness::Function, dfe::Function )
Conduct proportional selection in-place.
"""
function propsel!( pop::Population, dfe::Function, fitness_table::Dict{Int64,Float64} )
  fmax = 0.0
  for p in pop
    fitp = dfe_fitness( p, dfe, fitness_table )
    if fitp > fmax
      fmax = fitp
    end
  end
  if fmax == 0.0
    # all elements have fitness zero
    return
  end

  n = length(pop)
  selected = zeros(Int64, n)
  k = 0
  while k < n
    i = rand(1:n)
    w = dfe_fitness(pop[i], dfe, fitness_table ) / fmax
    if rand() < w
      selected[k + 1] = i
      k += 1
    end
  end
  pop[:] = [ pop[selected[i]] for i = 1:n ]
end

@doc """ function propsel_alt()
  Roulette wheel version of proportial selection
  Less efficient in almost all circumstances than the above version.
  Written as check of correctness for the above version
"""
function propsel_alt( pop::Population, dfe::Function)
  fitness_table = Dict{Int64,Float64}()
  N = length(pop)
  new_pop = zeros(Int64,N)
  fitdict = Dict{Int64,Float64}()
  popctr = pop_counter( pop )
  for p in pop
    fit = dfe_fitness(p, dfe, fitness_table )  # set fitness of p to be dfe(p).
  end
  sum_fitness = 0.0
  for p in keys(popctr)
    dfit = dfe_fitness( p, dfe, fitness_table )*popctr[p]
    sum_fitness += dfit
  end
  for p in keys(popctr)
    fitdict[p] = dfe_fitness( p, dfe, fitness_table )*popctr[p]/sum_fitness
  end
  for i = 1:N
    r = rand()
    sumfit = 0.0
    for p in keys(fitdict)
      sumfit += fitdict[p]
      if r < sumfit
        new_pop[i] = p
        break
      end
    end
  end
  new_pop
end

@doc """function propsel!(p::Population, fitness::Vector{Float64} )
Conduct proportional selection in-place.
"""
function propsel!( pop::Vector{Chromosome}, fitness::Vector{Float64}; maxfit::Float64=0.0  )
  if maxfit == 0.0
    maxfit = maximum(fitness)
    if maxfit == 0.0
      # all elements have fitness zero
      return
    end
  end
  n = length(pop)
  selected = zeros(Int64, n)
  k = 0
  while k < n
    i = rand(1:n)
    w = pop[i].fitness / maxfit
    if rand() < w
      selected[k + 1] = i
      k += 1
    end
  end
  pop[:] = [ pop[selected[i]] for i = 1:n ]
end


text:
"""
    CenteredBinaryRBM(a, b, w, offset_v = 0, offset_h = 0)

Construct a centered RBM with binary visible and hidden units
"""
function CenteredBinaryRBM(
    a::AbstractArray, b::AbstractArray, w::AbstractArray,
    offset_v::AbstractArray, offset_h::AbstractArray
)
    return CenteredRBM(RBMs.BinaryRBM(a, b, w), offset_v, offset_h)
end

function CenteredBinaryRBM(a::AbstractArray, b::AbstractArray, w::AbstractArray)
    return CenteredRBM(RBMs.BinaryRBM(a, b, w))
end


text:
@testset "643.maximum-average-subarray-i.jl" begin
    @test find_max_average([1, 12, -5, -6, 50, 3], 4) == 51 / 4
    @test find_max_average([13, 20, -1, -6, 28, 30, 10, -4, -10, 22, 26, -21, -1, -9, -9,
                            -27, -3, 21, -13, -4], 7) == 102 / 7
end


text:
@testset "utilities" begin
    @testset "affine/projective" begin
        @test projective([1, 2, 3]) == [1, 1, 2, 3]
        @test affine(projective([1, 2, 3])) == [1, 2, 3]
    end

    @testset "totaldegree" begin
        @PolyImpl.polyvar x y

        F = PolySystem([
            3.0x^2 + 2.0x * y - 1.0 * y^2,
            (-2.5 + 2im)x^2 + -3.0x * y + 4.0 * y^3])

        G, solutions = totaldegree(F)

        @test length(solutions) == 6
        @test degrees(G) == [2, 3]

        for sol in solutions
            @test norm(evaluate(G, sol)) ≈ 0.0 atol = 1e-12
        end

        @test any(x -> abs(1 - x) > 1e-8, norm.(solutions)) == true

        G, solutions = totaldegree(F, unit_roots=true)
        @test all(x -> norm(1 - abs2.(x)) < 1e-8, solutions)

        @test_throws ErrorException totaldegree(PolySystem([x-1, x-2]))
    end

    @testset "randomsystem" begin
        F = randomsystem(Complex64, 2, 3, mindegree=2, maxdegree=5)

        @test nvariables(F) == 3
        @test length(F) == 2
        @test F isa FixedPolySystem.PolySystem{Complex64}

        @test all(d -> 2 ≤ d ≤ 5, degrees(F))

        @test randomsystem(2, 3) isa FixedPolySystem.PolySystem{Complex128}

        G = randomsystem([2, 3, 4], [:x, :y])
        @test nvariables(G) == 2
        @test length(G) == 3
        @test degrees(G) == [2, 3, 4]
    end
end