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

name = "Minuit2"
version = v"6.18.04"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/root-project/root/archive/v6-18-04.tar.gz", "82421a5f0486a2c66170300285dce49a961e3459cb5290f6fa579ef617dc8b0a"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/root-*/math/minuit2

# Apply patch to add `project` command to CMake file
atomic_patch -p1 $WORKSPACE/srcdir/patches/cmake-add-project.patch

sed -i '/^add_library(Minuit2$/a SHARED' src/CMakeLists.txt
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS=-DMATHCORE_STANDALONE=1 ..
make -j ${nproc}
make install

if [[ "${target}" == *-mingw* ]]; then
    # We have to manually build the shared library
    mkdir -p "${libdir}"
    cd "${prefix}/lib"
    ar x libMinuit2.dll.a
    cc -shared -o "${libdir}/libMinuit2.dll" *.o
    rm *.o
fi
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module TimeWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractTimeScales = @jimport org.orekit.time.AbstractTimeScales
const BDTScale = @jimport org.orekit.time.BDTScale
const BiFunction = @jimport java.util.function.BiFunction
const BurstSelector = @jimport org.orekit.time.BurstSelector
const ChronologicalComparator = @jimport org.orekit.time.ChronologicalComparator
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const Comparator = @jimport java.util.Comparator
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const Date = @jimport java.util.Date
const DateComponents = @jimport org.orekit.time.DateComponents
const DateTimeComponents = @jimport org.orekit.time.DateTimeComponents
const DatesSelector = @jimport org.orekit.time.DatesSelector
const EOPHistory = @jimport org.orekit.frames.EOPHistory
const Enum = @jimport java.lang.Enum
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldTimeInterpolable = @jimport org.orekit.time.FieldTimeInterpolable
const FieldTimeShiftable = @jimport org.orekit.time.FieldTimeShiftable
const FieldTimeStamped = @jimport org.orekit.time.FieldTimeStamped
const FixedStepSelector = @jimport org.orekit.time.FixedStepSelector
const Function = @jimport java.util.function.Function
const GLONASSDate = @jimport org.orekit.time.GLONASSDate
const GLONASSScale = @jimport org.orekit.time.GLONASSScale
const GMSTScale = @jimport org.orekit.time.GMSTScale
const GNSSDate = @jimport org.orekit.time.GNSSDate
const GPSScale = @jimport org.orekit.time.GPSScale
const GalileoScale = @jimport org.orekit.time.GalileoScale
const IERSConventions = @jimport org.orekit.utils.IERSConventions
const IRNSSScale = @jimport org.orekit.time.IRNSSScale
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const LazyLoadedEop = @jimport org.orekit.frames.LazyLoadedEop
const LazyLoadedTimeScales = @jimport org.orekit.time.LazyLoadedTimeScales
const List = @jimport java.util.List
const Month = @jimport org.orekit.time.Month
const Object = @jimport java.lang.Object
const OffsetModel = @jimport org.orekit.time.OffsetModel
const Optional = @jimport java.util.Optional
const PreloadedTimeScales = @jimport org.orekit.time.PreloadedTimeScales
const QZSSScale = @jimport org.orekit.time.QZSSScale
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SatelliteSystem = @jimport org.orekit.gnss.SatelliteSystem
const Stream = @jimport java.util.stream.Stream
const TAIScale = @jimport org.orekit.time.TAIScale
const TAIUTCDatFilesLoader = @jimport org.orekit.time.TAIUTCDatFilesLoader
const TCBScale = @jimport org.orekit.time.TCBScale
const TCGScale = @jimport org.orekit.time.TCGScale
const TDBScale = @jimport org.orekit.time.TDBScale
const TTScale = @jimport org.orekit.time.TTScale
const TimeComponents = @jimport org.orekit.time.TimeComponents
const TimeInterpolable = @jimport org.orekit.time.TimeInterpolable
const TimeScalarFunction = @jimport org.orekit.time.TimeScalarFunction
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeScales = @jimport org.orekit.time.TimeScales
const TimeScalesFactory = @jimport org.orekit.time.TimeScalesFactory
const TimeShiftable = @jimport org.orekit.time.TimeShiftable
const TimeStamped = @jimport org.orekit.time.TimeStamped
const TimeVectorFunction = @jimport org.orekit.time.TimeVectorFunction
const TimeZone = @jimport java.util.TimeZone
const ToDoubleFunction = @jimport java.util.function.ToDoubleFunction
const ToIntFunction = @jimport java.util.function.ToIntFunction
const ToLongFunction = @jimport java.util.function.ToLongFunction
const UT1Scale = @jimport org.orekit.time.UT1Scale
const UTCScale = @jimport org.orekit.time.UTCScale
const UTCTAIBulletinAFilesLoader = @jimport org.orekit.time.UTCTAIBulletinAFilesLoader
const UTCTAIHistoryFilesLoader = @jimport org.orekit.time.UTCTAIHistoryFilesLoader
const UTCTAIOffset = @jimport org.orekit.time.UTCTAIOffset
const UTCTAIOffsetsLoader = @jimport org.orekit.time.UTCTAIOffsetsLoader
const UtcTaiOffsetLoader = @jimport org.orekit.time.UtcTaiOffsetLoader

export AbsoluteDate
export AbstractTimeScales
export BDTScale
export BurstSelector
export ChronologicalComparator
export DateComponents
export DateTimeComponents
export DatesSelector
export FieldAbsoluteDate
export FieldTimeInterpolable
export FieldTimeShiftable
export FieldTimeStamped
export FixedStepSelector
export GLONASSDate
export GLONASSScale
export GMSTScale
export GNSSDate
export GPSScale
export GalileoScale
export IRNSSScale
export LazyLoadedTimeScales
export Month
export OffsetModel
export PreloadedTimeScales
export QZSSScale
export TAIScale
export TAIUTCDatFilesLoader
export TCBScale
export TCGScale
export TDBScale
export TTScale
export TimeComponents
export TimeInterpolable
export TimeScalarFunction
export TimeScale
export TimeScales
export TimeScalesFactory
export TimeShiftable
export TimeStamped
export TimeVectorFunction
export UT1Scale
export UTCScale
export UTCTAIBulletinAFilesLoader
export UTCTAIHistoryFilesLoader
export UTCTAIOffset
export UTCTAIOffsetsLoader
export UtcTaiOffsetLoader
export add_default_utctai_offsets_loaders
export add_utctai_offsets_loader
export clear_utctai_offsets_loaders
export compare
export compare_to
export create_besselian_epoch
export create_from_week_components
export create_gps_date
export create_jd_date
export create_julian_epoch
export create_mjd_date
export describe_constable
export duration_from
export equals
export from_seconds
export get_arbitrary_epoch
export get_bdt
export get_beidou_epoch
export get_calendar_week
export get_capitalized_abbreviation
export get_capitalized_name
export get_ccsds_epoch
export get_class
export get_components
export get_date
export get_day
export get_day_number
export get_day_of_week
export get_day_of_year
export get_declaring_class
export get_eop_history
export get_field
export get_fifties_epoch
export get_first_known_leap_second
export get_future_infinity
export get_galileo_epoch
export get_glonass
export get_glonass_epoch
export get_gmst
export get_gps
export get_gps_epoch
export get_gst
export get_hour
export get_interval_number
export get_irnss
export get_irnss_epoch
export get_j2000_day
export get_j2000_epoch
export get_java_epoch
export get_jd0
export get_julian_epoch
export get_last_known_leap_second
export get_leap
export get_lower_case_abbreviation
export get_lower_case_name
export get_milli_in_week
export get_minute
export get_minutes_from_utc
export get_mjd
export get_mjd_ref
export get_modified_julian_epoch
export get_month
export get_month_enum
export get_name
export get_number
export get_offset
export get_offsets
export get_past_infinity
export get_qzss
export get_qzss_epoch
export get_rollover_reference
export get_sec_in_day
export get_second
export get_seconds_in_local_day
export get_seconds_in_utc_day
export get_slope
export get_start
export get_tai
export get_tcb
export get_tcg
export get_tdb
export get_time
export get_time_scales
export get_tt
export get_upper_case_abbreviation
export get_upper_case_name
export get_ut1
export get_utc
export get_utc_scale
export get_validity_start
export get_week_number
export get_year
export hash_code
export inside_leap
export interpolate
export is_after
export is_after_or_equal_to
export is_before
export is_before_or_equal_to
export is_between
export is_between_or_equal_to
export is_close_to
export is_equal_to
export load_data
export load_offsets
export minute_duration
export name
export notify
export notify_all
export of
export offset_from
export offset_from_tai
export offset_to_tai
export ordinal
export parse_ccsds_calendar_segmented_time_code
export parse_ccsds_day_segmented_time_code
export parse_ccsds_unsegmented_time_code
export parse_date
export parse_date_time
export parse_month
export parse_time
export reversed
export select_dates
export set_rollover_reference
export shifted_by
export still_accepts_data
export then_comparing
export then_comparing_double
export then_comparing_int
export then_comparing_long
export time_scales_offset
export to_absolute_date
export to_date
export to_string
export to_string_rfc3339
export value
export value_of
export values
export wait

include("absolute_date.jl")
include("abstract_time_scales.jl")
include("bdt_scale.jl")
include("burst_selector.jl")
include("chronological_comparator.jl")
include("date_components.jl")
include("date_time_components.jl")
include("dates_selector.jl")
include("field_absolute_date.jl")
include("field_time_interpolable.jl")
include("field_time_shiftable.jl")
include("field_time_stamped.jl")
include("fixed_step_selector.jl")
include("galileo_scale.jl")
include("glonass_date.jl")
include("glonass_scale.jl")
include("gmst_scale.jl")
include("gnss_date.jl")
include("gps_scale.jl")
include("irnss_scale.jl")
include("lazy_loaded_time_scales.jl")
include("month.jl")
include("offset_model.jl")
include("preloaded_time_scales.jl")
include("qzss_scale.jl")
include("tai_scale.jl")
include("taiutc_dat_files_loader.jl")
include("tcb_scale.jl")
include("tcg_scale.jl")
include("tdb_scale.jl")
include("time_components.jl")
include("time_interpolable.jl")
include("time_scalar_function.jl")
include("time_scale.jl")
include("time_scales.jl")
include("time_scales_factory.jl")
include("time_shiftable.jl")
include("time_stamped.jl")
include("time_vector_function.jl")
include("tt_scale.jl")
include("u_t1_scale.jl")
include("utc_scale.jl")
include("utc_tai_offset_loader.jl")
include("utctai_bulletin_a_files_loader.jl")
include("utctai_history_files_loader.jl")
include("utctai_offset.jl")
include("utctai_offsets_loader.jl")

end

ext: jl


text:
abstract type ADBackend end
struct ForwardDiffAD{chunk} <: ADBackend end
struct FluxTrackerAD <: ADBackend end

ADBackend() = ADBackend(ADBACKEND[])
ADBackend(T::Symbol) = ADBackend(Val(T))
function ADBackend(::Val{T}) where {T}
    if T === :forward_diff
        return ForwardDiffAD{CHUNKSIZE[]}
    else
        return FluxTrackerAD
    end
end

getchunksize(::Type{ForwardDiffAD{chunk}}) where chunk = Val(chunk)

"""
getADtype(alg)

Finds the autodifferentiation type of the algorithm `alg`.
"""
getADtype() = ADBackend()
getADtype(s::Sampler) = getADtype(typeof(s))
getADtype(s::Type{<:Sampler{TAlg}}) where {TAlg} = getADtype(TAlg)
getADtype(alg::Hamiltonian) = getADtype(typeof(alg))
getADtype(::Type{<:Hamiltonian{AD}}) where {AD} = AD

"""
gradient(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::Union{Nothing, Sampler}=nothing,
)

Computes the gradient of the log joint of `θ` for the model specified by
`(vi, sampler, model)` using whichever automatic differentation tool is currently active.
"""
@generated function gradient(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::TS=nothing,
) where {TS <: Union{Nothing, Sampler}}
    if TS == Nothing
        return quote
            ad_type = getADtype()
            if ad_type <: ForwardDiffAD
                gradient_forward(θ, vi, model, sampler, getchunksize(ad_type))
            else
                gradient_reverse(θ, vi, model, sampler)
            end
        end
    else
        ad_type = getADtype(TS)
        @assert any(T -> ad_type <: T, (ForwardDiffAD, FluxTrackerAD))
        if ad_type <: ForwardDiffAD
            chunk = getchunksize(ad_type)
            return :(gradient_forward(θ, vi, model, sampler, $chunk))
        else ad_type <: FluxTrackerAD
            return :(gradient_reverse(θ, vi, model, sampler))
        end
    end
end

"""
gradient_forward(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    spl::Union{Nothing, Sampler}=nothing,
    chunk_size::Int=CHUNKSIZE[],
)

Computes the gradient of the log joint of `θ` for the model specified by `(vi, spl, model)`
using forwards-mode AD from ForwardDiff.jl.
"""
function gradient_forward(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::Union{Nothing, Sampler}=nothing,
    ::Val{chunk_size}=Val(CHUNKSIZE[]),
) where chunk_size
    # Record old parameters.
    vals_old, logp_old = copy(vi.vals), copy(vi.logp)

    # Define function to compute log joint.
    function f(θ)
        vi[sampler] = θ
        return -runmodel!(model, vi, sampler).logp
    end

    # Set chunk size and do ForwardMode.
    chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
    config = ForwardDiff.GradientConfig(f, θ, chunk)
    ∂l∂θ = ForwardDiff.gradient!(similar(θ), f, θ, config)
    l = vi.logp.value

    # Replace old parameters to ensure this function doesn't mutate `vi`.
    vi.vals, vi.logp = vals_old, logp_old

    # Strip tracking info from θ to avoid mutating it.
    θ .= ForwardDiff.value.(θ)

    return l, ∂l∂θ
end

"""
gradient_reverse(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::Union{Nothing, Sampler}=nothing,
)

Computes the gradient of the log joint of `θ` for the model specified by
`(vi, sampler, model)` using reverse-mode AD from Flux.jl.
"""
function gradient_reverse(
    θ::AbstractVector{<:Real},
    vi::Turing.VarInfo,
    model::Model,
    sampler::Union{Nothing, Sampler}=nothing,
)
    vals_old, logp_old = copy(vi.vals), copy(vi.logp)

    # Specify objective function.
    function f(θ)
        vi[sampler] = θ
        return -runmodel!(model, vi, sampler).logp
    end

    # Compute forward and reverse passes.
    l_tracked, ȳ = Tracker.forward(f, θ)
    l, ∂l∂θ = Tracker.data(l_tracked), Tracker.data(ȳ(1)[1])

    # Remove tracking info from variables in model (because mutable state).
    vi.vals, vi.logp = vals_old, logp_old

    # Strip tracking info from θ to avoid mutating it.
    θ .= Tracker.data.(θ)

    # Return non-tracked gradient value
    return l, ∂l∂θ
end

import Base: <=
<=(a::Tracker.TrackedReal, b::Tracker.TrackedReal) = a.data <= b.data

function verifygrad(grad::AbstractVector{<:Real})
    if any(isnan, grad) || any(isinf, grad)
        @warn("Numerical error has been found in gradients.")
        @warn("grad = $(grad)")
        return false
    else
        return true
    end
end

import StatsFuns: binomlogpdf
binomlogpdf(n::Int, p::Tracker.TrackedReal, x::Int) = Tracker.track(binomlogpdf, n, p, x)
Tracker.@grad function binomlogpdf(n::Int, p::Tracker.TrackedReal, x::Int)
    return binomlogpdf(n, Tracker.data(p), x),
        Δ->(nothing, Δ * (x / p - (n - x) / (1 - p)), nothing)
end

import StatsFuns: poislogpdf
poislogpdf(v::Tracker.TrackedReal, x::Int) = Tracker.track(poislogpdf, v, x)
Tracker.@grad function poislogpdf(v::Tracker.TrackedReal, x::Int)
      return poislogpdf(Tracker.data(v), x),
          Δ->(Δ * (x/v - 1), nothing)
end

function binomlogpdf(n::Int, p::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(p)
    Δ = ForwardDiff.partials(p)
    return FD(binomlogpdf(n, val, x),  Δ * (x / val - (n - x) / (1 - val)))
end

function poislogpdf(v::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(v)
    Δ = ForwardDiff.partials(v)
    return FD(poislogpdf(val, x), Δ * (x/val - 1))
end


text:
using Flase
using Test
using StaticArrays

@testset "Dense Sheeps" begin
    sheep = Flase.DenseSheeps{10}()
    sheep10 = Flase.DenseSheeps{10}(10)
    zero_grid = SizedMatrix{10,10}(zeros(Int,10,10))
    @testset "Defaults" begin
        @test sheep.capacity == 1
        @test sheep.current_sheep[] == 0
        @test sheep.grid == zero_grid
        @test sheep10.current_sheep[] == 10
        @test sheep10.grid != zero_grid
    end # testset

    @testset "Get/Set" begin
        @test Flase._setNSheep!( sheep, 1, 2, 1 ) == 0
        @test isempty(sheep.diffusion_candidates) == true
        @test sheep.grid[1,2] == 1
        @test Flase._setNSheep!( sheep, 1, 2, 2 ) == 1
        @test sheep.diffusion_candidates[1] == (1,2)
    end # testset

    @testset "iteration" begin
        for (n, ind) in sheep10
            @test sheep10.grid[ind] == n
        end #for
        @test length(sheep10) == 10
    end # testset
end # testset


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using PythonOT

DocMeta.setdocmeta!(PythonOT, :DocTestSetup, :(using PythonOT); recursive=true)

makedocs(;
    modules=[PythonOT],
    authors="David Widmann",
    repo="https://github.com/devmotion/PythonOT.jl/blob/{commit}{path}#{line}",
    sitename="PythonOT.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/PythonOT.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md", "api.md"],
    strict=true,
    checkdocs=:exports,
    doctestfilters=[
        r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
        r"(Array{[a-zA-Z0-9]+,\s?1}|Vector{[a-zA-Z0-9]+})",
        r"(Array{[a-zA-Z0-9]+,\s?2}|Matrix{[a-zA-Z0-9]+})",
    ],
)

deploydocs(; repo="github.com/devmotion/PythonOT.jl", push_preview=true, devbranch="main")


text:
using GLAbstraction, GLPlot, Reactive

window = createdisplay(eyeposition=Vec3(4,4,3), w=1000, h=800)

function zdata(i, j, t)
    x 		= float32(i - 0.5)
	z 		= float32(j - 0.5)
	radius 	= sqrt((x * x) + (z * z))

	r = sin(10.0f0 * radius + t)
    Vec1(r + rand(1.0:0.01:1.1))
end
function zcolor(i, j, t)
	x 		= float32(i - 0.5)
	z 		= float32(j - 0.5)
	radius 	= sqrt((x * x) + (z * z))

	r = sin(10.0f0 * radius + t)
    g = cos(10.0f0 * radius + t)
    b = radius
    return Vec4(r,g,b, 1)
end

N         = 128
texdata   = [zdata(i/N, j/N, 15) for i=1:N, j=1:N]
color     = [zcolor(i/N, j/N, 15) for i=1:N, j=1:N] # Example on how to use react to change the color over time

obj = glplot(texdata, :zscale, primitive=CUBE(), color=color, xscale=0.05f0, yscale=0.05f0, xrange=(-4, 4), yrange=(-4, 4))

# You can look at surface.jl to find out how the primitives look like, and create your own.
# Also, it's pretty easy to extend the shader, which you can find under shader/instance_template.vert
# Its also planned, that you can just upload your own functions and uniforms, to further move computations to the shader.

# you can also simply update the texture, even though it's not nicely exposed by the API yet.
zscale = obj[:zscale]
tcolor = obj[:color]

lift(x-> begin
	update!(zscale, [zdata(i/N, k/N, sin(x/10)*15) for i=1:N, k=1:N])
	update!(tcolor, [zcolor(i/N, k/N, (sin(x)+1f0)*4) for i=1:N, k=1:N])
end, every(0.01))

renderloop(window)


text:
using SparseDiffTools
using VertexSafeGraphs
using Graphs
using Graphs: SimpleGraph
using Test

using Random
Random.seed!(123)

test_graphs = Array{VSafeGraph, 1}(undef, 0)

for _ in 1:20
    nv = rand(5:25)
    ne = rand(1: nv^2)
    inner = SimpleGraph(nv)
    graph = VSafeGraph(inner)
    for e in 1:ne
        v1 = rand(1:nv)
        v2 = rand(1:nv)
        while v1 == v2
            v2 = rand(1:nv)
        end
        add_edge!(graph, v1, v2)
    end
    push!(test_graphs, copy(graph))
end

for i in 1:20
    g = test_graphs[i]
    out_colors = SparseDiffTools.color_graph(g,SparseDiffTools.GreedyD1Color())
    for v = 1:nv(g)
        color = out_colors[v]
        for j in inneighbors(g, v)
            if out_colors[j] == color
                 @test false
            end
        end
    end
    @test true
end


text:
"""
    ResultsStruct

This struct gets updated in the GHPGHX.jl module and is used to record the results/outputs.
"""
Base.@kwdef mutable struct ResultsStruct
    # These get modified by the function below
    # Some of these values are required inputs, those that do not have "= value"
    X_Now::Array{Float64, 1}
    FX_Now::Array{Float64, 1}
    N_Bores::Array{Int64, 1}
    N_Bores_Final::Int64 = 0
    FX_Final::Float64 = 0.0
    Length_Boreholes::Float64 = 0.0
    StorageVolume::Float64 = 0.0
    Depth_DST::Float64 = 0.0
    Mdot_Borehole::Float64 = 0.0
    Radius_Field::Float64 = 0.0
    N_Circuits::Int64 = 0
    L_Circuit::Float64 = 0.0
    Mass_Circuit::Float64 = 0.0

    # These get reset each sizing iteration
    Tmax_GHX::Float64 = -9.99e20
    Tmin_GHX::Float64 = 9.99e20
    Power_GHXPump::Float64 = 0.0
    Qfluid_GHXPump::Float64 = 0.0
    Power_WSHP_H::Float64 = 0.0
    Power_WSHP_C::Float64 = 0.0
    LoadMet_WSHP_H::Float64 = 0.0
    LoadMet_WSHP_C::Float64 = 0.0
    Q_Rejected_Total::Float64 = 0.0
    Q_Absorbed_Total::Float64 = 0.0
    Q_GHX_Net::Float64 = 0.0
    Q_GHX_In::Float64 = 0.0
    Q_GHX_Top::Float64 = 0.0
    Q_GHX_Sides::Float64 = 0.0
    Q_GHX_Bottom::Float64 = 0.0
    Q_GHX_Stored::Float64 = 0.0
    Qf_GHX_Stored::Float64 = 0.0

    # Array results
    total_hours::Int64  # Required argument to be passed by InputsStruct
    P_GHXPump_Hourly::Array{Float64, 1} = zeros(total_hours)
    P_WSHPh_Hourly::Array{Float64, 1} = zeros(total_hours)
    P_WSHPc_Hourly::Array{Float64, 1} = zeros(total_hours)
    Qh_Hourly::Array{Float64, 1} = zeros(total_hours)
    Qc_Hourly::Array{Float64, 1} = zeros(total_hours)
    EWT::Array{Float64, 1} = zeros(total_hours)
end

"""
    init_sizing!(r::ResultsStruct, p::InputsStruct, size_iter::Int64)

Performs the initial sizing for the GHPGHX.jl `size_borefield(p)` function.
"""
function init_sizing!(r::ResultsStruct, p::InputsStruct, size_iter::Int64)
    #Changes each sizing iteration
    r.N_Bores[size_iter] = floor(r.X_Now[size_iter] / p.Depth_Bores) + 1
    r.Length_Boreholes = r.X_Now[size_iter] / r.N_Bores[size_iter]
    if p.SpacingType == 1
        r.StorageVolume = r.Length_Boreholes * r.N_Bores[size_iter] * pi * (p.BoreSpacing * 0.525)^2
    else
        r.StorageVolume = p.BoreSpacing^2 * r.Length_Boreholes * r.N_Bores[size_iter]
    end
    r.Depth_DST = 1.5 * r.Length_Boreholes
    r.Mdot_Borehole = p.Mdot_GHXPump / r.N_Bores[size_iter] * p.N_Series
    r.Radius_Field = (r.StorageVolume / (pi * r.Length_Boreholes))^0.5
    r.N_Circuits = r.N_Bores[size_iter] / p.N_Series
    r.L_Circuit = 2.0 * p.N_Series * r.Length_Boreholes + (p.N_Series - 1) / p.N_Series * r.Radius_Field
    r.Mass_Circuit = r.L_Circuit * pi * p.Ri_Pipe^2 * p.Rho_GHXFluid
    nothing
end

"""
    get_ghpghx_results_for_reopt(r::ResultsStruct, p::InputsStruct)

Extract the required results for input to REopt /job: average hourly dispatch performance and other summary performance metrics.
"""
function get_ghpghx_results_for_reopt(r::ResultsStruct, p::InputsStruct)
    results_dict = Dict{Any,Any}()
    results_dict["number_of_boreholes"] = r.N_Bores_Final
    results_dict["length_boreholes_ft"] = round(r.Length_Boreholes * p.METER_TO_FEET, digits=1)
    results_dict["yearly_heating_heatpump_electric_consumption_series_kw"] = zeros(8760)
    results_dict["yearly_cooling_heatpump_electric_consumption_series_kw"] = zeros(8760)
    results_dict["yearly_ghx_pump_electric_consumption_series_kw"] = zeros(8760)
    # Get average electric consumption
    for yr in 1:p.simulation_years
        results_dict["yearly_heating_heatpump_electric_consumption_series_kw"] += round.(r.P_WSHPh_Hourly[(yr-1)*8760+1:yr*8760] / p.simulation_years, digits=3)
        results_dict["yearly_cooling_heatpump_electric_consumption_series_kw"] += round.(r.P_WSHPc_Hourly[(yr-1)*8760+1:yr*8760] / p.simulation_years, digits=3)
        results_dict["yearly_ghx_pump_electric_consumption_series_kw"] += round.(r.P_GHXPump_Hourly[(yr-1)*8760+1:yr*8760] / p.simulation_years, digits=3)
    end
    results_dict["yearly_total_electric_consumption_series_kw"] =
        results_dict["yearly_heating_heatpump_electric_consumption_series_kw"] +
        results_dict["yearly_cooling_heatpump_electric_consumption_series_kw"] +
        results_dict["yearly_ghx_pump_electric_consumption_series_kw"]
    results_dict["yearly_total_electric_consumption_kwh"] = round(sum(results_dict["yearly_total_electric_consumption_series_kw"]), digits=1)
    results_dict["peak_heating_heatpump_thermal_ton"] = round(p.PeakTons_WSHP_H, digits=3)
    results_dict["peak_cooling_heatpump_thermal_ton"] = round(p.PeakTons_WSHP_C, digits=3)
    results_dict["peak_combined_heatpump_thermal_ton"] = round(p.PeakTons_WSHP_GHX, digits=3)
    results_dict["max_eft_f"] = round(maximum(r.EWT) * 1.8 + 32.0)
    results_dict["min_eft_f"] = round(minimum(r.EWT) * 1.8 + 32.0)
    # Calculate average COP for heating and cooling; estimate allocation of pump power to heating and cooling by thermal energy served
    heating_thermal_kwh = sum(p.HeatingThermalLoadKW)
    cooling_thermal_kwh = sum(p.CoolingThermalLoadKW)
    heating_thermal_frac = heating_thermal_kwh / (heating_thermal_kwh + cooling_thermal_kwh)
    cooling_thermal_frac = cooling_thermal_kwh / (heating_thermal_kwh + cooling_thermal_kwh)
    ghx_pump_electric_kwh = sum(results_dict["yearly_ghx_pump_electric_consumption_series_kw"])
    heating_ghx_pump_electric_kwh = heating_thermal_frac * ghx_pump_electric_kwh
    cooling_ghx_pump_electric_kwh = cooling_thermal_frac * ghx_pump_electric_kwh
    heating_heatpump_electric_kwh = sum(results_dict["yearly_heating_heatpump_electric_consumption_series_kw"])
    cooling_heatpump_electric_kwh = sum(results_dict["yearly_cooling_heatpump_electric_consumption_series_kw"])
    results_dict["heating_cop_avg"] = round(heating_thermal_kwh / (heating_heatpump_electric_kwh + heating_ghx_pump_electric_kwh), digits=3)
    results_dict["cooling_cop_avg"] = round(cooling_thermal_kwh / (cooling_heatpump_electric_kwh + cooling_ghx_pump_electric_kwh), digits=3)
    results_dict["solved_eft_error_f"] = round(r.FX_Final * 1.8, digits=3)
    return results_dict
end


text:
using BenchmarkTools
using BoxLeastSquares
using CSV
using ProgressLogging
using PyCall
using Random
using Statistics

ENV["OMP_NUM_THREADS"] = 1
@info "OMP NUM THREADS=$(ENV["OMP_NUM_THREADS"])"

apt = pyimport("astropy.timeseries")

# generate data
rng = Random.MersenneTwister(444219)
params = (
    depth=0.1,
    period=3,
    duration=0.2
)
function gendata(N)
    t = 20 .* rand(rng, N)
    y = @. 1 - params.depth * (mod(t, params.period) < params.duration)
    yerr = 0.005 .* (rand(rng, N) .+ 1)
    y .+= yerr .* randn(rng, N)
    return t, y, yerr, params
end

function periodogram_astropy(t, y, yerr, periods)
    model = apt.BoxLeastSquares(t, y, dy=yerr)
    periodogram = model.power(periods, params.duration)
    return model, periodogram
end

function periodogram_boxleastsquares(t, y, yerr, periods)
    return BLS(t, y, yerr; duration=params.duration, periods)
end

data_sizes = round.(Int, 10 .^ range(2, 4, length=11))
period_sizes = round.(Int, 10 .^ range(2, 4, length=3))

results = []
@progress for N_data in data_sizes
    t, y, yerr, params = gendata(N_data)
    @progress for N_per in period_sizes
        log_periods = range(log(params.period) - 1, log(params.period) + 1, length=N_per)
        periods = exp.(log_periods)

        bench_jl = @benchmark periodogram_boxleastsquares($t, $y, $yerr, $periods)
        bench_py = @benchmark periodogram_astropy($t, $y, $yerr, $periods)

        res = (N_data, N_per, t_astropy=median(bench_py).time/1e9, t_bls=median(bench_jl).time/1e9)
        @info "" res...
        push!(results, res)
    end
end

CSV.write(joinpath(@__DIR__, "benchmark_results.csv"), results)


text:
# Kalman filter code.
# Please see Lopes Tsay 2011 "Particle Filters and Bayesian Inference In Financial Econometrics"

function push_state_forward(G0, G1, Tau, m, V)
    a_t = G0 + G1 * m
    R_t = G1 * V * G1' + Tau
    return a_t, R_t
end

function calc_loglik_t(F0, F1, Σ, a_t, R_t, y_t)
    f_t = F0 + F1 * a_t
    Q_t = F1 * R_t * F1' + Σ
    e_t = (y_t - f_t)
    loglik = logpdf(MvNormal(Symmetric(Q_t)), e_t)[1]
    return e_t, Q_t, loglik
end

function update_state(F1, a_t, R_t, e_t, Q_t)
    # Kalman Gain
    A_t = R_t * F1' / Q_t
    m_t = a_t + A_t * e_t
    V_t = R_t - A_t * Q_t * A_t'
    return m_t, V_t
end

function iterate_kalman_filter_mv(F0, F1, G0, G1, Σ, Tau, m, V, y_t)
    # Push the state forward
    a_t, R_t = push_state_forward(G0, G1, Tau, m, V)
    # Calculate the likelihood
    e_t, Q_t, loglik = calc_loglik_t(F0, F1, Σ, a_t, R_t, y_t)
    # Calculate the Kalman Gain
    m_t, V_t = update_state(F1, a_t, R_t, e_t, Q_t)

    return m_t, V_t, loglik
end

struct KalmanModel{mutate,T<:Any}
    θ::T
end
function KalmanModel(; kwargs...)
    if haskey(kwargs, :mutate)
        mutate = kwargs[:mutate]
    else
        mutate = true
    end
    KalmanModel{mutate,typeof(kwargs.data)}(kwargs.data)
end

function Base.show(io::IO, x::KalmanModel{mutate}) where {mutate}
    print(io, "KalmanModel(θ=", x.θ, ")\n")
    print(io, "mutate=", mutate)
end

function kalman_filter_mv(F0, F1, G0, G1, Σ, Tau, μ0, Tau0, Y)
    m = KalmanModel(
        F0 = F0,
        F1 = F1,
        G0 = G0,
        G1 = G1,
        Σ = Σ,
        Tau = Tau,
        μ0 = μ0,
        Tau0 = Tau0,
        mutate = true,
    )
    loglik, state = m(Y)
    return loglik, state.ms, state.Vs, state.logliks
end

function init_kalman_state(::KalmanModel{true}, T, m, V, loglik)
    ms = repeat(m, outer = (1, T))
    Vs = repeat(V, outer = (1, 1, T))
    logliks = fill(loglik, T)
    return @NT(ms, Vs, logliks)
end

function init_kalman_state(::KalmanModel{false}, T, m, V, loglik)
    return (ms = m, Vs = reshape(V, size(V, 1), size(V, 2), 1))
end

function update_kalman_state!(::KalmanModel{true}, state, t, m, V, loglik_t)
    state.ms[:, t] = m
    state.Vs[:, :, t] = V
    state.logliks[t] = loglik_t
    return state
end
cat3(args...) = cat(args...; dims = Val(3))
function update_kalman_state!(
    ::KalmanModel{false},
    state,
    t,
    m::Array{Float64,2},
    V::Array{Float64,2},
    loglik_t::Float64,
)
    if t > 1
        ms = hcat(state.ms, m)
        Vs = cat3(state.Vs, reshape(V, size(V, 1), size(V, 2), 1))
    else
        ms = m
        Vs = reshape(V, size(V, 1), size(V, 2), 1)
    end
    return (ms = ms, Vs = Vs)
end

function (model::KalmanModel)(Y)
    θ = model.θ
    T = size(Y)[1]
    m = deepcopy(reshape(θ.μ0, :, 1))
    V = deepcopy(θ.Tau0)
    P = inv(V)
    loglik = zero(θ.μ0[1])

    state = init_kalman_state(model, T, m, V, loglik)
    if !(any(ismissing.(Y[1, :])))
        e_t, Q_t, loglik_t = calc_loglik_t(θ.F0(1), θ.F1(1), θ.Σ(1), m, V, @view(Y[1, :]))
        loglik += loglik_t
        m, V = update_state(θ.F1(1), m, V, e_t, Q_t)
        state = update_kalman_state!(model, state, 1, m, V, loglik_t)
    else
        loglik += zero(loglik)
    end

    if (T > 1)
        for t = 2:T
            # Push system forward
            if !any(ismissing.(Y[t, :]))
                m, V, loglik_t = iterate_kalman_filter_mv(
                    θ.F0(t),
                    θ.F1(t),
                    θ.G0(t),
                    θ.G1(t),
                    θ.Σ(t),
                    θ.Tau(t),
                    m,
                    V,
                    @view(Y[t, :])
                )
                loglik += loglik_t
            else
                m, V = push_state_forward(θ.G0(t), θ.G1(t), θ.Tau(t), m, V)
                loglik_t = zero(loglik)
            end
            state = update_kalman_state!(model, state, t, m, V, loglik_t)
        end
    end

    return loglik, state
end

Zygote.@nograd kalman_filter_mv

function kalman_smoother_mv(F0, F1, G0, G1, Σ, Tau, μ0, Tau0, Y)
    m = KalmanModel(
        F0 = F0,
        F1 = F1,
        G0 = G0,
        G1 = G1,
        Σ = Σ,
        Tau = Tau,
        μ0 = μ0,
        Tau0 = Tau0,
        mutate = true,
    )
    smooth(m, Y)
end

function init_smooth_state(::KalmanModel{true}, D, T, ms, Vs)
    ms_smoothed = deepcopy(ms)
    Vs_smoothed = deepcopy(Vs)

    ms_smoothed[:, 1:(end-1)] .= NaN
    Vs_smoothed[:, :, 1:(end-1)] .= NaN

    Covs_smoothed = zeros(eltype(ms_smoothed), D, D, T - 1)
    state = @NT(ms_smoothed, Vs_smoothed, Covs_smoothed)
end

function init_smooth_state(::KalmanModel{false}, D, T, ms, Vs)
    ms_smoothed = ms[:, end:end]
    Vs_smoothed = Vs[:, :, end:end]
    Covs_smoothed = Array{eltype(ms_smoothed)}(undef, D, D, 0)
    state = @NT(ms_smoothed, Vs_smoothed, Covs_smoothed)
end

function get_smooth_mean_cov(::KalmanModel{true}, state, t)
    ms_next = view(state.ms_smoothed, :, t + 1)
    Vs_next = view(state.Vs_smoothed, :, :, t + 1)
    return ms_next, Vs_next
end

function get_smooth_mean_cov(::KalmanModel{false}, state, t)
    ms_next = view(state.ms_smoothed, :, 1)
    Vs_next = view(state.Vs_smoothed, :, :, 1)
    return ms_next, Vs_next
end

function update_smooth_state!(
    ::KalmanModel{true},
    state,
    t,
    m_smoothed,
    V_smoothed,
    Cov_smoothed,
)
    state.ms_smoothed[:, t] = m_smoothed
    state.Vs_smoothed[:, :, t] = V_smoothed
    state.Covs_smoothed[:, :, t] = Cov_smoothed
    return state
end

function update_smooth_state!(
    ::KalmanModel{false},
    state,
    t,
    m_smoothed,
    V_smoothed,
    Cov_smoothed,
)
    ms_smoothed = hcat(m_smoothed, state.ms_smoothed)
    Vs_smoothed = cat3(V_smoothed, state.Vs_smoothed)
    Covs_smoothed = cat3(Cov_smoothed, state.Covs_smoothed)
    state = @NT(ms_smoothed, Vs_smoothed, Covs_smoothed)
end

function smooth(model::KalmanModel, Y)
    θ = model.θ
    T = size(Y)[1]

    loglik, filt_state = model(Y)
    ms = filt_state.ms
    Vs = filt_state.Vs
    D = size(ms, 1)

    state = init_smooth_state(model, D, T, ms, Vs)

    if (T > 1)
        for t = (T-1):-1:1
            ms_next, Vs_next = get_smooth_mean_cov(model, state, t)
            m_smoothed, V_smoothed, Cov_smoothed = backward_smooth_step_mv(
                θ.G0(t + 1),
                θ.G1(t + 1),
                θ.Tau(t + 1),
                view(ms, :, t),
                view(Vs, :, :, t),
                ms_next,
                Vs_next,
            )
            state =
                update_smooth_state!(model, state, t, m_smoothed, V_smoothed, Cov_smoothed)

        end
    end

    return state
end

Zygote.@nograd kalman_smoother_mv

function backward_smooth_step_mv(G0, G1, Tau, m_t, V_t, m_tp1_sm, V_tp1_sm)
    # One step ahead predictive
    a_tp1, R_tp1 = push_state_forward(G0, G1, Tau, m_t, V_t)
    B = (V_t * G1') / R_tp1
    m_t_sm = m_t + B * (m_tp1_sm - a_tp1)
    V_t_sm = V_t + B * (V_tp1_sm - R_tp1) * B'
    Cov_tp1_t = V_tp1_sm * B'

    return m_t_sm, V_t_sm, Cov_tp1_t
end

function can_mutate(::KalmanModel{mutate}) where {mutate}
    return mutate
end

function sample_smooth(m::KalmanModel, rng, ms, Vs)
    T = size(ms)[2]
    V_fill_in = zero(Vs[:, :, end])
    if can_mutate(m)
        Xs_smoothed = deepcopy(ms)
        Xs_smoothed[:, end] = rand(rng, MvNormal(ms[:, end], Symmetric(Vs[:, :, end])))
    else
        Xs_smoothed = rand(rng, MvNormal(ms[:, end], Symmetric(Vs[:, :, end])), 1)
    end

    if (T > 1)
        for t = (T-1):-1:1
            if can_mutate(m)
                Xs_next = view(Xs_smoothed, :, t + 1)
            else
                Xs_next = view(Xs_smoothed, :, 1)
            end
            m_smoothed, V_smoothed, _ = backward_smooth_step_mv(
                m.θ.G0(t + 1),
                m.θ.G1(t + 1),
                m.θ.Tau(t + 1),
                view(ms, :, t),
                view(Vs, :, :, t),
                Xs_next,
                V_fill_in,
            )
            if can_mutate(m)
                Xs_smoothed[:, t] = rand(rng, MvNormal(m_smoothed, Symmetric(V_smoothed)))
            else
                Xs_smoothed = hcat(
                    rand(rng, MvNormal(m_smoothed, Symmetric(V_smoothed)), 1),
                    Xs_smoothed,
                )
            end
        end
    end

    return Xs_smoothed
end

# Old interface
function draw_posterior_path(rng, ms, Vs, G0, G, Σ, Tau, μ0, Tau0)
    m = KalmanModel(G0 = G0, G1 = G, Σ = Σ, Tau = Tau, μ0 = μ0, Tau0 = Tau0, mutate = true)
    sample_smooth(m, rng, ms, Vs)
end


text:
using CSV, DataFrames, Random, TextAnalysis, Languages, Statistics, PyPlot, Flux, BSON

df_all=CSV.read("/Users/ygx/data/songdata.csv")
categorical!(df_all, :artist)
show(by(df_all, :artist, nrow))

artists=["Queen", "The Beatles", "Michael Jackson", "Eminem", "INXS"]
df=df_all[[x in artists for x in df_all[:artist]],:]
df_all=nothing
Random.seed!(1000);
df=df[shuffle(1:size(df, 1)),:]
df[1,:]

docs=Any[]
for i in 1:size(df,1)
    txt=df.text
    txt=replace(df[i,:].text, "\n" => " ")
    txt=replace(df[i,:].text, "'" => "")
    dm=TextAnalysis.DocumentMetadata(Languages.English(), df[i,:].song,"","")
    doc=StringDocument(txt, dm)
    push!(docs, doc)
end

crps=Corpus(docs)
orig_corpus=deepcopy(crps);
prepare!(crps, strip_non_letters | strip_punctuation | strip_case | strip_stopwords | strip_whitespace)

update_lexicon!(crps)
update_inverse_index!(crps)

m_dtm=DocumentTermMatrix(crps)
word_dict=m_dtm.column_indices

tk_idx(s) = haskey(word_dict, s) ? i=word_dict[s] : i=0

function pad_corpus(c, size)
    M=[]
    for doc in 1:length(c)
        tks = tokens(c[doc])
        if length(tks)>=size
            tk_indexes=[tk_idx(w) for w in tks[1:size]]
        end
        if length(tks)<size
            tk_indexes=zeros(Int64,size-length(tks))
            tk_indexes=vcat(tk_indexes, [tk_idx(w) for w in tks])
        end
        doc==1 ? M=tk_indexes' : M=vcat(M, tk_indexes')
    end
    return M
end

num_terms_in_songs=[length(tokens(crps[i])) for i in 1:length(crps)]

doc_pad_size=200
padded_docs = pad_corpus(crps, doc_pad_size)
X = padded_docs'

artist_dict = Dict()
for (n, a) in enumerate(unique(df.artist))
   artist_dict["$a"] = n
end
artist_dict

artist_indexes=[artist_dict[df[:artist][i]] for i in 1:size(df,1)]
y = Flux.onehotbatch(artist_indexes, 1:5)

X_train = X[:, 1:649]
y_train = y[:,1:649]
X_test = X[:, 650:end]
y_test = y[:, 650:end]

train_set = [(X_train, y_train)]

function load_embeddings(embedding_file)
    local LL, indexed_words, index
    indexed_words = Vector{String}()
    LL = Vector{Vector{Float32}}()
    open(embedding_file) do f
        index = 1
        for line in eachline(f)
            xs = split(line)
            word = xs[1]
            push!(indexed_words, word)
            push!(LL, parse.(Float32, xs[2:end]))
            index += 1
        end
    end
    return reduce(hcat, LL), indexed_words
end

embeddings, vocab = load_embeddings("glove.6B.300d.txt")
embed_size, max_features = size(embeddings)

#Function to return the index of the word 's' in the embedding (returns 0 if the word is not found)
function vec_idx(s)
    i=findfirst(x -> x==s, vocab)
    i==nothing ? i=0 : i
end

#Function to return the word vector for string 's'
wvec(s) = embeddings[:, vec_idx(s)]

max_features = 300
vocab_size = 8450

embedding_matrix=Flux.glorot_normal(max_features, vocab_size)

for term in m_dtm.terms
    if vec_idx(term)!=0
       embedding_matrix[:,word_dict[term]+1]=wvec(term)
    end
end

m = Chain(x -> embedding_matrix * Flux.onehotbatch(reshape(x, doc_pad_size*size(x,2)), 0:vocab_size-1),
    x -> reshape(x, max_features, doc_pad_size, trunc(Int64(size(x,2)/doc_pad_size))),
    x -> mean(x, dims=2),
    x -> reshape(x, max_features, :),
    Dense(max_features, 5),
    softmax
)

loss_h=[]
accuracy_train=[]
accuracy_test=[]
accuracy(x, y) = mean(Flux.onecold(x) .== Flux.onecold(y))
loss(x, y) = sum(Flux.crossentropy(m(x), y))
optimizer = Flux.Momentum(0.2)

Momentum(0.2, 0.9, IdDict{Any,Any}())

for epoch in 1:400
    Flux.train!(loss, Flux.params(m), train_set, optimizer)
    l = loss(X_train, y_train).data
    push!(loss_h, l)
    accuracy_trn=accuracy(m(X_train).data, y_train)
    push!(accuracy_train, accuracy_trn)
    accuracy_tst=accuracy(m(X_test).data, y_test)
    push!(accuracy_test, accuracy_tst)
    println("$epoch -> loss= $l accuracy train=$accuracy_trn accuracy test=$accuracy_tst")
end