text:
import IncrementalInference: getSample
import Base: convert, show
# import DistributedFactorGraphs: getManifold
import ApproxManifoldProducts: sample

using .Images

export ScatterAlignPose2, PackedScatterAlignPose2

export overlayScatter, overlayScatterMutate

"""
$TYPEDEF

This is but one incarnation for how radar alignment factor could work, treat it as a starting point.

Notes
- Stanard `cvt` argument is lambda function to convert incoming images to user convention of image axes,
  - **Geography map default** `cvt` flips image rows so that Pose2 +xy-axes corresponds to img[-x,+y]
    - i.e. rows down is "North" and columns across from top left corner is "East".
- Use rescale to resize the incoming images for lower resolution (faster) correlations
- Both images passed to the construct must have the same type some matrix of type `T`.

Example
-------
```julia
arp2 = ScatterAlignPose2(img1, img2, 2) # e.g. 2 meters/pixel
```

See also: [`overlayScanMatcher`](@ref)
"""
Base.@kwdef struct ScatterAlignPose2{ H1 <: Union{<:ManifoldKernelDensity, <:HeatmapGridDensity},
                                      H2 <: Union{<:ManifoldKernelDensity, <:HeatmapGridDensity} } <: IIF.AbstractManifoldMinimize
  """ reference image for scan matching. """
  cloud1::H1
  """ test image to scan match against the reference image. """
  cloud2::H2
  """ Common grid scale for both images -- i.e. units/pixel.
  Constructor uses two arguments `gridlength`*`rescale=1`=`gridscale`.
  Arg 0 < `rescale` ≤ 1 is also used to rescale the images to lower resolution for speed. """
  gridscale::Float64 = 1.0
  """ how many heatmap sampled particles to use for mmd ailgnment """
  sample_count::Int  = 100
  """ bandwidth to use for mmd """
  bw::Float64        = 1.0
end

# replace inner constructor with transform on image
ScatterAlignPose2(im1::AbstractMatrix{T},
                  im2::AbstractMatrix{T},
                  domain::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real}};
                  sample_count::Integer=75,
                  bw::Real=5e-5, # from a sensitivity analysis with marine radar data (50 or 100 samples)
                  rescale::Real=1,
                  N::Integer=1000,
                  cvt = (im)->reverse(Images.imresize(im,trunc.(Int, rescale.*size(im))),dims=1)
                ) where {T} = ScatterAlignPose2(;cloud1=HeatmapGridDensity(cvt(im1),domain,N=N),
                                                cloud2=HeatmapGridDensity(cvt(im2),domain,N=N),
                                                gridscale=float(rescale),
                                                sample_count,
                                                bw  )
#

getManifold(::IIF.InstanceType{<:ScatterAlignPose2}) = getManifold(Pose2Pose2)

Base.@kwdef struct _FastRetract{M_ <: AbstractManifold, T}
  M::M_ = SpecialEuclidean(2)
  pTq::T = ProductRepr(MVector(0,0.0), MMatrix{2,2}(1.0, 0.0, 0.0, 1.0))
  p::ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}} = ProductRepr(SA[0.0;0.0], SMatrix{2,2}(1.0, 0, 0, 1))
end

function (_user::_FastRetract)(pCq::AbstractVector{<:Real})
  retract!(_user.M, _user.pTq, _user.p, hat(_user.M, _user.p, pCq))

  return _user.pTq
end

function getSample( cf::CalcFactor{<:ScatterAlignPose2} )
  #
  pVi,  = sample(cf.factor.cloud1, cf.factor.sample_count)
  pts2, = sample(cf.factor.cloud2, cf.factor.sample_count)

  M = getManifold(Pose2)
  e0 = ProductRepr(SVector(0.0,0.0), SMatrix{2,2}(1.0, 0.0, 0.0, 1.0))

  # precalc SE2 points
  R0 = e0.parts[2]
  # source memory
  qVj_ = map(pt->ProductRepr(SVector(pt...), R0), pts2)
  # destination memory
  _pVj_! = map(pt->ProductRepr(MVector(pt...), MMatrix{2,2}(0.0,0.0,0.0,0.0)), pts2)
  _pVj! = map(x->x.parts[1], _pVj_!)

  function pVj(xyr)
    pTq! = _FastRetract()
    for i in eachindex(qVj_)
      Manifolds.compose!(M, _pVj_![i], pTq!(xyr), qVj_[i])
      # _pVj![i][:] = Manifolds.compose(M, pTq!(xyr), qVj_[i]).parts[1]
    end

    _pVj!
  end

  # # get the current relative transform estimate
  # pTq! = _FastRetract() # not really any faster yet
  # # not efficient, but okay for here
  # # pTq(xyr) = exp(M, e0, hat(M, e0, xyr))
  # # move other points with relative transform
  # pVj(xyr) = map(pt->Manifolds.compose(M, pTq!(xyr), pt).parts[1], qVj_)

  bw = SA[cf.factor.bw;]
  cost(xyr) = mmd(M.manifold[1], pVi, pVj(xyr), length(pVi), length(qVj_), cf._allowThreads; bw)

  # return mmd as residual for minimization
  res = Optim.optimize(cost, [10*randn(2); 0.1*randn()] )

  M, e0, hat(M, e0, res.minimizer)
end

function (cf::CalcFactor{<:ScatterAlignPose2})(Xtup, wPp, wPq)
  #

  # TODO move M to CalcFactor, follow IIF #1415
  M, e0, pXq = Xtup

  # get the current relative transform estimate
  wPq_ = Manifolds.compose(M, wPp, exp(M, e0, pXq))

  #TODO allocalte for vee! see Manifolds #412, fix for AD
  Xc = zeros(3)
  vee!(M, Xc, wPq, log(M, wPq, wPq_))
  return Xc
end

"""
    $SIGNATURES

Overlay the two images from `AlignRadarPose2` with the first (red) fixed and transform the second image (blue) according to `tf`.

Notes:
- `tf` is a Manifolds.jl type `::ProductRepr` (or newer `::ArrayPartition`) to represent a `SpecialEuclidean(2)` manifold point.
- `user_offset` translates both are the same.

See also: [`plotScatterAlign`](@ref)
"""
function overlayScatter(sap::ScatterAlignPose2,
                        trans::AbstractVector{<:Real}=[0;0.0],
                        rot::Real=0.0;
                        user_coords = [trans; rot],
                        score=Ref(0.0),
                        sample_count::Integer=sap.sample_count,
                        showscore::Bool=true,
                        findBest::Bool=true  )
  #

  # # sample points from the scatter field
  tfg = initfg()
  addVariable!(tfg, :x0, Pose2)
  addVariable!(tfg, :x1, Pose2)
  addFactor!(tfg, [:x0;:x1], sap; graphinit=false)
  meas = sampleFactor(tfg, :x0x1f1)[1]

  M, pts1, pts2_ = meas
  pts2 = map(pt->pt.parts[1], pts2_)
  e0 = identity_element(M)
  # R0 = e0.parts[2]

  # not efficient, but okay for here
  pTq(xyr=user_coords) = exp(M, e0, hat(M, e0, xyr))

  best_coords = zeros(3)
  if findBest
    # keep tfg separately so that optim can be more efficient
    # tfg = initfg()
    ev_(xyr) = calcFactorResidualTemporary( sap, (Pose2,Pose2),
                                            meas,
                                            (e0,pTq(xyr));
                                            tfg )[1]
    #
    score[] = ev_(user_coords)
    best = Optim.optimize(ev_, user_coords, BFGS())
    best_coords .= round.(best.minimizer,digits=3)

    if showscore
      @info "overlayScatter score" score[] string(best_coords) best.minimum # user_offset
    end
  end

  # transform points1 to frame of 2 -- take p as coordinate expansion point
  pts2T_u = map(pt->Manifolds.compose(M, inv(M, pTq()), pt).parts[1], pts2_)
  pts2T_b = map(pt->Manifolds.compose(M, inv(M, pTq(best_coords)), pt).parts[1], pts2_)

  @cast __pts1[i,j] := pts1[j][i]
  @cast __pts2[i,j] := pts2[j][i]
  @cast __pts2Tu[i,j] := pts2T_u[j][i]
  @cast __pts2Tb[i,j] := pts2T_b[j][i]

  (;pP1=__pts1, qP2=__pts2, pP2_u=__pts2Tu, pP2_b=__pts2Tb, best_coords, user_coords, score=round(score[],digits=5)) #, user_offset)
end

function overlayScatterMutate(sap_::ScatterAlignPose2;
                              sample_count::Integer = sap_.sample_count,
                              bw::Real = sap_.bw,
                              kwargs... )
  #
  sap = ScatterAlignPose2(sap_.cloud1, sap_.cloud2, sap_.gridscale, sample_count, float(bw))
  Caesar.overlayScatter(sap; kwargs...);
end

function Base.show(io::IO, sap::ScatterAlignPose2{H1,H2}) where {H1,H2}
  printstyled(io, "ScatterAlignPose2{", bold=true, color=:blue)
  println(io)
  printstyled(io, "    H1 = ", color=:magenta)
  println(io, H1)
  printstyled(io, "    H2 = ", color=:magenta)
  println(io, H2)
  printstyled(io, " }", color=:blue, bold=true)
  println(io)
  if H1 <: HeatmapGridDensity
    println(io, "  size(.data):     ", size(sap.cloud1.data))
    println(io, "    min/max:       ", round(minimum(sap.cloud1.data),digits=3), "/", round(maximum(sap.cloud1.data),digits=3))
  end
  if H2 <: HeatmapGridDensity
    println(io, "  size(.data):     ", size(sap.cloud2.data))
    println(io, "    min/max:       ", round(minimum(sap.cloud2.data),digits=3), "/", round(maximum(sap.cloud2.data),digits=3))
  end

  println(io, "  gridscale:       ", sap.gridscale)
  println(io, "  sample_count:    ", sap.sample_count)
  println(io, "  bw:              ", sap.bw)
  nothing
end

Base.show(io::IO, ::MIME"text/plain", sap::ScatterAlignPose2) = show(io, sap)
Base.show(io::IO, ::MIME"application/juno.inline", sap::ScatterAlignPose2) = show(io, sap)

## =========================================================================================
## Factor serialization below
## =========================================================================================

Base.@kwdef struct PackedScatterAlignPose2 <: PackedInferenceType
  cloud1::String # PackedHeatmapGridDensity # change to String
  cloud2::String # PackedHeatmapGridDensity # change to String
  gridscale::Float64 = 1.0
  sample_count::Int = 50
  bw::Float64 = 0.01
end

function convert(::Type{<:PackedScatterAlignPose2}, arp::ScatterAlignPose2)

  function _toDensityJson(dens::ManifoldKernelDensity)
    convert(PackedSamplableBelief,dens)
  end
  function _toDensityJson(dens::HeatmapGridDensity)
    cloud1 = convert(PackedHeatmapGridDensity,dens)
    JSON2.write(cloud1)
  end

  cloud1_ = _toDensityJson(arp.cloud1)
  cloud2_ = _toDensityJson(arp.cloud2)

  PackedScatterAlignPose2(;
    cloud1 = cloud1_,
    cloud2 = cloud2_,
    gridscale = arp.gridscale,
    sample_count = arp.sample_count,
    bw = arp.bw )
end

function convert(::Type{<:ScatterAlignPose2}, parp::PackedScatterAlignPose2)
  # first understand the schema friendly belief type to unpack
  _cloud1 = JSON2.read(parp.cloud1)
  _cloud2 = JSON2.read(parp.cloud2)
  # @show _cloud1
  # @show parp.cloud1
  #  _cloud2[Symbol("_type")]
  PackedT1 = DFG.getTypeFromSerializationModule(_cloud1[Symbol("_type")])
  PackedT2 = DFG.getTypeFromSerializationModule(_cloud2[Symbol("_type")])
  # re-unpack into the local PackedT (marshalling)
  # TODO check if there is perhaps optimization to marshal directly from _cloud instead - maybe `obj(;namedtuple...)`
  # pcloud1 = JSON2.read(parp.cloud1, PackedT1)
  # pcloud2 = JSON2.read(parp.cloud2, PackedT2)

  # outT
  outT1 = convert(SamplableBelief, PackedT1)
  outT2 = convert(SamplableBelief, PackedT2)

  # @info "deserialize ScatterAlignPose2" typeof(_cloud1) PackedT1 outT1

  # convert from packed schema friendly to local performance type
  cloud1 = convert(outT1, parp.cloud1)
  cloud2 = convert(outT2, parp.cloud2)

  # and build the final object
  ScatterAlignPose2(
    cloud1,
    cloud2,
    parp.gridscale,
    parp.sample_count,
    parp.bw )
end

#


text:
@testset "$TEST $G" begin

r1 = G(10,20)
r2 = DG(5,10)
@test typeof(r1) == G
@test typeof(r2) == DG
@test nv(r1) == 10
@test ne(r1) == 20
@test nv(r2) == 5
@test ne(r2) == 10

@test G(10, 20, seed=3) == G(10, 20, seed=3)
@test DG(10, 20, seed=3) == DG(10, 20, seed=3)
@test G(10, 20, seed=3) == erdos_renyi(10, 20, G, seed=3)
@test ne(G(10, 40, seed=3)) == 40
@test ne(DG(10, 80, seed=3)) == 80

er = erdos_renyi(100, 0.5, G)
@test er isa G
@test nv(er) == 100
@test 2300 < ne(er) < 2650        # 50*99/2 == 2475.0
@test is_directed(er) == false

er = erdos_renyi(100, 0.5, DG)
@test er isa DG
@test nv(er) == 100
@test 4600 < ne(er) < 5300                 # 50*99/2 == 4950.0
@test is_directed(er) == true

ws = watts_strogatz(10, 4, 0.2,  G)
@test nv(ws) == 10
@test ne(ws) == 20
@test is_directed(ws) == false

ws = watts_strogatz(10, 4, 0.2, DG)
@test nv(ws) == 10
@test ne(ws) == 20
@test is_directed(ws) == true

ba = barabasi_albert(10, 2,  G)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, 2,  G)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == false

ba = barabasi_albert(10, 4, 2, G)
@test nv(ba) == 10
@test ne(ba) == 12
@test is_directed(ba) == false

ba = barabasi_albert!(CompleteGraph(2, G), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 17
@test is_directed(ba) == false

ba = barabasi_albert!(CompleteGraph(2, G), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 17
@test is_directed(ba) == false

ba = barabasi_albert!(CompleteGraph(4, G), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, DG)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == true

ba = barabasi_albert(10, 2, 2, DG)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == true

ba = barabasi_albert(10, 4, 2, DG)
@test nv(ba) == 10
@test ne(ba) == 12
@test is_directed(ba) == true

ba = barabasi_albert!(CompleteDiGraph(2, DG), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == true

ba = barabasi_albert!(CompleteDiGraph(2, DG), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == true

ba = barabasi_albert!(CompleteDiGraph(4, DG), 10, 2)
@test nv(ba) == 10
@test ne(ba) == 24
@test is_directed(ba) == true

fm = static_fitness_model(20, rand(10), G)
@test nv(fm) == 10
@test ne(fm) == 20
@test is_directed(fm) == false

fm = static_fitness_model(20, rand(10), rand(10), DG)
@test nv(fm) == 10
@test ne(fm) == 20
@test is_directed(fm) == true

sf = static_scale_free(10, 20, 2.0, G)
@test nv(sf) == 10
@test ne(sf) == 20
@test is_directed(sf) == false

sf = static_scale_free(10, 20, 2.0, 2.0, DG)
@test nv(sf) == 10
@test ne(sf) == 20
@test is_directed(sf) == true
@test sf isa DG

rr = random_regular_graph(5, 0, G)
@test nv(rr) == 5
@test ne(rr) == 0
@test is_directed(rr) == false
@test rr isa G

rd = random_regular_digraph(10,0)
@test nv(rd) == 10
@test ne(rd) == 0
@test is_directed(rd)

rr = random_regular_graph(6, 3, G, seed=1)
@test all(degree(rr) .== 3)
@test nv(rr) == 6
@test ne(rr) == 9
@test is_directed(rr) == false

g = random_bipartite_regular_graph(6, 8, 4, 3, G, seed=1)
ks = degree(g)
@test all(ks[1:6] .== 4)
@test all(ks[7:14] .== 3)
@test nv(g) == 14
@test ne(g) == 24
@test is_directed(g) == false
@test g isa G

rr = random_regular_graph(100, 5, G)
@test nv(rr) == 100
@test ne(rr) == 250
@test is_directed(rr) == false
for v in vertices(rr)
    @test degree(rr, v) == 5
end
@test rr isa G

rr = random_configuration_model(10, repeat([2,4] ,5), G, seed=3)
@test nv(rr) == 10
@test ne(rr) == 15
@test is_directed(rr) == false
num2 = 0; num4 = 0
for v in vertices(rr)
    d = degree(rr, v)
    @test  d == 2 || d == 4
    d == 2 ? num2 += 1 : num4 += 1
end
@test num4 == 5
@test num2 == 5
@test rr isa G

ks1 = [1,3,2,1]
ks2 = [1,2,2,1,1]

g = random_bipartite_configuration_model(4, 5, ks1, ks2, G, seed=1)
ks = degree(g)
@test all(ks[1:4] .== ks1)
@test all(ks[5:9] .== ks2)
@test nv(g) == 9
@test ne(g) == 7
@test is_directed(g) == false
@test g isa G

rr = random_configuration_model(100, zeros(Int,100), G)
@test nv(rr) == 100
@test ne(rr) == 0
@test is_directed(rr) == false
@test rr isa G

rr = random_configuration_model(3, [2,2,2], G, check_graphical=true)
@test nv(rr) == 3
@test ne(rr) == 3
@test is_directed(rr) == false
@test rr isa G

rd = random_regular_digraph(100, 4, DG)
@test nv(rd) == 100
@test ne(rd) == 400
@test is_directed(rd)
@test std(out_degree(rd)) == 0

rd = random_regular_digraph(100, 4, DG, dir=:in)
@test nv(rd) == 100
@test ne(rd) == 400
@test is_directed(rd)
@test std(in_degree(rd)) == 0

rr = random_regular_graph(10, 8, G, seed=4)
@test nv(rr) == 10
@test ne(rr) == 40
@test is_directed(rr) == false
for v in vertices(rr)
    @test degree(rr, v) == 8
end

rd = random_regular_digraph(10, 8, DG, dir=:out, seed=4)
@test nv(rd) == 10
@test ne(rd) == 80
@test is_directed(rd)

g = stochastic_block_model(2., 3., [100,100], G)
@test  3.5 < mean(degree(g)) < 6.5
g = stochastic_block_model(3., 4., [100,100,100])
@test  9.5 < mean(degree(g)) < 12.5

end # testset

ext: jl


text:
@testset "environment_models" begin
    include("deterministic_distribution_model.jl")
    include("dynamic_distribution_model.jl")
    include("experience_sample_model.jl")
    include("prioritized_sweeping_sample_model.jl")
    include("time_based_sample_model.jl")
end


text:
"""
    Module CommonUtils

General purpose utilities usable by all element types

"""

module CommonUtils
using LinearAlgebra # for I matrix in geometricFactors
using SparseArrays  # for spdiagm

export meshgrid
export geometric_factors
export connect_mesh
export build_node_maps, build_periodic_boundary_maps, build_periodic_boundary_maps!
export rk45_coeffs

export unzip # convert array of tuples to tuples of arrays
unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
export eye,speye
eye(n) = diagm(ones(n))
speye(n) = spdiagm(0 => ones(n))

# 4th order 5-stage low storage Runge Kutta from Carpenter/Kennedy.
function rk45_coeffs()
    rk4a = [            0.0 ...
    -567301805773.0/1357537059087.0 ...
    -2404267990393.0/2016746695238.0 ...
    -3550918686646.0/2091501179385.0  ...
    -1275806237668.0/842570457699.0];

    rk4b = [ 1432997174477.0/9575080441755.0 ...
    5161836677717.0/13612068292357.0 ...
    1720146321549.0/2090206949498.0  ...
    3134564353537.0/4481467310338.0  ...
    2277821191437.0/14882151754819.0]

    rk4c = [ 0.0  ...
    1432997174477.0/9575080441755.0 ...
    2526269341429.0/6820363962896.0 ...
    2006345519317.0/3224310063776.0 ...
    2802321613138.0/2924317926251.0 ...
    1.0];
    return rk4a,rk4b,rk4c
end

include("./CommonUtils/meshgrid.jl")
include("./CommonUtils/geometric_factors.jl")
include("./CommonUtils/connect_mesh.jl")
include("./CommonUtils/build_node_maps.jl")
include("./CommonUtils/build_periodic_boundary_maps.jl")

end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
module PDFExtract

include("rectangle.jl")
include("pdtext.jl")
include("pddocument.jl")

include("jats.jl")
include("tex.jl")
include("table.jl")
include("pdf2xml.jl")

function download_jar(version::String)
    jarfile = normpath(joinpath(@__DIR__,"../deps/pdfextract-$version.jar"))
    if !isfile(jarfile)
        url = "https://github.com/paperai/pdfextract/releases/download/v$version/pdfextract-$version.jar"
        println("Downloading $url")
        download(url, jarfile)
    end
    jarfile
end

using Base.Filesystem
export mvjats
function mvjats(path::String)
    for file in readdir(path)
        endswith(file,".pdf") || continue
        pdffile = file
        xmlfile = file[1:end-4] * ".xml"
        isfile(joinpath(path,xmlfile)) || continue

        dirpath = joinpath(path, file[1:end-4])
        isdir(dirpath) || mkdir(dirpath)
        mv(joinpath(path,pdffile), joinpath(dirpath,pdffile))
        sleep(0.01)
        mv(joinpath(path,xmlfile), joinpath(dirpath,xmlfile))
        sleep(0.01)
    end
end

end


text:
using ZipFile

using AVRDUDE

ver = AVRDUDE.version

install_path = Pkg.dir("AVRDUDE", "deps", "avrdude", "v"*ver)

!isdir(install_path) && mkpath(install_path)

zip_loc = joinpath(install_path, ver*".zip")

# Grab the repo with all the flash tools
if !isfile(zip_loc)
    download("https://github.com/sjkelly/arduino-flash-tools/archive/v$(ver).zip",
             zip_loc)
end

# Determine the AVRDUDE version we need to use
@windows_only (avr_version = "tools_windows")
@osx_only (avr_version = "tools_darwin")

@linux_only begin
    if contains(Sys.MACHINE, "arm")
        avr_version = "tools_linux_arm"
    elseif WORD_SIZE == 32
        avr_version = "tools_linux_32"
    else
        avr_version = "tools_linux_64"
    end
end

avrdude_ziploc = "arduino-flash-tools-"*ver*"/"*avr_version*"/avrdude/"

r = ZipFile.Reader(zip_loc)
for f in r.files
    if startswith(f.name, avrdude_ziploc)
        zipf = replace(f.name, avrdude_ziploc, "")
        isempty(zipf) && continue
        if endswith(zipf, "/") || endswith(zipf, "\\")
            mkpath(joinpath(install_path, zipf))
        else
            open(joinpath(install_path, zipf), "w") do unzipped
                write(unzipped, readall(f))
            end
        end
    end
end

close(r)


text:
@testset "NonConvex1" begin
    p, q = synthetic_gradient(SynthSphere(50))
    z = zeros(Float64, size(p))
    Z = NonConvex1(z=z)(p, q)
    @test Z[38,38] ≈ -9.185495083459775e-17
    @test Z[1,1] ≈ 0.0
    @test Z[76,76] ≈ 0.00013642077595930305
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export VorNeighborhood

struct VorNeighborhood <: NeighborhoodAlgorithm
    base::Float64
end

function VorNeighborhood()
    return VorNeighborhood(1.1)
end

function optimize_neighborhood!(algo::VorNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where T
end

function neighborhood(algo::VorNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where T
    n = length(index.db)
    k = max(1, ceil(Int, log(algo.base, n)))
    knn = search(index, dist, item, KnnResult(k))
    N = Int32[]

    @inbounds for p in knn
        pobj = index.db[p.objID]
        covered = false
        for nearID in N
            d = convert(Float32, dist(index.db[nearID], pobj))
            if d <= p.dist
                covered = true
                break
            end
        end
        if !covered
            push!(N, p.objID)
        end

    end

    return knn, N
end


text:
using Base.Test

import JSONify
using JSON
using DataFrames

# Complete data
data  = DataFrame(a = ["a", "b", "c", "d"], b = rand(4), c = [1, 2, 3, 4])
x     = JSON.json(data)             # Convert DataFrame to JSON
data2 = DataFrame(JSON.parse(x))    # Parse x from JSON to DataFrame
@test data2 == data
@test eltypes(data2) == eltypes(data)

# Missing data
data[1, :a] = NA
x     = JSON.json(data)
data2 = DataFrame(JSON.parse(x))
@test isequal(data2, data)
@test eltypes(data2) == eltypes(data)


text:
##
using MatrixNetworks
using DelimitedFiles
using CSV
A = MatrixNetworks.readSMAT("cropsim-2010-harvestarea.smat")
xy = readdlm("cropsim-2010-harvestarea.xy")
##
df = CSV.read("cropsim-2010-harvestarea-metadata.csv"; header=false)
##
using Plots
## Show the magnitude of the harvest at each area.
scatter(xy[:,1],xy[:,2],markersize=1, markerstrokewidth=0, marker_z = df.Column4)
## Show the degrees
scatter(xy[:,1],xy[:,2],markersize=1, markerstrokewidth=0, marker_z = sum(A,dims=2))
##
using SparseArrays
i = rand(1:size(A,1))
scatter(xy[:,1],xy[:,2],markersize=1, markerstrokewidth=0)
neighs,vals = findnz(A[:,i])
println(A[:,i].nzval)
scatter!([xy[i,1]], [xy[i,2]], markersize=9)
scatter!(xy[neighs,1],xy[neighs,2], markerstrokewidth=0, markersize=5, marker_z = vals)

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
@testset "performance.jl" begin
    @testset "MSE == prior variance" begin
        imodel = PoissonExpModel(20, 0.025, 1)
        smodel = OUModel(1, 1)
        omodel = ExpGainModel(0, 0)
        filter = FullSF(1, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)
        mse = ComputeMSE(sim; num_timesteps = 100000, timestep = 0.01, burnin = 1000)

        @test 0.9 < mse < 1.1

        filter = DiagSF(1, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)
        mse = ComputeMSE(sim; num_timesteps = 100000, timestep = 0.01, burnin = 1000)

        @test 0.9 < mse < 1.1

        filter = BlockSF(1, 1, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)
        mse = ComputeMSE(sim; num_timesteps = 100000, timestep = 0.01, burnin = 1000)

        @test 0.9 < mse < 1.1
    end

    @testset "Benchmark ComputeMSE of FullSF with dim = 1024" begin
        dim = 1024

        imodel = PoissonExpModel(20, 0.025, dim)
        smodel = OUModel(100, 1)
        omodel = ExpGainModel(20., 0.01)
        filter = FullSF(dim, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)

        println("")
        println("Benchmarking ComputeMSE of FullSF with dim = 1024")
        display(@benchmark ComputeMSE($sim; num_timesteps = 1000, timestep = 0.001))
        println("")
        println("")
    end

    @testset "Benchmark ComputeMSE of BlockSF with 128 8x8 blocks" begin
        dim = 1024
        numblocks = 128
        blocksize = 8

        imodel = PoissonExpModel(20, 0.025, dim)
        smodel = OUModel(100, 1)
        omodel = ExpGainModel(20., 0.01)
        filter = BlockSF(numblocks, blocksize, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)

        println("")
        println("Benchmarking ComputeMSE of BlockSF with 128 8x8 blocks")
        display(@benchmark ComputeMSE($sim; num_timesteps = 1000, timestep = 0.001))
        println("")
        println("")
    end

    @testset "Benchmark ComputeMSE of DiagSF with dim = 1024" begin
        dim = 1024

        imodel = PoissonExpModel(20, 0.025, dim)
        smodel = OUModel(100, 1)
        omodel = ExpGainModel(20., 0.01)
        filter = DiagSF(dim, smodel, omodel)

        sim = Simulation(imodel, smodel, omodel, filter)

        println("")
        println("Benchmarking ComputeMSE of DiagSF with dim = 1024")
        display(@benchmark ComputeMSE($sim; num_timesteps = 1000, timestep = 0.001))
        println("")
        println("")
    end
end#performance.jl


text:
# New do-block syntax introduced in: https://github.com/JuliaLang/julia/pull/23718
function rewrite_do(expr::Expr)
    expr.head == :do || error("expression is not a do-block")
    call, body = expr.args
    Expr(:call, call.args[1], body, call.args[2:end]...)
end

iskwarg(x::Any) = isa(x, Expr) && (x.head === :parameters || x.head === :kw)

"""
    extract_kwargs(expr::Expr) -> Vector{Expr}

Extract the :parameters and :kw value into an array of :kw expressions
we don't evaluate any expressions for values yet though.
"""
function extract_kwargs(expr::Expr)
    kwargs = Expr[]
    for x in expr.args[2:end]
        if iskwarg(x)
            if x.head === :parameters
                for kw in x.args
                    push!(kwargs, kw)
                end
            else
                push!(kwargs, x)
            end
        end
    end
    return kwargs
end

"""
    splitdef(ex::Expr; throw::Bool=true) -> Union{Dict{Symbol,Any}, Nothing}

Split a function definition expression into its various components including:

- `:head`: Expression head of the function definition (`:function`, `:(=)`, `:(->)`)
- `:name`: Name of the function (not present for anonymous functions)
- `:params`: Parametric types defined on constructors
- `:args`: Positional arguments of the function
- `:kwargs`: Keyword arguments of the function
- `:rtype`: Return type of the function
- `:whereparams`: Where parameters
- `:body`: Function body (not present for empty functions)

All components listed may not be present in the returned dictionary with the exception of
`:head` which will always be present.

If the provided expression is not a function then an exception will be raised when
`throw=true`. Use `throw=false` avoid raising an exception and return `nothing` instead.

See also: [`combinedef`](@ref)
"""
function splitdef(ex::Expr; throw::Bool=true)
    def = Dict{Symbol,Any}()
    full_ex = ex  # Keep a reference to the full expression

    function invalid_def(section)
        if throw
            msg = "Function definition contains $section\n$(sprint(Meta.dump, full_ex))"
            Base.throw(ArgumentError(msg))
        else
            nothing
        end
    end

    if !(ex.head === :function || ex.head === :(=) || ex.head === :(->))
        return invalid_def("invalid function head `$(repr(ex.head))`")
    end

    def[:head] = ex.head

    if ex.head === :function && length(ex.args) == 1  # empty function definition
        def[:name] = ex.args[1]
        return def
    elseif length(ex.args) == 2  # Expect signature and body
        def[:body] = ex.args[2]
        ex = ex.args[1]  # Focus on the function signature
    else
        quan = length(ex.args) > 2 ? "too many" : "too few"
        return invalid_def("$quan of expression arguments for `$(repr(def[:head]))`")
    end

    # Where parameters
    if ex isa Expr && ex.head === :where
        def[:whereparams] = Any[]

        while ex isa Expr && ex.head === :where
            append!(def[:whereparams], ex.args[2:end])
            ex = ex.args[1]
        end
    end

    # Return type
    if def[:head] !== :(->) && ex isa Expr && ex.head === :(::)
        def[:rtype] = ex.args[2]
        ex = ex.args[1]
    end

    # Determine if the function is anonymous
    anon = (
        def[:head] === :(->) ||
        def[:head] === :function && !(ex isa Expr && ex.head === :call)
    )

    # Arguments and keywords
    if ex isa Expr && (anon && ex.head === :tuple || !anon && ex.head === :call)
        i = anon ? 1 : 2

        if length(ex.args) >= i
            if ex.args[i] isa Expr && ex.args[i].head === :parameters
                def[:kwargs] = ex.args[i].args

                if length(ex.args) > i
                    def[:args] = ex.args[(i + 1):end]
                end
            else
                def[:args] = ex.args[i:end]
            end
        end
    elseif ex isa Expr && anon && ex.head === :block
        # Note: Short-form anonymous functions (:->) will use a block expression when the
        # arguments are divided by semi-colons but do not use commas:
        #
        # (;) -> ...
        # (x;) -> ...
        # (x;y) -> ...
        # (;x) -> ...  # Note: this is an exception to this rule

        for arg in ex.args
            arg isa LineNumberNode && continue

            if !haskey(def, :args)
                def[:args] = [arg]
            elseif !haskey(def, :kwargs)
                def[:kwargs] = [arg]
            else
                return invalid_def("an invalid block expression as arguments")
            end
        end

        !haskey(def, :kwargs) && (def[:kwargs] = [])

    elseif def[:head] === :(->)
        def[:args] = [ex]
    else
        return invalid_def("invalid or missing arguments")
    end

    # Function name and type parameters
    if !anon
        ex = ex.args[1]

        if ex isa Expr && ex.head === :curly
            def[:params] = ex.args[2:end]
            ex = ex.args[1]
        end

        def[:name] = ex
    end

    return def
end

"""
    combinedef(def::Dict{Symbol,Any}) -> Expr

Create a function definition expression from various components. Typically used to construct
a function using the result of [`splitdef`](@ref).

For more details see the documentation on [`splitdef`](@ref).
"""
function combinedef(def::Dict{Symbol,Any})
    # Determine the name of the function including parameterization
    name = if haskey(def, :params)
        Expr(:curly, def[:name], def[:params]...)
    elseif haskey(def, :name)
        def[:name]
    else
        nothing
    end

    # Empty generic function definitions must not contain additional keys
    empty_extras = (:args, :kwargs, :rtype, :whereparams)
    if !haskey(def, :body) && any(k -> haskey(def, k), empty_extras)
        throw(ArgumentError(string(
            "Function definitions without a body must not contain keys: ",
            join(string.('`', repr.(setdiff(empty_extras, keys(def))), '`'), ", ", ", or "),
        )))
    end

    # Combine args and kwargs
    args = Any[]
    haskey(def, :kwargs) && push!(args, Expr(:parameters, def[:kwargs]...))
    haskey(def, :args) && append!(args, def[:args])

    # Create a partial function signature including the name and arguments
    sig = if name !== nothing
        :($name($(args...)))  # Equivalent to `Expr(:call, name, args...)` but faster
    elseif def[:head] === :(->) && length(args) == 1 && !haskey(def, :kwargs)
        args[1]
    else
        :(($(args...),))  # Equivalent to `Expr(:tuple, args...)` but faster
    end

    # Add the return type
    if haskey(def, :rtype)
        sig = Expr(:(::), sig, def[:rtype])
    end

    # Add any where parameters. Note: Always uses the curly where syntax
    if haskey(def, :whereparams)
        sig = Expr(:where, sig, def[:whereparams]...)
    end

    func = if haskey(def, :body)
        Expr(def[:head], sig, def[:body])
    else
        Expr(def[:head], name)
    end

    return func
end


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    perimeters = Dict{Int, Int}()
    for a ∈ 3:n÷3
        for b ∈ a:n÷2
            k = a * a + b * b
            if √k % 1 == 0
                p = a + b + trunc(Int, √k)
                if haskey(perimeters, p)
                    perimeters[p] += 1
                else
                    push!(perimeters, p => 1)
                end
            end
        end
    end
    return findmax(perimeters)[2]
end

compute(1_000)

@benchmark compute(1_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
@testset "objective types" begin
    @testset "autodiff" begin
        # Should throw, as :wah is not a proper autodiff choice
        @test_throws ErrorException OnceDifferentiable(x->x, rand(10); autodiff=:wah)

        for T in (OnceDifferentiable, TwiceDifferentiable)
            odad1 = T(x->5.0, rand(1); autodiff = :finite)
            odad2 = T(x->5.0, rand(1); autodiff = :forward)
            Optim.gradient!(odad1, rand(1))
            Optim.gradient!(odad2, rand(1))
            #    odad3 = T(x->5., rand(1); autodiff = :reverse)
            @test Optim.gradient(odad1) == [0.0]
            @test Optim.gradient(odad2) == [0.0]
            #    @test odad3.g == [0.0]
        end

        for a in (1.0, 5.0)
            xa = rand(1)
            odad1 = OnceDifferentiable(x->a*x[1], xa; autodiff = :finite)
            odad2 = OnceDifferentiable(x->a*x[1], xa; autodiff = :forward)
        #    odad3 = OnceDifferentiable(x->a*x[1], xa; autodiff = :reverse)
            Optim.gradient!(odad1, xa)
            Optim.gradient!(odad2, xa)
            @test Optim.gradient(odad1) ≈ [a]
            @test Optim.gradient(odad2) == [a]
        #    @test odad3.g == [a]
        end
        for a in (1.0, 5.0)
            xa = rand(1)
            odad1 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :finite)
            odad2 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :forward)
        #    odad3 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :reverse)
            Optim.gradient!(odad1, xa)
            Optim.gradient!(odad2, xa)
         @test Optim.gradient(odad1) ≈ 2.0*a*xa
            @test Optim.gradient(odad2) == 2.0*a*xa
        #    @test odad3.g == 2.0*a*xa
        end
        for dtype in (OnceDifferentiable, TwiceDifferentiable)
            for autodiff in (:finite, :forward)
                differentiable = dtype(x->sum(x), rand(2); autodiff = autodiff)
                Optim.value(differentiable)
                Optim.value!(differentiable, rand(2))
                Optim.value_gradient!(differentiable, rand(2))
                Optim.gradient!(differentiable, rand(2))
                dtype == TwiceDifferentiable && Optim.hessian!(differentiable, rand(2))
            end
        end
    end
    @testset "value/grad" begin
        a = 3.0
        x_seed = rand(1)
        odad1 = OnceDifferentiable(x->a*x[1]^2, x_seed)
        Optim.value_gradient!(odad1, x_seed)
        @test Optim.gradient(odad1) ≈ 2 .* a .* (x_seed)
        @testset "call counters" begin
            @test Optim.f_calls(odad1) == 1
            @test Optim.g_calls(odad1) == 1
            @test Optim.h_calls(odad1) == 0
            Optim.value_gradient!(odad1, x_seed .+ 1.0)
            @test Optim.f_calls(odad1) == 2
            @test Optim.g_calls(odad1) == 2
            @test Optim.h_calls(odad1) == 0
        end
        @test Optim.gradient(odad1) ≈ 2 .* a .* (x_seed .+ 1.0)
    end

end


text:
using Test
using QCMaterial

@testset "hamiltonian.generate_ham_1d_hubbard" begin
    nsite = 2
    n_qubit = 2*nsite
    U = 1.0
    t = -0.01
    μ = 0.0

    ham = generate_ham_1d_hubbard(t, U, nsite, μ)

    ham_ref = FermionOperator()
    up1 = up_index(1)
    down1 = down_index(1)
    up2 = up_index(2)
    down2 = down_index(2)

    #クーロン斥力
    ham_ref += FermionOperator("$(up1)^ $(down1)^ $(up1) $(down1)", -U)
    ham_ref += FermionOperator("$(up2)^ $(down2)^ $(up2) $(down2)", -U)

    #ホッピング積分
    ham_ref += FermionOperator("$(up2)^ $(up1)", t)
    ham_ref += FermionOperator("$(up1)^ $(up2)", t)
    ham_ref += FermionOperator("$(down2)^ $(down1)", t)
    ham_ref += FermionOperator("$(down1)^ $(down2)", t)

    #ケミカルポテンシャルの項
    ham_ref += FermionOperator("$(up1)^  $(up1) ", -μ)
    ham_ref += FermionOperator("$(down1)^ $(down1)", -μ)
    ham_ref += FermionOperator("$(up2)^  $(up2) ", -μ)
    ham_ref += FermionOperator("$(down2)^ $(down2)", -μ)

    @test ham == ham_ref
end


text:
# Bio.Seq
# =======
#
# Module for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Seq

export
    Nucleotide,
    DNANucleotide,
    RNANucleotide,
    DNA_A,
    DNA_C,
    DNA_G,
    DNA_T,
    DNA_M,
    DNA_R,
    DNA_W,
    DNA_S,
    DNA_Y,
    DNA_K,
    DNA_V,
    DNA_H,
    DNA_D,
    DNA_B,
    DNA_N,
    DNA_Gap,
    ACGT,
    RNA_A,
    RNA_C,
    RNA_G,
    RNA_U,
    RNA_M,
    RNA_R,
    RNA_W,
    RNA_S,
    RNA_Y,
    RNA_K,
    RNA_V,
    RNA_H,
    RNA_D,
    RNA_B,
    RNA_N,
    RNA_Gap,
    ACGU,
    iscompatible,
    isambiguous,
    ispurine,
    ispyrimidine,
    Sequence,
    BioSequence,
    DNASequence,
    RNASequence,
    AminoAcidSequence,
    CharSequence,
    NucleotideSequence,
    SeqRecord,
    seqname,
    sequence,
    metadata,
    @dna_str,
    @rna_str,
    @aa_str,
    @char_str,
    @biore_str,
    @prosite_str,
    matched,
    captured,
    alphabet,
    gap,
    complement,
    complement!,
    reverse_complement,
    reverse_complement!,
    mismatches,
    ispalindromic,
    hasambiguity,
    isrepetitive,
    ambiguous_positions,
    gc_content,
    SequenceGenerator,
    randdnaseq,
    randrnaseq,
    randaaseq,
    canonical,
    neighbors,
    eachkmer,
    each,
    Composition,
    composition,
    NucleotideCounts,
    Kmer,
    DNAKmer,
    RNAKmer,
    DNACodon,
    RNACodon,
    translate,
    ncbi_trans_table,
    AminoAcid,
    AA_A,
    AA_R,
    AA_N,
    AA_D,
    AA_C,
    AA_Q,
    AA_E,
    AA_G,
    AA_H,
    AA_I,
    AA_L,
    AA_K,
    AA_M,
    AA_F,
    AA_P,
    AA_S,
    AA_T,
    AA_W,
    AA_Y,
    AA_V,
    AA_O,
    AA_U,
    AA_B,
    AA_J,
    AA_Z,
    AA_X,
    AA_Term,
    AA_Gap,
    FASTAReader,
    FASTAWriter,
    FASTASeqRecord,
    FASTQReader,
    FASTQWriter,
    FASTQSeqRecord,
    TwoBitReader,
    TwoBitWriter,
    Alphabet,
    DNAAlphabet,
    RNAAlphabet,
    NucleotideAlphabets,
    AminoAcidAlphabet,
    CharAlphabet,
    NucleotideAlphabet,
    ExactSearchQuery,
    ApproximateSearchQuery,
    approxsearch,
    approxsearchindex,
    approxrsearch,
    approxrsearchindex,
    ReferenceSequence,
    Demultiplexer,
    demultiplex,
    seqmatrix,
    majorityvote

import Bio
using
    BufferedStreams,
    Iterators,
    IndexableBitVectors,
    Bio.Ragel,
    Bio.StringFields,
    Combinatorics

import ..Ragel: tryread!
export tryread!

"""
    alphabet(typ)

Return an iterator of symbols of `typ`.

`typ` is one of `DNANucleotide`, `RNANucleotide`, or `AminoAcid`.
"""
function alphabet end

"""
    gap(typ)

Return the gap symbol of `typ`.

`typ` is one of `DNANucleotide`, `RNANucleotide`, `AminoAcid`, or `Char`.
"""
function gap end

gap(::Type{Char}) = '-'

include("sequence.jl")
include("nucleotide.jl")
include("aminoacid.jl")
include("alphabet.jl")
include("bitindex.jl")
include("bioseq.jl")
include("hash.jl")
include("randseq.jl")
include("kmer.jl")
include("nmask.jl")
include("refseq.jl")
include("eachkmer.jl")
include("composition.jl")
include("geneticcode.jl")
include("seqrecord.jl")
include("demultiplexer.jl")

# Parsing file types
include("fasta/fasta.jl")
include("fastq/fastq.jl")
include("twobit/twobit.jl")

include("search/exact.jl")
include("search/approx.jl")
include("search/re.jl")

include("deprecated.jl")

end  # module Bio.Seq


text:
export twri_objective

function twri_objective(model_full::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing},
                        options::Options, optionswri::TWRIOptions)
    # Load full geometry for out-of-core geometry containers
    dObs.geometry = Geometry(dObs.geometry)
    source.geometry = Geometry(source.geometry)

    # Limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        model = limit_model_to_receiver_area(source.geometry, dObs.geometry, model, options.buffer_size)
    else
        model = model_full
    end

    # Set up Python model structure
    modelPy = devito_model(model, options)
    dtComp = get_dt(model; dt=options.dt_comp)

    # Extrapolate input data to computational grid
    qIn = time_resample(source.data[1],source.geometry,dtComp)[1]
    obsd = typeof(dObs.data[1]) == SegyIO.SeisCon ? convert(Array{Float32,2}, dObs.data[1][1].data) : dObs.data[1]
    dObserved = time_resample(obsd, dObs.geometry, dtComp)[1]

    isnothing(y) ? Y = nothing : Y = time_resample(y.data[1], y.geometry, dtComp)[1]

    # Set up coordinates
    src_coords = setup_grid(source.geometry, model.n)  # shifts source coordinates by origin
    rec_coords = setup_grid(dObs.geometry, model.n)    # shifts rec coordinates by origin

    ~isempty(options.frequencies) ? freqs = options.frequencies : freqs = nothing
    ~isempty(options.frequencies) ? (wfilt, freqs) =  filter_w(qIn, dtComp, freqs) : wfilt = nothing
    obj, gradm, grady = pycall(ac."wri_func", PyObject,
                               modelPy, src_coords, qIn, rec_coords, dObserved, Y,
                               t_sub=options.subsampling_factor, space_order=options.space_order,
                               grad=optionswri.params, grad_corr=optionswri.grad_corr, eps=optionswri.eps,
                               alpha_op=optionswri.comp_alpha, w_fun=optionswri.weight_fun,
                               freq_list=freqs, wfilt=wfilt)

    if (optionswri.params in [:m, :all])
        gradm = remove_padding(gradm, modelPy.padsizes; true_adjoint=options.sum_padding)
        options.limit_m==true && (gradm = extend_gradient(model_full, model, gradm))
    end
    if ~isnothing(grady)
        grady = time_resample(grady, dtComp, dObs.geometry)
        grady = judiVector(dObs.geometry, grady)
    end

    return obj, gradm, grady
end

function filter_w(qIn, dt, freqs)
    ff = FFTW.fftfreq(length(qIn), 1/dt)
    df = ff[2] - ff[1]
    inds = [findmin(abs.(ff.-f))[2] for f in freqs]
    DFT = joDFT(length(qIn); DDT=Float32)
    R = joRestriction(size(DFT,1), inds; RDT=Complex{Float32}, DDT=Complex{Float32})
    qfilt = DFT'*R'*R*DFT*qIn
    return qfilt, freqs
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# THE SIX LOGICAL FUNCTIONS
#
# Bit shifting and rotation (used by the six SHA-XYZ logical functions:
#
#   NOTE:  The naming of R and S appears backwards here (R is a SHIFT and
#   S is a ROTATION) because the SHA2-256/384/512 description document
#   (see http://csrc.nist.gov/cryptval/shs/sha256-384-512.pdf) uses this
#   same "backwards" definition.

# 32-bit Rotate-right (equivalent to S32 in SHA-256) and rotate-left
rrot(b,x,width) = ((x >> b) | (x << (width - b)))
lrot(b,x,width) = ((x << b) | (x >> (width - b)))

# Shift-right (used in SHA-256, SHA-384, and SHA-512):
R(b,x)   = (x >> b)
# 32-bit Rotate-right (used in SHA-256):
S32(b,x) = rrot(b,x,32)
# 64-bit Rotate-right (used in SHA-384 and SHA-512):
S64(b,x) = rrot(b,x,64)
# 64-bit Rotate-left (used in SHA3)
L64(b,x) = lrot(b,x,64)

# Two of six logical functions used in SHA-256, SHA-384, and SHA-512:
Ch(x,y,z)  = ((x & y) ⊻ (~x & z))
Maj(x,y,z) = ((x & y) ⊻ (x & z) ⊻ (y & z))

# Four of six logical functions used in SHA-256:
Sigma0_256(x) = (S32(2,  UInt32(x)) ⊻ S32(13, UInt32(x)) ⊻ S32(22, UInt32(x)))
Sigma1_256(x) = (S32(6,  UInt32(x)) ⊻ S32(11, UInt32(x)) ⊻ S32(25, UInt32(x)))
sigma0_256(x) = (S32(7,  UInt32(x)) ⊻ S32(18, UInt32(x)) ⊻ R(3 ,   UInt32(x)))
sigma1_256(x) = (S32(17, UInt32(x)) ⊻ S32(19, UInt32(x)) ⊻ R(10,   UInt32(x)))

# Four of six logical functions used in SHA-384 and SHA-512:
Sigma0_512(x) = (S64(28, UInt64(x)) ⊻ S64(34, UInt64(x)) ⊻ S64(39, UInt64(x)))
Sigma1_512(x) = (S64(14, UInt64(x)) ⊻ S64(18, UInt64(x)) ⊻ S64(41, UInt64(x)))
sigma0_512(x) = (S64( 1, UInt64(x)) ⊻ S64( 8, UInt64(x)) ⊻ R( 7,   UInt64(x)))
sigma1_512(x) = (S64(19, UInt64(x)) ⊻ S64(61, UInt64(x)) ⊻ R( 6,   UInt64(x)))

# Let's be able to bswap arrays of these types as well
bswap!(x::Vector{<:Integer}) = map!(bswap, x, x)

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# These methods are used solely to help the auto-arima algorithm
# Kwiatkowski, D.;
# Phillips, P. C. B.;
# Schmidt, P.;
# Shin, Y. (1992)
# "Testing the null hypothesis
#  of stationarity against the
#  alternative of a unit root"

function calc_lag_kpss(lags::Bool, n::Int)
    return lags ? trunc(Int, 4*(n/100)^0.25) : trunc(Int, 12*(n/100)^0.25)
end

function σ²_estimator(l::Int, n::Int, ϵ::Vector{Fl}) where Fl
    return ((1/n)*sum(ϵ[t]^2 for t in 1:n)) + ((2/n)*sum((1-(s/(l + 1))) * (sum(ϵ[t]*ϵ[t - s] for t in s + 1:n)) for s in 1:l))
end

function p_value_weighted_average(η::Fl, p1::Fl, cv1::Fl,
                                  p2::Fl, cv2::Fl) where Fl
    return ((p1*(cv1 - cv2)) - ((cv1 - η)*(p1 - p2)))/(cv1 - cv2)
end

function p_value_from_η(η::Fl, crit_vals::Vector{Fl}) where Fl
    if η >= crit_vals[1]
        p_value = 0.01
    elseif η < crit_vals[1] && η >= crit_vals[2]
        p_value = p_value_weighted_average(η, 0.01, crit_vals[1], 0.025, crit_vals[2])
    elseif η < crit_vals[2] && η >= crit_vals[3]
        p_value = p_value_weighted_average(η, 0.025, crit_vals[2], 0.05, crit_vals[3])
    elseif η < crit_vals[3] && η <= crit_vals[4]
        p_value = p_value_weighted_average(η, 0.05, crit_vals[3], 0.1, crit_vals[4])
    else
        p_value = 0.1
    end
    return p_value
end

function calc_η_level(y::Vector{Fl}, lags::Bool) where Fl
    n = length(y)
    ϵ = y .- (mean(y))
    l = calc_lag_kpss(lags, n)
    η = (1/n^2)*sum(sum((((sum(ϵ[i] for i in 1:t))^2)/σ²_estimator(l, n, ϵ)) for t in 1:n))
    return η
end

function kpss_test(y::Vector{Fl}; lags::Bool = true) where Fl
    η = calc_η_level(y, lags)
    crit_vals = [0.739, 0.574, 0.463, 0.347]
    p_value = p_value_from_η(η, crit_vals)
    return p_value
end


text:
using Plots
using Random
using LaTeXStrings
include("../cmd.jl")
include("../simmd.jl")
include("../extramd.jl")
include("../pextramd.jl")
include("../AutoDiffUtils.jl")
include("../Potentials.jl")

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

m = 1
n = 1

true_x = [0.1]
# true_x[randperm(m)[1:20]] .= 0.0
true_y = [0.1]
# true_y[randperm(m)[1:20]] .= 0.0

A = ones(m, n)

l3 = 3.0

f(x, y) = 0.1 * 4^l3 * sum((x .- true_x)' * A * (y .- true_y))
g(x, y) = -f(x, y)

x0 = [0.4]
y0 = [0.3]

let
    max_iter = 5000
    global scatter_cmd = plot()
    global xHist_cmd = zeros(m, max_iter + 1)
    global yHist_cmd = zeros(n, max_iter + 1)
    global m
    global n
    global A
    global f
    global g
    global x0
    global true_x
    global y0
    global true_y

    α = 1.0

    ψ_x = shannon_entropy(α)
    ψ_y = shannon_entropy(α)

    x = ψ_x[2](x0)
    y = ψ_y[2](y0)

    opt = CMD(x, ψ_x, f, y, ψ_y, g, cg!)

    grad_computations = zeros(Int, max_iter + 1)
    for k = 1 : max_iter
      xHist_cmd[:, k], yHist_cmd[:, k] = iterate_primal(opt)
      # grad_computations[k+1] = grad_computations[k] + step!(opt)
      step!(opt)

      opt.x̄ .= ψ_x[2](max.(iterate_primal(opt)[1], 0.0))
      opt.ȳ .= ψ_y[2](max.(iterate_primal(opt)[2], 0.0))
      # @show grad_computations[k+1] - grad_computations[k]
    end
    xHist_cmd[:, end], yHist_cmd[:, end]= iterate_primal(opt)

    @show grad_computations[end]/max_iter

    scatter!(scatter_cmd, xHist_cmd[1, :], yHist_cmd[1, :], xlims=(0.0, 0.6), ylims=(0.0, 0.6))
end

let
    max_iter = 5000
    global scatter_extramd = plot()
    global xHist_extramd = zeros(m, max_iter + 1)
    global yHist_extramd = zeros(n, max_iter + 1)
    global m
    global n
    global A
    global f
    global g
    global x0
    global true_x
    global y0
    global true_y

    α = 1.0

    ψ_x = positive_projection(α)
    ψ_y = positive_projection(α)

    x = ψ_x[2](x0)
    y = ψ_y[2](y0)

    opt = ExtraMD(x, ψ_x, f, y, ψ_y, g)

    iters = zeros(Int, max_iter + 1)
    for k = 1 : max_iter
      xHist_extramd[:, k], yHist_extramd[:, k] = iterate_primal(opt)
      # iters[k+1] = iters[k] + step!(opt)
      step!(opt)

      opt.x̄ .= ψ_x[2](max.(iterate_primal(opt)[1], 0.0))
      opt.ȳ .= ψ_y[2](max.(iterate_primal(opt)[2], 0.0))
    end
    @show xHist_extramd[:, end], yHist_extramd[:, end]= iterate_primal(opt)

     scatter!(scatter_extramd, xHist_extramd[1, :], yHist_extramd[1, :], xlims=(0.0, 0.6), ylims=(0.0, 0.6))
end

out_plot = scatter(xHist_cmd[1, :], yHist_cmd[1, :], xlims=(0.0, 0.6), ylims=(0.0, 0.6), markercolor=silver, aspect_ratio=:equal, size=(500,500), label=L"\mathrm{CMW}", thickness_scaling=2.00, markerstrokewidth=0.1)
scatter!(out_plot, xHist_extramd[1, :], yHist_extramd[1, :], xlims=(0.0, 0.6), ylims=(0.0, 0.6), markercolor=rust, aspect_ratio=:equal, label=L"\mathrm{Extra gradient}", markerstrokewidth=0.1)

savefig(out_plot, "./out/figures/bilinear_positive_3to$(Int(l3)).pdf")
run(`pdfcrop --margins '5 5 5 5' ./out/figures/bilinear_positive_3to$(Int(l3)).pdf ./out/figures/bilinear_positive_3to$(Int(l3)).pdf`)

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text: eval(load_intervalmatrices_overapproximation())


text:
"""
Tools for generating Julia packages for _educational_ purposes.

# Extended Help

## Overview

Have you ever wanted to write a note-set, or a textbook, or a
series of homework assignments with tools within the Julia
ecosystem? Your instincts are correct! Mediums like
`Documenter.jl` websites, `Pluto.jl` notebooks, and `Markdown.jl`
and PDF formats each have their own strengths, and all can help
you write informative, and educational content.

Wouldn't it be nice if there was a tool which allowed you to
write _one_ file, and generate content within all of the following
mediums automatically? Well, it already exists, and it's called
`Literate.jl`.

Alright, but wouldn't it be nice if there was a tool to help you
create packages with structures that `Literate.jl` can benefit
from? That's where `EducationalPkg` comes in! The following
use cases are a few examples that might benefit from
`EducationalPkg`.

1. Writing a note-set which is compiled with `Documenter.jl`
2. Developing a bunch of homework assignments with `Pluto.jl`
3. Showing several investigations, or experiments that you
want accessible in a variety of mediums

## Usage

```julia
julia> import EducationalPkg
edu} # enter `}` in the REPL to enter `EducationalPkg` mode!
edu} generate FutureNobelPrize
```

"""
module EducationalPkg

import Pkg

using Pluto
using Literate
using Documenter
using DocStringExtensions

include(joinpath("EduREPL", "EduREPL.jl"))

module EduREPLInit
    using ..EduREPL
    function __init__()
        if isdefined(Base, :active_repl)
            EduREPL.runEduREPL()
        end
    end
end

end # module


text:
push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, CorrNoise

makedocs(modules = [CorrNoise],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename = "CorrNoise.jl",
    pages = Any[
        "Manual" => "index.md",
    ])

deploydocs(
    repo = "github.com/ziotom78/CorrNoise.jl.git",
)


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
pyplot()
Random.seed!(11)
using BSON: @save
include("../scripts/paper_settings.jl");
include("../../utils/sindy_turing_utils.jl")
# include("../../utils/sindy_statistics_utils.jl")
# include("../../utils/main_algorithm_utils.jl")

################################################################################
#####
##### Rossler example
#####
include(EXAMPLE_FOLDERNAME*"example_rossler.jl")

# Define the multivariate forcing function
num_ctr = 10;
    U_starts = rand(3, num_ctr) .* tspan[2]
    U_widths = [0.0, 0.0, 0.5];
    amplitudes = [1.0, 1.0, 1000.0]
my_U_func_time2(t) = U_func_time_multivariate(t, u0,
                        U_widths, U_starts,
                        F_dim=[1,2,3],
                        amplitudes=amplitudes)

# Get data
sol = solve_rossler_system(U_func_time=my_U_func_time2)
dat = Array(sol)
plot3d(dat[1,:],dat[2,:],dat[3,:],label="Data")
# plot(ts,dat[1,:]);plot!(ts,dat[2,:]);plot!(ts,dat[3,:])
# plot(ts,dat[3,:])
numerical_grad = numerical_derivative(dat, ts)

## Also get baseline true/ideal cases
# Uncontrolled
dat_raw = Array(solve_rossler_system())
numerical_grad_raw = numerical_derivative(dat_raw, ts)
true_grad_raw = zeros(size(dat))
# for i in 1:size(dat,2)
#     true_grad_raw[:,i] = rossler_system(dat_raw[:,i], Float64.(p), [0])
# end

# True control signal
U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time2(t)
end
plot(ts,dat[1,:]);plot!(ts,dat[2,:]);plot!(ts,dat[3,:])
    plot!(ts,U_true')

###
### Fit a naive model
###
sindy_library = Dict("cross_terms"=>2,"constant"=>nothing);

# Upgrade: Use AIC
# val_list = calc_permutations(5,3)
# val_list = combinations(1:5, 3)
val_list = Iterators.product(1:3,1:3,1:3) # DEBUG
(sindy_unctr,best_criterion,all_criteria,all_models) =
    sindyc_ensemble(dat, numerical_grad, sindy_library, val_list,
                    selection_criterion=my_aicc,
                    sparsification_mode="num_terms",
                    selection_dist=Normal(0.0,20),
                    use_clustering_minimization=true)
print_equations(sindy_unctr)

# Integrate this (poor) model
# Callback aborts if it blows up
condition(u,t,integrator) = any(abs.(u).>1e4)
cb = DiscreteCallback(condition, terminate!)
prob_unctr = ODEProblem(sindy_unctr, u0, tspan, [0], callback=cb)
sindy_dat_unctr = Array(solve(prob_unctr, Tsit5(), saveat=ts));

let d = sindy_dat_unctr, d2 = dat
    plot(d[1,:]);plot!(d[2,:]);plot!(d[3,:])
    plot!(d2[1,:]);plot!(d2[2,:]);plot!(d2[3,:])
end
dat_unctr = sindy_dat_unctr

###
### Calculate the residuals
###
# NOTE: do not use Turing, for speed
sindy_grad1 = sindy_unctr(dat, 0)
residual = numerical_grad .- sindy_grad1
noise_guess = 10*abs(median(residual[3,:]))

scatter(residual[3,:])
    hline!([-noise_guess noise_guess])

ctr_guess = process_residual(residual, noise_guess)

###
### Subsample based on residual
###
accepted_ind = subsample_using_residual(residual,
            noise_guess, min_length=4, shorten_length=0)

num_pts = 1000
start_ind = 1
subsample_ind = accepted_ind[start_ind:num_pts+start_ind-1]
dat_sub = dat[:,subsample_ind]
grad_sub = numerical_grad[:,subsample_ind]

# Plot the subsampled points
i = 3;
    plot(dat[i,:], label="Data; NOT residual")
    scatter!(subsample_ind, dat_sub[i,:], color=:blue, label="Subsampled points")
    plot!(U_true[i,:], label="Control Signal")

# Better SINDy model
# val_list = calc_permutations(7,3)
# val_list = combinations(1:5, 3)
val_list = Iterators.product(1:3,1:3,1:3) # DEBUG
(sindy_sub,best_criterion,all_criteria,all_models,best_index) =
    sindyc_ensemble(dat_sub, grad_sub, sindy_library, val_list,
                    selection_criterion=my_aicc,
                    sparsification_mode="num_terms",
                    selection_dist=Normal(0.0,0.1*noise_guess),
                    use_clustering_minimization=false)
println("Best index is $best_index")
print_equations(sindy_sub)
print("True equations")
print_equations(core_dyn_true)
scatter(sum.(val_list), all_criteria)
    title!("AIC for various sparsities")
    xlabel!("Number of nonzero terms")

# Generate the attractor using this system
condition(u,t,integrator) = any(abs.(u).>1e3)
cb = DiscreteCallback(condition, terminate!)
prob_ctr = ODEProblem(sindy_sub, u0, tspan, [0], callback=cb)
sindy_dat_ctr = solve(prob_ctr, Tsit5(), saveat=ts);
dat_ctr = Array(sindy_dat_ctr)

plot3d(sindy_dat_ctr[1,:], sindy_dat_ctr[2,:], sindy_dat_ctr[3,:])

# Also generate the derivatives
sindy_grad_ctr = sindy_sub(dat, 0)

# And then the guessed control signal

ctr_guess2 = process_residual(numerical_grad.-sindy_grad_ctr, noise_guess)

#####
##### Save Rossler data
#####
this_dat_name = DAT_FOLDERNAME*"TMP_dat_library_of_examples_rossler"

# Uncontrolled data
fname = this_dat_name*"uncontrolled.bson"
@save fname dat_raw numerical_grad_raw

# Original ODE and controller variables
fname = this_dat_name*"ode_vars.bson";
@save fname dat U_true #dat_grad

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_bayes.bson";
@save fname dat_unctr noise_guess residual accepted_ind ctr_guess

# SINDY variables 2: after control
fname = this_dat_name*"ctr_vars_sindy.bson";
@save fname dat_ctr sindy_grad_ctr ctr_guess2

# Bayesian variables 2: after control
# fname = this_dat_name*"ctr_vars_turing.bson";
# # @save fname dat_ctr sample_trajectory_noise_ctr sample_gradients_ctr turing_dat_ctr ctr_guess2
# @save fname dat_ctr


text:
module Evaluation

using Statistics
using Distributions

include("utils.jl")
include("stats.jl")
include("show.jl")
include("alias.jl")
include("crit_diagrams.jl")

end


text:
# import CLIMAParameters
import CLIMAParameters
const CP = CLIMAParameters

path_to_params = joinpath(@__DIR__, "toml", "typed_parameters.toml")

@testset "parameter types from file" begin

    # read parameters needed for tests
    toml_dict_64 = CP.create_toml_dict(
        Float64;
        override_file = path_to_params,
        dict_type = "name",
    )

    toml_dict_32 = CP.create_toml_dict(
        Float32;
        override_file = path_to_params,
        dict_type = "name",
    )

    param_names = [
        "int_array_param",
        "int_param",
        "string_param",
        "string_array_param",
        "untyped_param",
        "ft_array_param",
        "light_speed",
    ]

    param_pairs_64 = CP.get_parameter_values(toml_dict_64, param_names)
    nt = (; param_pairs_64...)
    @test typeof(nt.string_param) == String
    @test eltype(nt.string_array_param) == String
    @test typeof(nt.untyped_param) == String

    @test typeof(nt.int_param) == Int
    @test eltype(nt.int_array_param) == Int

    @test typeof(nt.light_speed) == Float64
    @test eltype(nt.ft_array_param) == Float64

    param_pairs_32 = CP.get_parameter_values(toml_dict_32, param_names)

    nt = (; param_pairs_32...)

    @test typeof(nt.string_param) == String
    @test eltype(nt.string_array_param) == String
    @test typeof(nt.untyped_param) == String

    @test typeof(nt.int_param) == Int
    @test eltype(nt.int_array_param) == Int

    @test typeof(nt.light_speed) == Float32
    @test eltype(nt.ft_array_param) == Float32

    @test_throws ErrorException CP.get_parameter_values(
        toml_dict_32,
        "badtype_param",
    )

end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
using Plots

function step_func(x)
    return x > 0
end

function sigmoid(x)
    return 1 / (1 + exp(-x))
end

function relu(x)
    return max(0, x)
end

function identity_func(x)
    return x
end

function softmax(a)
    c = maximum(a)
    exp_a = exp.(a .- c)  # オーバーフロー対策
    sum_exp_a = sum(exp_a)
    y = exp_a ./ sum_exp_a

    return y
end

function init_network()
    network = Dict(
        "W1" => [0.1 0.3 0.5; 0.2 0.4 0.6],
        "b1" => [0.1 0.2 0.3],
        "W2" => [0.1 0.4; 0.2 0.5; 0.3 0.6],
        "b2" => [0.1 0.2],
        "W3" => [0.1 0.3; 0.2 0.4],
        "b3" => [0.1 0.2]
    )
end

function forward(network, x)
    W1, W2, W3 = network["W1"], network["W2"], network["W3"]
    b1, b2, b3 = network["b1"], network["b2"], network["b3"]

    a1 = x * W1 .+ b1
    z1 = sigmoid.(a1)
    a2 = z1 * W2 .+ b2
    z2 = sigmoid.(a2)
    a3 = z2 * W3 .+ b3
    y = identity_func.(a3)

    return y
end

function main()
    network = init_network()
    x = [1.0 0.5]
    y = forward(network, x)
    println(y)
end

function activationFuncTest()
    x = -5.0:0.1:5.0
    y1 = step_func.(x)
    plot(x, y1)
    y2 = sigmoid.(x)
    plot!(x, y2)

    y3 = relu.(x)
    plot(x, y3)

    a = [0.3 2.9 4.0]
    y4 = softmax(a)
    println(y4)
    println(sum(y4))

end

# activationFuncTest()
main()


text:
function updateνH⁺(L, Σ, v, ϵ)
    m, d = size(L)
    @assert m == length(v)
    H = (L' * inv(Σ) * L + ϵ * I)
    H⁺ = inv(H)
    ν = H⁺ * L' * inv(Σ) * v
    ν, H⁺
end
function updateC(L, Σ, v, C = 0.0)
    m, d = size(L)
    C += 0.5 * dot(v, Σ\v)
    C += m/2*log(2pi) + 0.5*logdet(Σ)
    C
end
function updateνH⁺C(L, Σ, v, ϵ)
    updateνH⁺(L, Σ, v, ϵ)..., updateC(L, Σ, v)
end

function partialbridgeodeνH!(::R3, t, νt, Ht, P, (ν, H⁺, C))
    #m, d = size(L)
    #@assert m == length(v)
    Ht[end] = H = inv(H⁺)
    νt[end] = ν

    b̃(t, y, P) = B(t, P)*y + β(t, P)
    dH⁺(t, y, P) = B(t, P)*y + (B(t, P)*y)' - a(t, P)
#    dH(t, y, P) = - B(t, P)'*y - y * B(t,P) + y*a(t, P)*y'
#    dF(t, y, (H,P)) = -B(t, P)'*y + H*a(t, P)*y  + H*β(t, P)
    dC(t, (F, H), P) = dot(β(t, P), F)  + 0.5*dot(F, a(t, P)*F) - 0.5*tr(H*a(t, P))
#    function dS(t, (ν, H⁺, C), P)
#        H = inv(H⁺)
#        b̃(t, ν, P), dH⁺(t, H⁺, P), dC(t, (H*ν, H), P)
#    end

    for i in length(t)-1:-1:1
        dt = t[i] - t[i+1]
        H⁺ = kernelr3(dH⁺, t[i+1], H⁺, dt, P)
        #H = kernelr3(dH, t[i+1], H, dt, P)
        #F = kernelr3(dF, t[i+1], F, dt, (H,P))
        #ν, H⁺, C = kernelr3dot(dS, t[i+1], (ν, H⁺, C), dt, P)

        F = H*ν
        C += dC(t[i+1], (F, H), P)*dt

        νt[i] = ν = kernelr3(b̃, t[i+1], ν, dt, P)
        Ht[i] = H = inv(H⁺)

    end

    νt, Ht, C
end

function updateFHC(L, Σ, v, F, H, ϵ = 0.0, C = 0.0)
    m, d = size(L)
    H += (L' * inv(Σ) * L + ϵ * I)
    F += L' * inv(Σ) * v
    F, H, updateC(L, Σ, v, C)
end

function partialbridgeodeHνH!(::R3, t, Ft, Ht, P, (F, H, C))
    Ht[end] = H
    Ft[end] = F

    dH(t, y, P) = - B(t, P)'*y - y * B(t,P) + y*a(t, P)*y'
    dF(t, y, (H,P)) = -B(t, P)'*y + H*a(t, P)*y  + H*β(t, P)

    for i in length(t)-1:-1:1
        dt = t[i] - t[i+1]
        C += β(t[i+1], P)'*F*dt + 0.5*F'*a(t[i+1], P)*F*dt - 0.5*tr(H*a(t[i+1], P))*dt
        H = kernelr3(dH, t[i+1], H, dt, P)
        F = kernelr3(dF, t[i+1], F, dt, (H, P))
        Ft[i] = F
        Ht[i] = H
    end

    Ft, Ht, C
end

struct Lyap
end
function partialbridgeodeνH!(::Lyap, t, νt, Ht, P, (νend, Hend⁺, C))
    Ht[end] = H = inv(Hend⁺)
    νt[end] = νend
    H⁺ = Hend⁺
    ν = νend
    b̃(t, y, P) = B(t, P)*y + β(t,P)

    for i in length(t)-1:-1:1
        dt = t[i] - t[i+1]
        ν = kernelr3(b̃, t[i+1], ν, dt, P)
        H⁺ = lyapunovpsdbackward_step(t[i+1], H⁺, -dt, P)
        F = Ht[i+1]*νt[i+1]
        C += β(t[i+1], P)'*F*dt  + 0.5*F'*a(t[i+1], P)*F*dt - 0.5*tr(H*a(t[i+1], P))*dt
        νt[i] = ν
        Ht[i] = H = inv(H⁺)
    end
    ν, H⁺, C
end

"""
    PartialBridgeνH

Guided proposal process for diffusion bridge using backward recursion.

    PartialBridgeνH(tt, P, Pt, L, v,ϵ Σ)

    Guided proposal process for a partial diffusion bridge of `P` to `v` on
    the time grid `tt` using guiding term derived from linear process `Pt`.

    PartialBridgeνH(tt, P, Pt, ν, Hend⁺)

    Guided proposal process on the time grid `tt` using guiding term derived from
    linear process `Pt` with backwards equation initialized at `ν, Hend⁺`.

"""
struct PartialBridgeνH{T,TP,TPt,Tν,TH,TC} <: ContinuousTimeProcess{T}
    Target::TP
    Pt::TPt
    tt::Vector{Float64}
    ν::Vector{Tν}
    H::Vector{TH}
    C::TC
    PartialBridgeνH(P::TP, Pt::TPt, tt, νt::Vector{Tν}, Ht::Vector{TH}, C::TC=0.0) where {TP,TPt,Tν,TH,TC} =
        new{Bridge.valtype(P),TP,TPt,Tν,TH,TC}(P, Pt, tt, νt, Ht, C)

    # 6-7 arg
    function PartialBridgeνH(tt_, P, Pt, L, v::Tv, ϵ, Σ = Bridge.outer(zero(v))) where {Tv}
        tt = collect(tt_)
        N = length(tt)
        m, d = size(L)
        TH = typeof(SMatrix{d,d}(1.0I))
        Ht = zeros(TH, N)
        Tν = typeof(@SVector zeros(d))
        νt = zeros(Tν, N)
        ν, H⁺, C = updateνH⁺C(L, Σ, v, ϵ)
        _, _, C = partialbridgeodeνH!(R3(), tt, νt, Ht, Pt, (ν, H⁺, C))
        PartialBridgeνH(P, Pt, tt, νt, Ht, C)
    end
end
# 5 arg
function partialbridgeνH(tt_, P, Pt, νend::Tv, Hend⁺::TH) where {Tv,TH}
    tt = collect(tt_)
    N = length(tt)
    Ht = zeros(TH, N)
    νt = zeros(Tv, N)
    ν, H⁺, C = partialbridgeodeνH!(Lyap(), tt, νt, Ht, Pt, (νend, Hend⁺, C))
    PartialBridgeνH(P, Pt, tt, νt, Ht, C), ν, H⁺, C
end

function _b((i,t)::IndexedTime, x, P::PartialBridgeνH)
    b(t, x, P.Target) + a(t, x, P.Target)*(P.H[i]*(P.ν[i] - x))
end

r((i,t)::IndexedTime, x, P::PartialBridgeνH) = P.H[i]*(P.ν[i] - x)
H((i,t)::IndexedTime, x, P::PartialBridgeνH) = P.H[i]

σ(t, x, P::PartialBridgeνH) = σ(t, x, P.Target)
a(t, x, P::PartialBridgeνH) = a(t, x, P.Target)
Γ(t, x, P::PartialBridgeνH) = Γ(t, x, P.Target)
constdiff(P::PartialBridgeνH) = constdiff(P.Target) && constdiff(P.Pt)

lptilde(x, P::PartialBridgeνH) = -0.5*((P.ν[1] - x)'*P.H[1] * (P.ν - x)) - P.C

function llikelihood(::LeftRule, Xcirc::SamplePath, Po::PartialBridgeνH; skip = 0)
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    for i in 1:length(tt)-1-skip #skip last value, summing over n-1 elements
        s = tt[i]
        x = xx[i]
        r = Bridge.r((i,s), x, Po)

        som += dot( _b((i,s), x, target(Po)) - _b((i,s), x, auxiliary(Po)), r ) * (tt[i+1]-tt[i])
        if !constdiff(Po)
            H = H((i,s), x, Po)
            som -= 0.5*tr( (a((i,s), x, target(Po)) - atilde((i,s), x, Po))*(H) ) * (tt[i+1]-tt[i])
            som += 0.5*( r'*(a((i,s), x, target(Po)) - atilde((i,s), x, Po))*r ) * (tt[i+1]-tt[i])
        end
    end
    som
end


text:
using Parameters
using ModelParameters
using Statistics

abstract type EnsembleNode <: NetworkNode end

Base.@kwdef mutable struct BaseEnsemble{N<:NetworkNode, P, R<:Real} <: EnsembleNode
    @network_node

    instances::Array{N} = []
    weights::Array{P} = [
        Param(1.0; bounds=[0.0,1.0])
    ]

    # Default to weighted sum
    comb_method::Function = (X, weights) -> sum([x * w for (x, w) in zip(X, weights)])

    outflow::Array{R} = []
end

function param_info(node::BaseEnsemble; kwargs...)::Tuple
    values = Float64[w.val for w in node.weights]
    bounds = [w.bounds for w in node.weights]
    param_names = ["weights_n$(i)" for i in 1:length(bounds)]

    return param_names, values, bounds
end

function BaseEnsemble(nodes::Array{<:NetworkNode}; weights::Array{Float64},
                      bounds::Union{Nothing, Array{Tuple}}=nothing,
                      comb_method::Union{Nothing, Function}=nothing)::BaseEnsemble
    if isnothing(bounds)
        num_nodes = length(nodes)
        @assert(length(weights) == num_nodes, "Number of nodes do not match provided number of weights")

        bounds = repeat([(0.0, 1.0)], num_nodes, 1)
    else
        @assert(length(bounds) == length(weights), "Number of bounds do not match provided number of weights")
    end

    p_weights = [Param(w; bounds=b) for (w, b) in zip(weights, bounds)]

    n1 = nodes[1]
    if !isnothing(comb_method)
        tmp = BaseEnsemble{NetworkNode, Param, Float64}(; name=n1.name, area=n1.area, instances=nodes, weights=p_weights, comb_method)
    else
        tmp = BaseEnsemble{NetworkNode, Param, Float64}(; name=n1.name, area=n1.area, instances=nodes, weights=p_weights)
    end

    return tmp
end

function run_node!(ensemble::BaseEnsemble, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing)
    for inst in ensemble.instances
        run_node!(inst, climate; inflow=inflow, extraction=extraction, exchange=exchange)
    end

    ensemble.outflow = ensemble.comb_method([inst.outflow for inst in ensemble.instances], ensemble.weights)
end

function run_timestep!(node, rain, et, ts; inflow=0.0, extraction=0.0, exchange=0.0)
    for inst in node.instances
        run_timestep!(inst, rain, et, ts; inflow=inflow, extraction=extraction, exchange=exchange)
    end

    Qt = node.comb_method([inst.outflow[ts] for inst in node.instances], node.weights)

    append!(node.outflow, Qt)
end

function reset!(ensemble::BaseEnsemble)
    for inst in ensemble.instances
        reset!(inst)
    end

    ensemble.outflow = []
end

# Functions to support calibration of weights
function update_params!(node::BaseEnsemble, weights...)
    for (i, w) in enumerate(weights)
        node.weights[i] = Param(w, bounds=node.weights[i].bounds)
    end
end

# Calibrate each individual instance?
# function calibrate!(ensemble::BaseEnsemble, climate::Climate)
#     ensemble.weights
# end


text:
# lower triangle mat to vec
function mat_to_vec_scale!(vec::Vector{Float64}, mat::AbstractMatrix{Float64})
    k = 1
    m = size(mat, 1)
    for i in 1:m # fill vector from lower triangle, adding off-diagonals twice
        for j in 1:(i - 1)
            vec[k] = 2.0 * mat[i, j]
            k += 1
        end
        vec[k] = mat[i, i]
        k += 1
    end
    return vec
end

# vec to lower triangle mat
function vec_to_mat_noscale!(mat::AbstractMatrix{Float64}, vec::Vector{Float64})
    k = 1
    m = size(mat, 1)
    for i in 1:m, j in 1:i # fill lower triangle from vector
        mat[i, j] = vec[k]
        k += 1
    end
    return mat
end

# # get (i,j) index in symmetric matrix from lower triangle index
# vec_to_mat_idx(k::Int) = div(1 + isqrt(8 * k - 7), 2)

# get lower triangle vector index from (i,j) index in symmetric matrix
function mat_to_vec_idx(i::Int, j::Int)
    return (i < j) ? (div((j - 1) * j, 2) + i) : (div((i - 1) * i, 2) + j)
end


text:
function create_log_node(com)
    com.input[:logs] || return
    back_idx = length(com.backtrack_vec)
    parent_idx = 0
    if length(com.backtrack_vec) > 0
        status = com.backtrack_vec[back_idx].status
        vidx = com.backtrack_vec[back_idx].vidx
        lb = com.backtrack_vec[back_idx].lb
        ub = com.backtrack_vec[back_idx].ub
        best_bound = com.backtrack_vec[back_idx].best_bound

        parent_idx = com.backtrack_vec[back_idx].parent_idx
    else # for initial solve
        status = :Closed
        vidx = 0
        lb = 0
        ub = 0
        best_bound = com.best_bound
    end

    tree_log_node = TreeLogNode(
        back_idx,
        status,
        true, # feasible
        best_bound,
        -1, # initial step nr
        vidx,
        lb,
        ub,
        Dict{Int,Vector{Int}}(), # var_states
        Dict{Int,Vector{Tuple{Symbol,Int,Int,Int}}}(), # var_changes
        Dict{Int,Float64}(), # activity
        Vector{TreeLogNode{typeof(best_bound)}}(), # children
    )

    push!(com.logs, tree_log_node)
    @assert length(com.logs) == back_idx

    if parent_idx > 0
        push!(com.logs[parent_idx].children, tree_log_node)
    end

end

function update_log_node!(com, back_idx; feasible = nothing)
    com.input[:logs] || return
    tree_log_node = com.logs[back_idx]

    # never set an infeasible one to feasible
    # + infeasible nodes are closed
    if feasible !== nothing
        tree_log_node.feasible = feasible
        close_node!(com, back_idx)
    end

    tree_log_node.status = com.backtrack_vec[back_idx].status
    tree_log_node.vidx = com.backtrack_vec[back_idx].vidx
    tree_log_node.lb = com.backtrack_vec[back_idx].lb
    tree_log_node.ub = com.backtrack_vec[back_idx].ub
    tree_log_node.best_bound = com.backtrack_vec[back_idx].best_bound
    tree_log_node.step_nr = com.backtrack_vec[back_idx].step_nr

    step_nr = tree_log_node.step_nr

    variables = com.search_space
    if tree_log_node.status == :Closed
        for var in variables
            # easier to visualize if we save it every step even if unchanged
            # increases time and file size but that is not critical if keep_logs is `true` anyway.
            # Hopefully :D
            tree_log_node.var_states[var.idx] = sort!(values(var))
            if step_nr > 0 && num_changes(var, step_nr) > 0
                tree_log_node.var_changes[var.idx] = copy(view_changes(var, step_nr))
            end
            tree_log_node.activity[var.idx] = var.activity
        end
    end
end

function bfs_list(start_node::CS.TreeLogNode)
    to_process = Vector{CS.TreeLogNode}()
    depths = Vector{Int}()
    nodes_list = Vector{CS.TreeLogNode}()
    num_children = Vector{Int}()

    push!(to_process, start_node)
    push!(depths, 0)
    node = deepcopy(start_node)
    node.children = TreeLogNode[]

    while !isempty(to_process)
        current_node = popfirst!(to_process)
        depth = popfirst!(depths)

        node = deepcopy(current_node)
        push!(num_children, length(node.children))
        node.children = TreeLogNode[]
        push!(nodes_list, node)

        for child_node in current_node.children
            push!(to_process, child_node)
            push!(depths, depth + 1)
        end
    end
    return nodes_list, num_children
end

"""
    sanity_check_logs(log)

Check that
- there are at least 2 closed nodes
- no step nr > 0 is used more than once
"""
function sanity_check_log(log)
    nodes_list, num_children = bfs_list(log)
    nclosed = 0
    step_nrs = Int[]
    for i in 1:length(nodes_list)
        node = nodes_list[i]
        num_child = num_children[i]
        if node.status == :Closed
            nclosed += 1
        end
        if node.step_nr > 0
            push!(step_nrs, node.step_nr)
        end
    end
    passed = nclosed >= 2 && allunique(step_nrs)
    !passed && write("FAILED.json", JSON.json(log))
    return passed
end

function same_logs(log1, log2)
    nodes_list1, num_children1 = bfs_list(log1)
    nodes_list2, num_children2 = bfs_list(log2)
    if length(nodes_list1) != length(nodes_list2)
        println("Different length")
        return false
    end
    for i in 1:length(nodes_list1)
        node1 = nodes_list1[i]
        node2 = nodes_list2[i]
        num_child1 = num_children1[i]
        num_child2 = num_children2[i]
        if num_child1 != num_child2
            println("Different number of children at: ", i)
            return false
        end

        if node1.id != node2.id ||
           node1.status != node2.status ||
           node1.best_bound != node2.best_bound ||
           node1.step_nr != node2.step_nr ||
           node1.vidx != node2.vidx ||
           node1.lb != node2.lb ||
           node1.ub != node2.ub ||
           node1.var_states != node2.var_states ||
           node1.feasible != node2.feasible
            println("Not identical at i=", i)
            println("node1: ")
            println(
                "id: ",
                node1.id,
                " status: ",
                node1.status,
                " best_bound: ",
                node1.best_bound,
                " step_nr: ",
                node1.step_nr,
            )
            println(
                "vidx: ",
                node1.vidx,
                " lb:",
                node1.lb,
                " ub: ",
                node1.ub,
                " var_states: ",
                node1.var_states,
            )
            println("node2: ")
            println(
                "id: ",
                node2.id,
                " status: ",
                node2.status,
                " best_bound: ",
                node2.best_bound,
                " step_nr: ",
                node2.step_nr,
            )
            println(
                "vidx: ",
                node2.vidx,
                " lb:",
                node2.lb,
                " ub: ",
                node2.ub,
                " var_states: ",
                node2.var_states,
            )
            return false
        end
    end
    return true
end

function get_logs(com::CS.CoM)
    logs = Dict{Symbol,Any}()
    nvars = length(com.search_space)

    logs[:init] = Vector{Vector{Int}}(undef, nvars)
    for var in com.init_search_space
        logs[:init][var.idx] = values(var)
    end

    logs[:tree] = com.logs[1]
    return logs
end

function add_variable_mapping(log, args...)
    var_mapping = Dict{Symbol,Any}()
    for var in args
        # transpose such that JSON has it row by row instead of column by column
        var_mapping[var.first] = var_idx.(var.second)'
    end
    log[:variable_mapping] = var_mapping
end

"""
    save_logs(com::CS.CoM, filepath)

Save the tree structure and some additional problem information in a json file `filepath`.
Can be only used if `keep_logs` is set to `true` in the [`solve!`](@ref) call.
"""
function save_logs(com::CS.CoM, filepath, vars...)
    logs = get_logs(com)
    add_variable_mapping(logs, vars...)
    write(filepath, JSON.json(logs))
end


text:
using SPHShockTube1D
# @static if VERSION < v"0.7.0-DEV.2005"
#     using Base.Test
# else
#     using Test
# end

# write your own tests here
(x, ρ, u, v, p) = SPHShockTube1D.runsimulation()
gr()
plot(x, ρ)


text:
using InfoZIP
using Base.Test

using Compat.write
using Compat.readstring

# Command line "unzip" interface.

function unzip_tool_is_missing()
    try
        readstring(`unzip`)
        return false
    catch
        println("WARNING: unzip tool not found!")
        return true
    end
end

# Unzip file to dict using external "unzip" tool.

function test_unzip_file(z)
    r = Dict()
    for f in readlines(`unzip -Z1 $z`)
        f = chomp(f)
        if basename(f) != ""
            r[f] = readstring(`unzip -qc $z $f`)
        end
    end
    return r
end

# Unzip zip data to dict using external "unzip" tool.

function test_unzip(zip)
    mktemp((tmp,io)-> begin
        write(io, zip)
        close(io)
        return test_unzip_file(tmp)
    end)
end

dict = Dict("hello.txt"     => "Hello!\n",
            "foo/text.txt"  => "text\n")

# In memory ZIP from Dict...
@test dict == test_unzip(create_zip(dict))

@test dict == Dict(open_zip(create_zip(dict)))

@test open_zip(create_zip(dict))["hello.txt"] == "Hello!\n"

@test open_zip(create_zip("empty" => ""))["empty"] == ""

# In memory ZIP from pairs...
@test dict == test_unzip(create_zip("hello.txt"     => "Hello!\n",
                                   "foo/text.txt"  => "text\n"))

# In memory ZIP from tuples...
@test dict == test_unzip(create_zip(("hello.txt",     "Hello!\n"),
                                    ("foo/text.txt",  "text\n")))

# In memory ZIP from tuples...
@test dict == test_unzip(create_zip([("hello.txt",     "Hello!\n"),
                                    ("foo/text.txt",  "text\n")]))

# In memory ZIP from arrays...
@test dict == test_unzip(create_zip(["hello.txt", "foo/text.txt"],
                                    ["Hello!\n", "text\n"]))

# In memory ZIP using "do"...
zip_data = UInt8[]
open_zip(zip_data) do z
    z["hello.txt"] = "Hello!\n"
    z["foo/text.txt"] = "text\n"
end
@test dict == Dict(open_zip(zip_data))

# ZIP to file from Dict...
unzip_dict = ""
z = tempname() * ".zip"
#try
    create_zip(z, dict)

    @test unzip_tool_is_missing() || dict == test_unzip_file(z)

    @test open_zip(z)["hello.txt"] == "Hello!\n"

    @test dict == Dict(open_zip(z))

    create_zip(z, "foo" => "bar")

    @test open_zip(z)["foo"] == "bar"
    @test !haskey(open_zip(z), "hello.txt")

#finally
#    rm(z)
##end

# ZIP to file from Pairs...
unzip_dict = ""
z = tempname() * ".zip"
try
    create_zip(z, "hello.txt"     => "Hello!\n",
                  "foo/text.txt"  => "text\n")
    @test unzip_tool_is_missing() || dict == test_unzip_file(z)
    @test open_zip(z)["foo/text.txt"] == "text\n"
finally
    rm(z)
end

# Incremental ZIP to file...

f = tempname() * ".zip"
try
    open_zip(f) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
    end
    @test dict == Dict(open_zip(f))
    open_zip(f) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt"] == "text\n" # read again
    end

    # Add one file...
    open_zip(f) do z
        z["newfile"] = "new!\n"
    end
    open_zip(f) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["newfile"] == "new!\n"
        @test z["foo/text.txt"] == "text\n"
    end

    # Read and write (read first)...
    open_zip(f) do z
        z["hello.txt"] *= "World!\n"
        z["foo/bar.txt"] = "bar\n"
        z["foo/text.txt.copy"] = z["foo/text.txt"] * "Copy\n"
    end
    open_zip(f) do z
        @test z["hello.txt"] == "Hello!\nWorld!\n"
        @test z["foo/bar.txt"] == "bar\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n" # read again
        @test z["newfile"] == "new!\n"
    end

    # Read and write (write first)...
    open_zip(f) do z
        z["foo/bar.txt"] = "bar\n"
        z["foo/text.txt.copy"] = z["foo/text.txt"] * "Copy\n"
    end
    open_zip(f) do z
        @test z["hello.txt"] == "Hello!\nWorld!\n"
        @test z["foo/bar.txt"] == "bar\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n"
        @test z["newfile"] == "new!\n"
    end

finally
    rm(f)
end

# Write new file, then read...
f = tempname() * ".zip"
try
    open_zip(f) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test z["hello.txt"] == "Hello!\n"
    end
    @test dict == Dict(open_zip(f))

finally
    rm(f)
end

# Write new file, then iterate...
f = tempname() * ".zip"
try
    open_zip(f) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test [(n,v) for (n,v) in z] == [("hello.txt","Hello!\n"),("foo/text.txt","text\n")]
    end
    @test dict == Dict(open_zip(f))

finally
    rm(f)
end

# Write new file, then read and write...
f = tempname() * ".zip"
try
    open_zip(f) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["hello.txt"] == "Hello!\n"
        z["foo2/text.txt"] = "text\n"
    end
    open_zip(f) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo2/text.txt"] == "text\n"
        @test_throws KeyError z["foo3/text.txt"]
    end

finally
    rm(f)
end

# Incremental ZIP to buffer...

buf = UInt8[]

    open_zip(buf) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
    end
    @test dict == Dict(open_zip(buf))
    open_zip(buf) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt"] == "text\n" # read again
    end

    # Add one file...
    open_zip(buf) do z
        z["newfile"] = "new!\n"
    end
    open_zip(buf) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["newfile"] == "new!\n"
        @test z["foo/text.txt"] == "text\n"
    end

    # Read and write (read first)...
    open_zip(buf) do z
        z["hello.txt"] *= "World!\n"
        z["foo/bar.txt"] = "bar\n"
        z["foo/text.txt.copy"] = z["foo/text.txt"] * "Copy\n"
    end
    open_zip(buf) do z
        @test z["hello.txt"] == "Hello!\nWorld!\n"
        @test z["foo/bar.txt"] == "bar\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n" # read again
        @test z["newfile"] == "new!\n"
    end

    # Read and write (write first)...
    open_zip(buf) do z
        z["foo/bar.txt"] = "bar\n"
        z["foo/text.txt.copy"] = z["foo/text.txt"] * "Copy\n"
    end
    open_zip(buf) do z
        @test z["hello.txt"] == "Hello!\nWorld!\n"
        @test z["foo/bar.txt"] == "bar\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo/text.txt.copy"] == "text\nCopy\n"
        @test z["newfile"] == "new!\n"
    end

# Write new buffer, then read...
buf = UInt8[]

    open_zip(buf) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test z["hello.txt"] == "Hello!\n"
    end
    @test dict == Dict(open_zip(buf))

# Write new buffer, then iterate...
buf = UInt8[]

    open_zip(buf) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test [(n,v) for (n,v) in z] == [("hello.txt","Hello!\n"),("foo/text.txt","text\n")]
    end
    @test dict == Dict(open_zip(buf))

# Write new buffer, then read and write...
buf = UInt8[]

    open_zip(buf) do z
        z["hello.txt"] = "Hello!\n"
        z["foo/text.txt"] = "text\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["hello.txt"] == "Hello!\n"
        z["foo2/text.txt"] = "text\n"
    end
    open_zip(buf) do z
        @test z["hello.txt"] == "Hello!\n"
        @test z["foo/text.txt"] == "text\n"
        @test z["foo2/text.txt"] == "text\n"
    end

# Unzip file created by command-line "zip" tool...

testzip = joinpath(Pkg.dir("InfoZIP"),"test","test.zip")
d = Dict(open_zip(testzip))
@test sum(d["test.png"]) == 462242
delete!(d, "test.png")
@test dict == d

# unzip()...

mktempdir() do d
    InfoZIP.unzip(testzip, d)
    @test readstring(joinpath(d, "hello.txt")) == "Hello!\n"
    @test readstring(joinpath(d, "foo/text.txt")) == "text\n"
end

mktempdir() do d
    InfoZIP.unzip(create_zip(dict), d)
    @test readstring(joinpath(d, "hello.txt")) == "Hello!\n"
    @test readstring(joinpath(d, "foo/text.txt")) == "text\n"
end

mktempdir() do d
    cd(d) do
        write("hello.txt", "Hello!\n")
        mkdir("foo")
        write("foo/text.txt", "text\n")
        @test dict == Dict(open_zip(create_zip(["hello.txt", "foo/text.txt"])))
    end
end


text:
@testset "Basis Set" begin
    @set printstyle none
    bfH = Fermi.GaussianBasis.BasisFunction(0,
    [0.9817067282501578, 0.9494640079478909, 0.2959064596969707],
    [3.425250914, 0.6239137298, 0.168855404])

    # Read a simple basis function
    x = Fermi.GaussianBasis.read_basisset("sto-3g", "H")[1]
    @test begin
        x.l == bfH.l &&
        x.coef == bfH.coef &&
        x.exp == bfH.exp
    end

    # Test error: Invalid Basis set
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.read_basisset("invalid_basis", "H")

    # Test error: Invalid Atom
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.read_basisset("sto-3g", "U")

    # Test some basis set non sense
    x = """W    3   1.00
    0.3425250914D+01       0.1543289673D+0
    0.6239137298D+00       0.5353281423D+00
    0.1688554040D+00       0.4446345422D+00"""
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.basis_from_string(x)
    x = """S    3   1.00
    0.3425250914D+01       0.1543289673D+0
    0.6239137298D+00       Ops
    0.1688554040D+00       0.4446345422D+00"""
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.basis_from_string(x)

    x =  """SP D   3   1.00
          0.5033151319D+01      -0.9996722919D-01       0.1559162750D+00
          0.1169596125D+01       0.3995128261D+00       0.6076837186D+00
          0.3803889600D+00       0.7001154689D+00       0.3919573931D+00"""
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.two_basis_from_string(x)

    x =  """Some   3   1.00
          0.5033151319D+01      -0.9996722919D-01       0.1559162750D+00
          0.1169596125D+01       0.3995128261D+00       0.6076837186D+00
          0.3803889600D+00       0.7001154689D+00       0.3919573931D+00"""
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.two_basis_from_string(x)

    x =  """SP   3   1.00
          0.5033151319D+01      -0.9996722919D-01       0.1559162750D+00
          0000000000000000       0.3995128261D+00       0.6076837186D+00
          0.3803889600D+00       0.7001154689D+00       0.3919573931D+00"""
    @test_throws Fermi.Options.FermiException Fermi.GaussianBasis.two_basis_from_string(x)

    @set basis sto-3g

    bs = Fermi.GaussianBasis.BasisSet()
    x = Fermi.GaussianBasis.string_repr(bs)
    @test begin
          occursin("sto-3g Basis Set", x) &&
          occursin(r"Number of shells:\s+?5", x) &&
          occursin(r"Number of basis:\s+?7", x) &&
          occursin("O: 1s 2s 1p", x) &&
          occursin("H: 1s", x)
    end
    x = Fermi.GaussianBasis.string_repr(bs[1,1])
    @test begin
          occursin("S shell with 1 basis built from 3 primitive gaussians", x) &&
          occursin(r"χ₀₀\s+?=\s+?\d+?\.\d+?⋅Y₀₀⋅exp\(-\d+?\.\d+?⋅r²\)", x) &&
          occursin(r"[+-]\s+?\d+?\.\d+?⋅Y₀₀⋅exp\(-\d+?\.\d+?⋅r²\)", x)
    end
    x = Fermi.GaussianBasis.string_repr(bs[1][3])
    @test begin
          occursin("P shell with 3 basis built from 3 primitive gaussians", x) &&
          occursin(r"χ.{2,3}\s+?=\s+?\d+?\.\d+?⋅Y.{2,3}⋅r¹⋅exp\(-\d+?\.\d+?⋅r²\)", x) &&
          occursin(r"[+-]\s+?\d+?\.\d+?⋅Y.{2,3}⋅r¹⋅exp\(-\d+?\.\d+?⋅r²\)", x)
    end
end


text:
lazy_dir = joinpath(@__DIR__, "LazilyInitializedFields")
lazy_file = joinpath(lazy_dir, "LazilyInitializedFields.jl")
mkpath(lazy_dir)
download("https://raw.githubusercontent.com/KristofferC/LazilyInitializedFields.jl/master/src/LazilyInitializedFields.jl",
         lazy_file)


text:
using ReverseDiffExpressions, LoopVectorization
using ReverseDiffExpressions.LoopSetDerivatives
using Test

@testset "ReverseDiffExpressions.jl" begin

    AmulBq = :(for m ∈ 1:M, n ∈ 1:N
                C[m,n] = zero(eltype(B))
                for k ∈ 1:K
                C[m,n] += A[m,k] * B[k,n]
                end
                end);

    lsAmulB = LoopVectorization.LoopSet(AmulBq);

    tracked_vars = Set([:A, :B]);
    ∂lsAmulB = ReverseDiffExpressions.LoopSetDerivatives.∂LoopSet(lsAmulB, tracked_vars);
    ∂lsAmulB.fls
    ∂lsAmulB.rls
    # Write your own tests here.
end


text:
# TODO: create FoldsBase and use it here

struct Break{T}
    value::T
end

# unbreak(v) = v
# unbreak(b::Break) = b.value

function foldl_modules_rec(rf, init, m::Module)
    acc = rf(init, m)
    acc isa Break && return acc

    for n in sort!(names(m; all = true))
        sub = try
            getproperty(m, n)
        catch
            continue
        end
        sub isa Module || continue
        parentmodule(sub) === m || continue

        acc = rf(acc, sub)
        acc isa Break && return acc
    end

    return acc
end


text:
#---------------------------------------------------------------------#
#This code computes the element mass and differentiation matrices
#Written by F.X. Giraldo on April 24, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function element_matrices(ψ,dψ,Np,Nq,ωq,DFloat)

    #Initialize Matrices
    M=zeros(DFloat,Np,Np)
    D=zeros(DFloat,Np,Np)
    Dw=zeros(DFloat,Np,Np)
    F=zeros(DFloat,Np,Np)

    #Construct Mass and Differentiation Matrices
    @inbounds for k=1:Nq, j=1:Np, i=1:Np
        M[i,j]+=ωq[k]*ψ[i,k]*ψ[j,k]
        D[i,j]+=ωq[k]*ψ[i,k]*dψ[j,k]
        Dw[i,j]+=ωq[k]*dψ[i,k]*ψ[j,k]
    end
    F[1,1]=-1; F[Np,Np]=+1

    return (M,D,Dw,F)
end