text:
module MPITestImages

using DocStringExtensions
using Pkg.Artifacts
using FileIO
using Images

const remotefiles = [
	"Phantom1",
	"Phantom2",
	"Phantom3",
	"Phantom3_small",
	"phantom_2delta_notInLine_3px",
	"phantom_3delta_InLine_3px",
	"phantom_3delta_notInLine_3px",
	"phantom_brain_30",
	"phantom_checkers_4px",
	"phantom_delta_1px",
	"phantom_delta_3px",
	"phantom_delta_3px_shift1",
	"phantom_delta_3px_shift2",
	"phantom_vessel",
	"phantom_vessel2"
]

onTheFlyImages = Vector{Symbol}()

"""
Add a function symbol to the list of known test image generation functions.
"""
function addOnTheFlyImage(fun::Symbol)
	if !(fun in onTheFlyImages)
		if !(fun in remotefiles)
			push!(onTheFlyImages, fun)
		else
			@warn "The on-the-fly image `$(string(fun))` is already linked with a remote file."
		end
	else
		error("The on-the-fly image `$(string(fun))` has already been added.")
	end
end

export testimage_gen
"""
Macro for annotating functions that can be used to generate test images.
"""
macro testimage_gen(expr::Expr)
	addOnTheFlyImage(expr.args[1].args[1]) # Note: use dump to debug expr
	return expr
end

export testimage
"""
		$(SIGNATURES)

Retrieve a test image with the given `name` and the matching parameters.

Note: The name must correspond either to a remote file or a function
			name annotated by the `testimage_gen` macro. If both exist, precedence
			is given to the function.

# Examples
```jldoctest
julia> image = testimage("delta_image", (8, 8), 2; sizeOfPoint=(3, 2), distanceOfPoints=(x -> 0, x -> 4), pivot=(3, 3))
8×8 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function testimage(name::String, args...; kwargs...)
	if Symbol(name) in onTheFlyImages # Prioritize on-the-fly images over file-based ones
		f = getfield(@__MODULE__, Symbol(name))
		return f(args...; kwargs...)
	elseif name in remotefiles
		rootpath = artifact"testimages"

		try
			return changeScale(name, rootpath, args...)
		catch err
			if isa(err, MethodError)
				@warn "No scale provided. Using default image size."
			end
		end

		return load(joinpath(rootpath, "phantoms", name*".png"))
	else
		error("The given name `$name` did not match a known test image.")
	end
end

"""
Loads the specified image from remote source and scales it accordingly.
"""
function changeScale(name::String, rootpath::String, size::Tuple{Integer, Integer}, args...)
	image = load(joinpath(rootpath, "phantoms", name*".png"))

	return imresize(image, size)
end

export TestImage
"""
Struct describing a testimage.
"""
struct TestImage
	name::String
	data::AbstractArray
	args::Tuple
	kwargs::Base.Pairs

	function TestImage(name::String, args...; kwargs...)
		data = testimage(name, args...; kwargs...)
		return new(name, data, args, kwargs)
	end
end

export name
name(img::TestImage) = img.name

export data
data(img::TestImage) = img.data

include("OnTheFly.jl")

end


text:
"""
    output_type(problem)

Get type returned by `solve(problem)`.
"""
function output_type end

output_type(::BTPDE{T}) where {T} = Vector{Complex{T}}
output_type(::HADC{T}) where {T} = Vector{T}
output_type(::Karger{T}) where {T} = Complex{T}
output_type(::Laplace{T}) where {T} = NamedTuple
output_type(::MatrixFormalism{T}) where {T} = Vector{Complex{T}}
output_type(::AnalyticalLaplace{T}) where {T} = NamedTuple
output_type(::AnalyticalMatrixFormalism{T}) where {T} = Complex{T}

ext: jl


text:
mutable struct ProbabilityVector{T, d} <: QuantumVariable{T}
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int, Int}
    value::Convex.ValueOrNothing
    vexity::Convex.Vexity
    function ProbabilityVector{T}(d::Int) where {T <: Real}
        this = new{T, d}(:ProbabilityVector, 0, (d,1), nothing, Convex.AffineVexity())
        this.id_hash = objectid(this)
        Convex.id_to_variables[this.id_hash] = this
        this
    end
    ProbabilityVector(d::Int) = ProbabilityVector{Float64}(d)
end

Convex.constraints(p::ProbabilityVector) = [ p >= 0, sum(p) == 1 ]
Convex.sign(::ProbabilityVector) = Convex.Positive()
Convex.vartype(::ProbabilityVector) = Convex.ContVar


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   General structures related to orbits.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export Orbit, KeplerianElements, OrbitStateVector

"""
    Orbit

Abstract type of an orbit representation.

"""
abstract type Orbit{T} end

"""
    KeplerianElements{T1,T2}

This structure defines the orbit in terms of the Keplerian elements.

# Fields

* `t`: Epoch.
* `a`: Semi-major axis [m].
* `e`: Eccentricity [ ].
* `i`: Inclination [rad].
* `Ω`: Right ascension of the ascending node [rad].
* `ω`: Argument of perigee [rad].
* `f`: True anomaly [rad].

"""
@with_kw_noshow struct KeplerianElements{T} <: Orbit{T}
    t::T
    a::T
    e::T
    i::T
    Ω::T
    ω::T
    f::T
end

function KeplerianElements(t::T1, a::T2, e::T3, i::T4, Ω::T5, ω::T6, f::T7) where
    {T1<:Number, T2<:Number, T3<:Number, T4<:Number, T5<:Number, T6<:Number, T7<:Number}
    T = promote_type(T1, T2, T3, T4, T5, T6, T7)
    return KeplerianElements{T}(t,a,e,i,Ω,ω,f)
end

"""
    OrbitStateVector{T}

Store the state vector representation of an orbit.

# Fields

* `t`: Epoch [Julian Day].
* `r`: Position vector [m].
* `v`: Velocity vector [m/s].
* `a`: Acceleration vector [m/s²].

"""
@with_kw_noshow struct OrbitStateVector{T} <: Orbit{T}
    t::T            = 0
    r::SVector{3,T} = SVector{3,T}(0,0,0)
    v::SVector{3,T} = SVector{3,T}(0,0,0)
    a::SVector{3,T} = SVector{3,T}(0,0,0)
end


text:
using BinaryBuilder, Pkg.Types

name = "AFLplusplus"
version = v"2.68"
sources = [
    ArchiveSource("https://github.com/AFLplusplus/AFLplusplus/archive/2.68c.tar.gz", "862e155c97737770baa26ffedf324a7fa255b757c85b0c9a6f312264f2ca29c5"),
    ArchiveSource("http://download.qemu-project.org/qemu-3.1.0.tar.xz", "6a0508df079a0a33c2487ca936a56c12122f105b8a96a44374704bef6c69abfc", unpack_target="AFLplusplus-2.68c/qemu_mode"),
    DirectorySource("bundled"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/AFLplusplus-2.68c

# Patch the makefile to not use llvm-config. The llvm-config binary from Clang_jll cannot be used when cross compiling
atomic_patch -p1 ${WORKSPACE}/srcdir/llvm_mode_makefile.patch

export REAL_CC="clang"
export REAL_CXX="clang++"
export AFL_CC="$prefix/tools/clang"
export AFL_CXX="$prefix/tools/clang++"
export AFL_PATH="$libdir/afl"
export PREFIX=$prefix

cd llvm_mode
make
make install

install_license ${WORKSPACE}/srcdir/AFLplusplus-2.68c/LICENSE

# QEMU mode is only working on Linux glibc platforms
if [[ "${target}" != *-freebsd* ]] && [[ "${target}" != *-apple-* ]] && [[ "${target}" != *-*-musl ]]; then
    cd ../qemu_mode/qemu-3.1.0/

	# apply patches from AFL's qemu_mode
	atomic_patch -p1 ../patches/elfload.diff
	atomic_patch -p1 ../patches/mips-fpu.diff
	atomic_patch -p1 ../patches/bsd-elfload.diff
	atomic_patch -p1 ../patches/cpu-exec.diff
	atomic_patch -p1 ../patches/syscall.diff
	atomic_patch -p1 ../patches/translate-all.diff
	atomic_patch -p1 ../patches/tcg.diff
	atomic_patch -p1 ../patches/i386-translate.diff
	atomic_patch -p1 ../patches/arm-translate.diff
	atomic_patch -p1 ../patches/arm-translate-a64.diff
	atomic_patch -p1 ../patches/i386-ops_sse.diff
	atomic_patch -p1 ../patches/i386-fpu_helper.diff
	atomic_patch -p1 ../patches/softfloat.diff
	atomic_patch -p1 ../patches/configure.diff
	atomic_patch -p1 ../patches/tcg-runtime.diff
	atomic_patch -p1 ../patches/tcg-runtime-head.diff
	atomic_patch -p1 ../patches/translator.diff
	atomic_patch -p1 ../patches/__init__.py.diff
	atomic_patch -p1 ../patches/make_strncpy_safe.diff
	atomic_patch -p1 ../patches/mmap_fixes.diff

	./configure --disable-system \
		--enable-linux-user --disable-gtk --disable-sdl --disable-vnc --enable-capstone=internal \
		--target-list="i386-linux-user,x86_64-linux-user,arm-linux-user,aarch64-linux-user,ppc-linux-user" --enable-pie --python="/usr/bin/python3"

	make -j${nproc}

	install i386-linux-user/qemu-i386 $bindir/afl-qemu-trace-i386
	install x86_64-linux-user/qemu-x86_64 $bindir/afl-qemu-trace-x86_64
	install arm-linux-user/qemu-arm $bindir/afl-qemu-trace-arm
	install aarch64-linux-user/qemu-aarch64 $bindir/afl-qemu-trace-aarch64
	install ppc-linux-user/qemu-ppc $bindir/afl-qemu-trace-ppc
fi
"""

# Windows is not supported by AFL, ARMv7l musl won't build
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("i686", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("x86_64", "freebsd"),
    Platform("x86_64", "macos"),
]
platforms = expand_cxxstring_abis(platforms)

products = [
    ExecutableProduct("afl-clang-fast", :afl_clang_fast),
    ExecutableProduct("afl-clang-fast++", :afl_clang_fast_cxx),
]

dependencies = [
	Dependency("Clang_jll", v"9.0.1"; compat="9.0.1"),
	Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
	Dependency("Pixman_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Neurthino

@test Neurthino.tracklength(0, 2.0) ≈ 2.0 atol=0.01
@test Neurthino.tracklength(pi/2, 1.0) ≈ 112.88 atol=0.01

@test Neurthino.tracklength(pi*1/3, Neurthino.EARTH_RADIUS / 2) ≈ 4149.99 atol=0.01
@test Neurthino.tracklength(0.7, Neurthino.EARTH_RADIUS / 2) ≈ 3595.04 atol=0.01
@test Neurthino.tracklength(pi/6, Neurthino.EARTH_RADIUS / 2) ≈ 3409.97 atol=0.01
@test Neurthino.tracklength(1.91, Neurthino.EARTH_RADIUS / 2) ≈ 6678.27 atol=0.01
@test Neurthino.tracklength(2.44, Neurthino.EARTH_RADIUS / 2) ≈ 8463.26 atol=0.01

@test Neurthino.tracklength(2.62, 1860) ≈ 9872.5 atol=1.0
@test Neurthino.tracklength(2.97, 1860) ≈ 10769.02 atol=1.0
@test Neurthino.tracklength(2.27, 5600) ≈ 6839.15 atol=1.0

path = Neurthino.prempath(pi, 1.5)

@test path.baseline[1] ≈ 128.69 atol=0.1

@test path.density[1] ≈ 3.38 atol=0.1
@test path.density[99] ≈ 3.38 atol=0.1
@test path.density[50] ≈ 13.09 atol=0.1
@test path.density[20] ≈ 5.38 atol=0.1
@test path.density[80] ≈ 5.38 atol=0.1


text:
using HorizonSideRobots

function paint_perimeter(robot::Robot)
    c = 0
    while !isborder(robot, West)
        move!(robot, West)
        c += 1
    end

    sides = [Nord, Ost, Sud, West]
    for side in sides
        while !isborder(robot, side)
            putmarker!(robot)
            move!(robot, side)
        end
    end

    while !ismarker(robot)
        putmarker!(robot)
        move!(robot, Nord)
    end

    while c > 0
        move!(robot, Ost)
        c -= 1
    end
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
module Imagine

using ImagineInterface, Unitful, Distributed
import ImagineWorker

#ENV["JULIA_PROJECT"] = joinpath(@__DIR__, "..")

@static if (Sys.iswindows() && isfile("C:\\Windows\\System32\\nicaiu.dll"))
    using NIDAQ #just for the trigger pulse
    const NPERSISTENT_WORKERS = 4 #one each for AI, AO, DI, DO, one for PFI pulse
    devs = NIDAQ.devices()
    if isempty(devs)
	error("No NI devices detected")
    end
    DEFAULT_DEVICE = devs[1] * "/" #"Dev1/"
    FREE_WORKERS = Int[]
    USED_WORKERS = Int[]
    WORKERS = Int[]

    export run_imagine,
            set_device,
            record_signals,
            output_signals

    include("manage_workers.jl")
    include("io.jl")
    include("run.jl")
end

@static if !Sys.iswindows()
    @warn("This package only works on Windows")
end

@static if (Sys.iswindows() && !isfile("C:\\Windows\\System32\\nicaiu.dll"))
    @warn("nicaiu.dll not found")
end

end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ e708e17d-617e-4a3f-bca1-164d22c3dffc
begin
    import Pkg
    Pkg.activate("..")
end

# ╔═╡ b0dbee5a-9881-48be-9fe8-56cc4f09f4dd
begin
    using Revise
    using Hydroelastics

    using MeshCat, Colors
    using GeometryBasics
    using CoordinateTransformations
    using LinearAlgebra, StaticArrays
    using PlutoUI
end

# ╔═╡ 8bebb189-d470-4077-88dd-a9ad97b94b76
@bind d Slider(0.0:0.01:1.0; show_value = true, default = 0.5)

# ╔═╡ 27b6e401-897c-459a-b6f6-575f3f89b0b5
begin
    vis = Visualizer()
    IJuliaCell(vis)
end

# ╔═╡ 90ff5a44-0464-4e10-8494-99a91b301ed4
sphere = make_icosphere(1)

# ╔═╡ 51f5ad18-2172-41b9-ba71-96e0e1871eff
cube = make_cube()

# ╔═╡ c9706e48-450e-4d1f-8541-737d7e0aa470
begin
    object1 = translate(sphere, @SVector [0.031, -d, 0.052])
    object2 = translate(sphere, @SVector [0, d, 0])
end

# ╔═╡ e32777f6-d54f-4ef1-be3f-0c57c1d01cd0
compute_force(object1, object2)

# ╔═╡ 539ed9a3-71d4-4da6-bc46-743718018abd
begin
    local material1 = MeshPhongMaterial(color = RGBA(0.3, 0.1, 0.7, 1.0))
    setobject!(vis[:object1], render_obj(object1), material1)
    local material2 = MeshPhongMaterial(color = RGBA(0.8, 0.3, 0.4, 1.0))
    setobject!(vis[:object2], render_obj(object2), material2)
end

# ╔═╡ f63d1c23-78cc-489c-8acc-62c7d0aca7f3
begin
    import Hydroelastics: HalfPlane, Point, intersect_halfplanes
    import Plots: plot, quiver!, plot!, scatter!

    halfplanes = [
        HalfPlane(Point(0.05328478907262187, 0.032931810730249304), Point(-0.5257311121191336, 0.85065080835204)),
        HalfPlane(Point(-0.1780353002260096, 0.2880671669629754), Point(0.85065080835204, 0.5257311121191336)),
        HalfPlane(Point(-0.3302474442330619, 1.1330056073001064e-17), Point(-3.4307778215552914e-17, -1.0)),
        HalfPlane(Point(-0.0, 1.0685339887498948), Point(-1.0, -0.0)),
        HalfPlane(Point(0.007597871376374749, -0.004695742752749542), Point(-0.5257311121191337, -0.8506508083520398)),
        HalfPlane(Point(0.1547801932600118, 0.2504396134799765), Point(0.85065080835204, -0.5257311121191336)),
        HalfPlane(Point(0.3501096768180673, -1.2011485143393161e-17), Point(3.4307778215552914e-17, 1.0)),
        HalfPlane(Point(1.919491578940921e-16, 1.0685339887498948), Point(-1.0, 1.7963785889362146e-16)),
    ]
    halfplanes = [
        HalfPlane(Point(-0.16469227056401645, -0.2664776914569725), Point(-0.8506508083520399, 0.5257311121191335)),
        HalfPlane(Point(0.028287259696359788, -0.040547756753354335), Point(0.8201447309019166, 0.5721561153861965)),
        HalfPlane(Point(-0.16620458460952442, 0.09547326678562677), Point(-0.49810107091297723, -0.8671189786617204)),
        HalfPlane(Point(0.2970629430715192, -0.23999299562154927), Point(0.6284282897093981, 0.77786752387082)),
        HalfPlane(Point(-0.05197402131154934, 0.06002998491806123), Point(-0.7560128937297079, -0.6545567236797993)),
        HalfPlane(Point(0.1879473775300143, 0.30410524493997143), Point(0.85065080835204, -0.5257311121191336)),
        HalfPlane(Point(0.1619918586413753, -0.0851641599977757), Point(0.46534112719498627, 0.8851314226378474)),
        HalfPlane(Point(0.2970629430715192, -0.23999299562154938), Point(0.6284282897093982, 0.7778675238708199)),
    ]

    plot()
    for hp in halfplanes
        x, y = hp.p
        x1, y1 = hp.pq
        quiver!([x - 0.5 * x1], [y - 0.5 * y1]; quiver = ([x1], [y1]))
    end

    points = intersect_halfplanes(halfplanes)
    if !isnothing(points)
        scatter!(points[1, :], points[2, :])
    end

    plot!(legend = :none)
end

# ╔═╡ Cell order:
# ╠═e708e17d-617e-4a3f-bca1-164d22c3dffc
# ╠═b0dbee5a-9881-48be-9fe8-56cc4f09f4dd
# ╠═8bebb189-d470-4077-88dd-a9ad97b94b76
# ╠═e32777f6-d54f-4ef1-be3f-0c57c1d01cd0
# ╠═27b6e401-897c-459a-b6f6-575f3f89b0b5
# ╠═90ff5a44-0464-4e10-8494-99a91b301ed4
# ╠═51f5ad18-2172-41b9-ba71-96e0e1871eff
# ╠═c9706e48-450e-4d1f-8541-737d7e0aa470
# ╠═539ed9a3-71d4-4da6-bc46-743718018abd
# ╠═f63d1c23-78cc-489c-8acc-62c7d0aca7f3

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
using IncrementalInference
using Test

##

@testset "saving to and loading from FileDFG" begin
##

fg = generateCanonicalFG_Kaess()
addVariable!(fg, :x4, ContinuousScalar)
addFactor!(fg, [:x2;:x3;:x4], LinearRelative(Normal()), multihypo=[1.0;0.6;0.4])

saveFolder = "/tmp/dfg_test"
saveDFG(fg, saveFolder)
# VERSION above 1.0.x hack required since Julia 1.0 does not seem to havfunction `splitpath`

retDFG = initfg()
retDFG = loadDFG!(retDFG, saveFolder)
Base.rm(saveFolder*".tar.gz")

@test symdiff(ls(fg), ls(retDFG)) == []
@test symdiff(lsf(fg), lsf(retDFG)) == []

@show getFactor(fg, :x2x3x4f1).solverData.multihypo
@show getFactor(retDFG, :x2x3x4f1).solverData.multihypo

# check for match
@test getFactor(fg, :x2x3x4f1).solverData.multihypo - getFactor(retDFG, :x2x3x4f1).solverData.multihypo |> norm < 1e-10
@test getFactor(fg, :x2x3x4f1).solverData.certainhypo - getFactor(retDFG, :x2x3x4f1).solverData.certainhypo |> norm < 1e-10

##
end

@testset "saving to and loading from FileDFG with nullhypo, eliminated, solveInProgress" begin
##

fg = generateCanonicalFG_Kaess()
getSolverParams(fg).attemptGradients = true

addVariable!(fg, :x4, ContinuousScalar)
addFactor!(fg, [:x2;:x3;:x4], LinearRelative(Normal()), multihypo=[1.0;0.6;0.4])
addFactor!(fg, [:x1;], Prior(Normal(10,1)), nullhypo=0.5)

solveTree!(fg)

#manually change a few fields to test if they are preserved
fa = getFactor(fg, :x2x3x4f1)
fa.solverData.eliminated = true
fa.solverData.solveInProgress = 1
fa.solverData.nullhypo = 0.5

saveFolder = "/tmp/dfg_test"
saveDFG(fg, saveFolder)

retDFG = initfg()
loadDFG!(retDFG, saveFolder)
Base.rm(saveFolder*".tar.gz")

@test issetequal(ls(fg), ls(retDFG))
@test issetequal(lsf(fg), lsf(retDFG))

@show getFactor(fg, :x2x3x4f1).solverData.multihypo
@show getFactor(retDFG, :x2x3x4f1).solverData.multihypo

# check for match
@test isapprox(getFactor(fg, :x2x3x4f1).solverData.multihypo, getFactor(retDFG, :x2x3x4f1).solverData.multihypo)
@test isapprox(getFactor(fg, :x2x3x4f1).solverData.certainhypo, getFactor(retDFG, :x2x3x4f1).solverData.certainhypo)

fb = getFactor(retDFG, :x2x3x4f1)
@test fa == fb

fa = getFactor(fg, :x1f2)
fb = getFactor(retDFG, :x1f2)

@test fa == fb

##
end


text:
using Random
using Parameters
#using Debugger
using DataStructures

Random.seed!(1)

function actUCS(mdp::ActMdp, s::Vector{Float64})::Tuple{Float64, Vector{Float64}, Vector{Float64}}
	# frontier(state, pastCost)
	frontier = PriorityQueue{Vector{Float64} , Float64}()
	# explored[s] = (pastCost, prev_s, prev_a, prev_cost)
	explored = Dict{Vector{Float64}, Tuple{Float64, Vector{Float64}, Float64, Float64}}()
	# previous[sp] = (s, a, cost)
	previous = Dict{Vector{Float64}, Tuple{Vector{Float64}, Float64, Float64}}()

	#sizehint!(frontier, 1500000)
	sizehint!(explored, 1_500_000)
	sizehint!(previous, 1_500_000)

	frontier[s] = 0.0
	previous[s] = ([0., 0., -1.0], 0.0, 0.0)

	while true
		s, pastCost = dequeue_pair!(frontier)
		prev_s, prev_a, prev_cost = previous[s]
		explored[s] = (pastCost, prev_s, prev_a, prev_cost)
		if isEnd(mdp, s) > 0
			global minCost = pastCost
			break
		end
		for res in succAndCost(mdp, s)
			a, sp, cost = res
			if haskey(explored, sp)
				continue
			end
			newCost = pastCost + cost
			if !haskey(frontier, sp) || newCost < frontier[sp]
				frontier[sp] = newCost
				previous[sp] = (s, a, cost)
			end
		end
	end

	#sf = deepcopy(s)
	sf = copy(s)

	# recover history
	println("n_states explored : ", length(explored))
	history = Tuple{Float64, Float64, Vector{Float64}}[]
	sequence = Float64[]
	while s[3] > 0 # !isnothing(s)
		pastCost, prev_s, prev_a, prev_cost = explored[s]
		if prev_s[3] > 0 # !isnothing(prev_s)
			push!(history, (prev_a, pastCost, prev_s))
			push!(sequence, prev_a)
		end
		s = prev_s
	end

	return minCost, reverse(sequence), sf
end

function actDP(mdp::ActMdp, s::Vector{Float64})
	futureCost = Dict{Vector{Float64}, Tuple{Float64, Float64, Vector{Float64}, Float64}}()

	sizehint!(futureCost, 2500000)

	function recurse(mdp::ActMdp, s::Vector{Float64})::Float64
		if haskey(futureCost, s)
			return futureCost[s][1]
		end
		if isEnd(mdp, s) > 0
			return 0
		end
		succ = []
		for res in succAndCost(mdp, s)
			a, sp, cost = res
			push!(succ, (cost + recurse(mdp, sp), a, sp, cost))
		end
		futureCost[s] = minimum(succ)
		#futureCost[s] = minimum([(cost + recurse(sp), a, sp, cost) for a, sp, cost in succAndCost(mdp, s)])
		return futureCost[s][1]
	end

	minCost = recurse(mdp, s)

	# recover history
	println("n_states explored : ", length(futureCost))
	history = []
	sequence = []
	while isEnd(mdp, s) == 0
		totalCost, a, sp, cost = futureCost[s]
		push!(history, (a, totalCost, s))
		push!(sequence, a)
		s = sp
	end

	#println("cost: ", minCost)
	#println("history: ", history)
	return minCost, sequence, history[end][3]
end


text:
"""

The exported names are
$(EXPORTS)
"""
module Algebra

using NLsolve
using DocStringExtensions
using SparseArrays
using LinearAlgebra
using Test
using SparseMatricesCSR

using Gridap.Helpers

import Base: convert, size, getindex, show, count, *
import LinearAlgebra: mul!
import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals

export length_to_ptrs!
export rewind_ptrs!
export allocate_vector
export allocate_matrix
export allocate_matrix_and_vector
export allocate_in_domain
export allocate_in_range
export add_entries!
export muladd!
export nz_counter
export nz_allocation
export create_from_nz
export LoopStyle
export Loop
export DoNotLoop
export ArrayBuilder
export SparseMatrixBuilder
export MinMemory
export MinCPU
export get_array_type

export nz_index
export is_entry_stored
export finalize_coo!
export sparse_from_coo
export add_entry!
export copy_entries!
export allocate_coo_vectors
export push_coo!

export LinearSolver
export SymbolicSetup
export NumericalSetup
export solve
export solve!
export symbolic_setup
export numerical_setup
export numerical_setup!
export test_linear_solver

export LUSolver
export BackslashSolver

export NonlinearOperator
export residual!
export residual
export jacobian!
export jacobian
export hessian
export hessian!
export residual_and_jacobian!
export residual_and_jacobian
export allocate_residual
export allocate_jacobian
export allocate_residual_and_jacobian
export zero_initial_guess
export test_nonlinear_operator

export NonlinearSolver
export test_nonlinear_solver

export NewtonRaphsonSolver
export NLSolver

export AffineOperator
export get_matrix
export get_vector

include("AlgebraInterfaces.jl")

include("SparseMatrixCSC.jl")

include("SparseMatrixCSR.jl")

include("SymSparseMatrixCSR.jl")

include("NonlinearOperators.jl")

include("NonlinearSolvers.jl")

include("LinearSolvers.jl")

include("NLSolvers.jl")

end # module


text:
using DashBootstrapComponents, DashHtmlComponents

popover_children = [
    dbc_popoverheader("Popover header"),
    dbc_popoverbody("And here's some amazing content. Cool!"),
];

popovers = html_div([
    dbc_button(
        "Click",
        id="click-target",
        color="danger",
        className="mr-1",
        n_clicks=0,
    ),
    dbc_popover(popover_children, id="click", target="click-target", trigger="click"),
    dbc_button(
        "Focus",
        id="focus-target",
        color="danger",
        className="mr-1",
        n_clicks=0,
    ),
    dbc_popover(popover_children, id="focus", target="focus-target", trigger="focus"),
    dbc_button(
        "Hover",
        id="hover-target",
        color="danger",
        className="mr-1",
        n_clicks=0,
    ),
    dbc_popover(popover_children, id="hover", target="hover-target", trigger="hover"),
    dbc_button("Legacy", id="legacy-target", color="danger", n_clicks=0),
    dbc_popover(
        popover_children,
        id="legacy",
        target="legacy-target",
        trigger="legacy",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module ONS
using JSON, Requests

struct Contact
    email::String
    name::String
    telephone::String
end
Contact(d::Dict) = Contact(get(d, "email", ""),
get(d, "name", ""),
get(d, "telephone", ""))

struct MetaDataDescription
    title::String
    summary::String
    keywords::Vector{String}
    metaDescription::String
    nationalStatistics::Bool
    contact::Contact
    releaseDate::String
    nextRelease::String
    edition::String
    datasetId::String
    datasetUri::String
    cdid::String
    unit::String
    preUnit::String
    source::String
    date::String
    number::String
    keyNote::String
    sampleSize::String
end
MetaDataDescription(d::Dict) = MetaDataDescription(
    get(d, "title", ""),
    get(d, "summary", ""),
    get(d, "keywords", []),
    get(d, "metaDescription", ""),
    get(d, "nationalStatistics", false),
    haskey(d, "contact") ? Contact(d["contact"]) :  empty(Contact),
    get(d, "releaseDate", ""),
    get(d, "nextRelease", ""),
    get(d, "edition", ""),
    get(d, "datasetId", ""),
    get(d, "datasetUri", ""),
    get(d, "cdid", ""),
    get(d, "unit", ""),
    get(d, "preUnit", ""),
    get(d, "source", ""),
    get(d, "date", ""),
    get(d, "number", ""),
    get(d, "keyNote", ""),
    get(d, "sampleSize", "")
)

struct MetaData
    uri::String
    typ::String
    description::MetaDataDescription
    searchBoost::Vector{String}
end
MetaData(d::Dict) = MetaData(
    get(d, "uri", ""),
    get(d, "type", ""),
    haskey(d, "description") ? MetaDataDescription(d["description"]) : empty(MetaDataDescription),
    get(d, "searchBoost", String[])
)
function Base.display(r::MetaData)
    !isempty(r.description.title) ? println(r.description.title) : println("No title")
    println("  dataset    : ", r.description.datasetId)
    println("  timeseries : ", r.description.cdid)
    println()
end

function Base.display(R::Vector{MetaData})
    for (i, r) in enumerate(R)
        println("$i.")
        display(r)
        println()
    end
end

struct Records
    items::Vector{MetaData}
    itemsPerPage::Int
    startIndex::Int
    totalItems::Int
end
Records(d::Dict) = Records(
    haskey(d, "items") ? MetaData.(d["items"]) : MetaData[],
    get(d, "itemsPerPage", typemin(Int)),
    get(d, "startIndex", typemin(Int)),
    get(d, "totalItems", typemin(Int))
)

function Base.display(R::Records)
    println("Total items: ", R.totalItems)
    println("Displaying: ", R.startIndex + (1:R.itemsPerPage))
    println()
    display(R.items)
end
Base.getindex(R::Records, i) = R.items[i]

struct Period
    date::String
    value::String
    label::String
    year::String
    month::String
    quarter::String
    sourceDataset::String
    updateDate::String
end
Period(d::Dict) = Period(
    get(d, "date", ""),
    get(d, "value", ""),
    get(d, "label", ""),
    get(d, "year", ""),
    get(d, "month", ""),
    get(d, "quarter", ""),
    get(d, "sourceDataset", ""),
    get(d, "updateDate", "")
)

struct Version
    uri::String
    updateDate::String
    correctionNotice::String
    label::String
end
Version(d::Dict) = Version(
    get(d, "uri", ""),
    get(d, "updateDate", ""),
    get(d, "correctionNotice", ""),
    get(d, "label", ""),
)

struct DataSetRef
    uri::String
end
DataSetRef(d::Dict) = DataSetRef(get(d, "uri", ""))

struct Data
    years::Vector{Period}
    quarters::Vector{Period}
    months::Vector{Period}
    relatedDataset::Vector{DataSetRef}
    relatedDocuments::Vector{DataSetRef}
    relatedData::Vector{DataSetRef}
    versions::Vector{Version}
    typ::String
    uri::String
    description::MetaDataDescription
end
Data(d::Dict) = Data(
    haskey(d, "years") ? Period.(d["years"]) : [],
    haskey(d, "quarters") ? Period.(d["quarters"]) : [],
    haskey(d, "months") ? Period.(d["months"]) : [],
    haskey(d, "relatedDataset") ? DataSetRef.(d["relatedDataset"]) : [],
    haskey(d, "relatedDocuments") ? DataSetRef.(d["relatedDocuments"]) : [],
    haskey(d, "relatedData") ? DataSetRef.(d["relatedData"]) : [],
    haskey(d, "versions") ? Version.(d["versions"]) : [],
    get(d, "type", ""),
    get(d, "uri", ""),
    haskey(d, "description") ? MetaDataDescription(d["description"]) : empty(MetaDataDescription),
)

function Base.display(d::Data)
    display(d.description.title)
    if !isempty(d.years)
        println("  Annual data from $(first(d.years).year) to $(last(d.years).year)")
    end
    if !isempty(d.quarters)
        println("  Quarterly data from $(first(d.quarters).year):$(first(d.quarters).quarter) to $(last(d.quarters).year):$(last(d.quarters).quarter)")
    end
    if !isempty(d.months)
        println("  Quarterly data from $(first(d.months).year):$(first(d.months).month) to $(last(d.months).year):$(last(d.months).month)")
    end
end

function Base.getindex(d::Data, freq::Symbol)
    if freq == :yearly
        return [parse(x.value) for x in d.years]
    elseif freq == :years
        return [parse(x.year) for x in d.years]
    elseif freq == :quarterly
        return [parse(x.value) for x in d.quarters]
    elseif freq == :quarterly
        return [parse(x.value) for x in d.quarters]
    elseif freq == :monthly
        return [parse(x.value) for x in d.months]
    else
        return Float64[]
    end
end

empty(::Type{Int}) = typemin(Int)
empty(::Type{Bool}) = false
empty(::Type{String}) = ""
empty(::Type{Vector{T}}) where T = T[]

for t in (Records, MetaData, Contact, MetaDataDescription, Period, Version, Data)
    sig = :(empty(::Type{$t}))
    body = :($(Symbol(t))())
    for i = 1:nfields(t)
        fname = fieldname(t, i)
        ftype = fieldtype(t, i)
        push!(body.args, :(empty($ftype)))
    end
    @eval $sig = $body
end

function list_datasets()
    req = get("https://api.ons.gov.uk/dataset")
    d = JSON.Parser.parse(String(req.data))
    return Records(d)
end

function list_timeseries(start = 0)
    req = get("https://api.ons.gov.uk/timeseries?start=$start&limit=100")
    d = JSON.Parser.parse(String(req.data))
    return Records(d)
end

function search_datasets(q, start = 0)
    q = replace(q, " ", "+")
    req = get("https://api.ons.gov.uk/search?q=$q&limit=100&start=$(start*100)")
    d = JSON.Parser.parse(String(req.data))
    return Records(d)
end

function get_dataset(q)
    req = get("https://api.ons.gov.uk/dataset/$q/timeseries")
    d = JSON.Parser.parse(String(req.data))
    return Records(d)
end
function get_dataset(md::MetaData)
    if !isempty(md.description.datasetId)
        get_dataset(md.description.datasetId)
    else
        error("No datasetId found")
    end
end

function get_data(md::MetaData)
    dataset = md.description.datasetId
    timeseries = md.description.cdid
    req = get("https://api.ons.gov.uk/dataset/$dataset/timeseries/$timeseries/data")
    d = JSON.Parser.parse(String(req.data))
    return Data(d)
end

export search_datasets, list_datasets, list_timeseries, get_dataset, get_data
end


text:
"""
Julia translation of code from the following:
------------------------------------------------------------------------------
readable.c: My random number generator, ISAAC.
(c) Bob Jenkins, March 1996, Public Domain
You may use this code in any way you wish, and it is free.  No warrantee.
------------------------------------------------------------------------------
"""
# maximum length of message here is set to 4096
const MAXMSG = 4096

# cipher modes for encryption versus decryption modes
@enum CipherMode mEncipher mDecipher mNone

# external results
mutable struct IState
    randrsl::Array{UInt32, 1}
    randcnt::UInt32
    mm::Array{UInt32, 1}
    aa::UInt32
    bb::UInt32
    cc::UInt32
    function IState()
        this = new()
        this.randrsl = zeros(UInt32, 256)
        this.randcnt = UInt32(0)
        this.mm = zeros(UInt32, 256)
        this.aa = this.bb = this.cc = UInt32(0)
        this
    end
end

"""
    isaac
Randomize the pool
"""
function isaac(istate)
    istate.cc +=  1             # cc gets incremented once per 256 results
    istate.bb += istate.cc      # then combined with bb

    for (j, c) in enumerate(istate.mm)   # Julia NB: indexing ahead so use i for c indexing
        i = j - 1
        xmod4 = i % 4
        if(xmod4 == 0)
            istate.aa ⊻= istate.aa << 13
        elseif(xmod4 == 1)
            istate.aa ⊻= istate.aa >> 6
        elseif(xmod4 == 2)
            istate.aa ⊻= istate.aa << 2
        else
            istate.aa ⊻= istate.aa >> 16
        end
        istate.aa += istate.mm[(i + 128) % 256 + 1]
        y = istate.mm[(c >> 2) % 256 + 1] + istate.aa + istate.bb
        istate.mm[j] = y
        istate.bb = istate.mm[(y >> 10) % 256 + 1] + c
        istate.randrsl[j] = istate.bb
    end
    # not in original readable.c
    istate.randcnt = 0
end

"""
    mix
Mix the bytes in a reversible way.
"""
function mix(arr)              # Julia NB: use E for e in c code here
   (a,b,c,d,E,f,g,h) = arr
   a⊻=b<<11; d+=a; b+=c;
   b⊻=c>>2;  E+=b; c+=d;
   c⊻=d<<8;  f+=c; d+=E;
   d⊻=E>>16; g+=d; E+=f;
   E⊻=f<<10; h+=E; f+=g;
   f⊻=g>>4;  a+=f; g+=h;
   g⊻=h<<8;  b+=g; h+=a;
   h⊻=a>>9;  c+=h; a+=b;
   (a,b,c,d,E,f,g,h)
end

"""
    randinit
Make a random UInt32 array.
If flag is true, use the contents of randrsl[] to initialize mm[].
"""
function randinit(istate, flag::Bool)
    istate.aa = istate.bb = istate.cc = 0
    mixer = Array{UInt32,1}(8)
    fill!(mixer, 0x9e3779b9)             # the golden ratio
    for i in 1:4                         # scramble it
        mixer = mix(mixer)
    end
    for i in 0:8:255                     # fill in mm[] with messy stuff
        if(flag)                         # use all the information in the seed
            mixer = [mixer[j] + istate.randrsl[i+j] for j in 1:8]
        end
        mixer = mix(mixer)
        istate.mm[i+1:i+8] .= mixer
    end
    if(flag)                             # do a second pass to seed all of mm
        for i in 0:8:255
            mixer = [mixer[j] + istate.mm[i+j] for j in 1:8]
            mixer = mix(mixer)
            istate.mm[i+1:i+8] .= mixer
        end
    end
    isaac(istate)                        # fill in the first set of results
    istate.randcnt = 0
end

"""
    Get a random 32-bit value 0..MAXINT
"""
function irandom(istate)
    retval::UInt32 = istate.randrsl[istate.randcnt+1]
    istate.randcnt += 1
    if(istate.randcnt > 255)
        isaac(istate)
        istate.randcnt = 0
    end
    retval
end

"""
    Get a random character in printable ASCII range
"""
iranda(istate) = UInt8(irandom(istate) % 95 + 32)

"""
    Do XOR cipher on random stream.
    Output: UInt8 array
"""
vernam(istate, msg) = [UInt8(iranda(istate) ⊻ c) for c in msg]

"""
    Seed ISAAC with a string
"""
function iseed(istate, seed, flag)
    fill!(istate.mm, 0)
    fill!(istate.randrsl, 0)
    len = min(length(seed), length(istate.randrsl))
    istate.randrsl[1:len] .= seed[1:len]
    randinit(istate, flag)   # initialize ISAAC with seed
end

tohexstring(arr::Array{UInt8,1}) = join([hex(i, 2) for i in arr])

function test(istate, msg, key)
    # Vernam ciphertext & plaintext
    vctx = zeros(UInt8, MAXMSG)
    vptx = zeros(UInt8, MAXMSG)
    # Encrypt: Vernam XOR
    iseed(istate, Vector{UInt8}(key), true)
    vctx = vernam(istate, Vector{UInt8}(msg))
    # Decrypt: Vernam XOR
    iseed(istate, Vector{UInt8}(key), true)
    vptx = vernam(istate, vctx)
    # Program output
    println("Message: $msg")
    println("Key    : $key")
    println("XOR    : $(tohexstring(vctx))")
    # Output Vernam decrypted plaintext
    println("XOR dcr: $(join(map(c -> Char(c), vptx)))")
    0
end

"""
Test the above.
"""
const msg = "a Top Secret secret"
const key = "this is my secret key"
test(IState(), msg, key)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Squared Exponential Function with ARD

"""
    SEArd <: StationaryARD{WeightedSqEuclidean}

ARD Squared Exponential kernel (covariance)
```math
k(x,x') = σ²\\exp(- (x - x')ᵀL⁻²(x - y)/2)
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and signal standard deviation ``σ`` where
``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct SEArd <: StationaryARD{WeightedSqEuclidean}
    "Inverse squared length scale"
    iℓ2::Vector{Float64}
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        SEArd(ll::Vector{Float64}, lσ::Float64)

    Create `SEArd` with length scale `exp.(ll)` and signal standard deviation `exp(lσ)`.
    """
    SEArd(ll::Vector{Float64}, lσ::Float64) = new(exp.(-2 .* ll), exp(2 * lσ), [])
end

function set_params!(se::SEArd, hyp::VecF64)
    length(hyp) == num_params(se) || throw(ArgumentError("SEArd only has $(num_params(se)) parameters"))
    @views @. se.iℓ2 = exp(-2 * hyp[1:(end-1)])
    se.σ2 = exp(2 * hyp[end])
end

get_params(se::SEArd) = [-log.(se.iℓ2) / 2 ; log(se.σ2) / 2]
get_param_names(k::SEArd) = [get_param_names(k.iℓ2, :ll); :lσ]
num_params(se::SEArd) = length(se.iℓ2) + 1

Statistics.cov(se::SEArd, r::Float64) = se.σ2*exp(-r / 2)

@inline dk_dll(se::SEArd, r::Float64, wdiffp::Float64) = wdiffp*cov(se,r)
@inline function dKij_dθp(se::SEArd, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    if p <= dim
        return dk_dll(se, distij(metric(se),X,i,j,dim), distijk(metric(se),X,i,j,p))
    elseif p==dim+1
        return dk_dlσ(se, distij(metric(se),X,i,j,dim))
    else
        return NaN
    end
end
@inline function dKij_dθp(se::SEArd, X::MatF64, data::StationaryARDData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(se,X,i,j,p,dim)
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m8_opt_rho13_5(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m8_opt_rho13_5!(A_copy)
end

@inline function exp_sastre_m8_opt_rho13_5!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 Ba9 Bb9 B9 T2k11
    # Computing Bb2 = x*I+x*A
    coeff1=-0.11490308438229017
    coeff2=0.15730697013710307
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.11490308438229017
    coeff2=0.15730697013710307
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.004105729412589247
    coeff2=-0.006768382232040512
    coeff3=0.9999471251526917
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=-0.0028647512623116655
    coeff2=0.0009954035208072304
    coeff3=0.001211033078915847
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.027019231244596563
    coeff2=0.06170763799314927
    coeff3=-0.0018954361862412274
    coeff4=1.0067849094795742
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.013620620860785195
    coeff2=0.07638110624385147
    coeff3=0.14911649238384225
    coeff4=0.9051158545637227
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0032451095696861268
    coeff2=-0.13204423967508136
    coeff3=0.9186197117187667
    coeff4=0.42211811124544796
    coeff5=-0.01259382817947474
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.011274427044379665
    coeff2=0.001517105603937254
    coeff3=0.001158435856997462
    coeff4=-0.0047467287622741695
    coeff5=-7.137183226929846e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.011539085894536814
    coeff2=0.0935630221089867
    coeff3=0.06407916439324149
    coeff4=-0.285986613635287
    coeff5=0.044725605639523214
    coeff6=1.0008083736633122
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.02459649940876415
    coeff2=0.016824149148287185
    coeff3=0.03458106185410125
    coeff4=-0.047464943469750574
    coeff5=0.11471274073059659
    coeff6=1.0130277788963542
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=3.0826131069323353
    coeff2=0.2781411205676187
    coeff3=0.15593200742501928
    coeff4=7.925112972476668
    coeff5=1.001825018602058
    coeff6=-0.10421074361186333
    coeff7=0.16056529744528691
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.025970161938770917
    coeff2=0.05228210309031749
    coeff3=0.10608150003582509
    coeff4=0.049985371773829376
    coeff5=0.28766922644726883
    coeff6=1.936009578758007
    coeff7=1.0130307509989027
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9222548073365477
    coeff2=0.10238368416498037
    coeff3=0.1032554965257702
    coeff4=0.030058902112840666
    coeff5=-0.11384698074301966
    coeff6=0.011734571936531983
    coeff7=-0.05247830183077817
    coeff8=0.9972225209602774
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9222548073365477
    coeff2=0.10238368416498037
    coeff3=0.1032554965257702
    coeff4=0.030058902112840666
    coeff5=-0.11384698074301966
    coeff6=0.011734571936531983
    coeff7=-0.05247830183077817
    coeff8=0.9972225209602774
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing Ba9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.3035081839334926
    coeff2=0.042128860424959606
    coeff3=0.005532874783353146
    coeff4=0.05011906670785025
    coeff5=0.03702751170872807
    coeff6=0.06620554678336489
    coeff7=0.15748603557212862
    coeff8=-0.03947847759964156
    coeff9=1.058382881818128
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=-0.2738445106644693
    coeff2=-0.033002263370480726
    coeff3=0.004912635491924672
    coeff4=0.04449492500250622
    coeff5=0.03287724382394223
    coeff6=0.05879027030590609
    coeff7=0.13983408425497476
    coeff8=-0.03505336576885022
    coeff9=0.9397507089701836
    # Smart lincomb recycle B4
    memslots6 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots6,true,I*coeff1,true,true)
    # Computing B9 with operation: mult
    mul!(memslots3,memslots2,memslots6)
    # Deallocating Ba9 in slot 2
    # Deallocating Bb9 in slot 6
    # Computing T2k11 = x*I+x*A+x*B2+x*B3+x*B5+x*B6+x*B7+x*B8+x*B9
    coeff1=0.08072764899462631
    coeff2=0.02311006099515489
    coeff3=0.05527883142327172
    coeff4=-1.5378493201072977e-5
    coeff5=1.292606058383822e-5
    coeff6=6.578027159862427e-7
    coeff7=-1.3533497256574092e-7
    coeff8=1.9207658227240783e-6
    coeff9=0.9800057467529937
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots7 .+ coeff6.*memslots8 .+ coeff7.*memslots9 .+ coeff8.*memslots10 .+ coeff9.*memslots3
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    # Deallocating B9 in slot 3
    return memslots1 # Returning T2k11
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using NCDatasets
using Test

fname = tempname()

ds = NCDataset(fname,"c")
# Dimensions

ds.dim["time"] = 3

# Declare variables
ncLAT = defVar(ds,"LAT", Float64, ("time",), attrib = [
    "standard_name"             => "latitude",
    "long_name"                 => "Latitude coordiante",
    "units"                     => "degree_north",
    "ancillary_variables"       => "QC_LAT",
    "axis"                      => "Y",
])

ncQC_LAT = defVar(ds,"QC_LAT", Int8, ("time",), attrib = [
    "standard_name"             => "latitude status_flag",
    "long_name"                 => "Quality flag for latitude",
    "_FillValue"                => Int8(10),
    "flag_values"               => Int8[0, 1, 2, 3, 4, 6, 9],
    "flag_meanings"             => "no_qc_performed good_data probably_good_data probably_bad_data bad_data spike missing_value",
])

ncDEPTH = defVar(ds,"DEPTH", Float64, (), attrib = [
    "standard_name"             => "depth",
    "long_name"                 => "Depth coordinate",
    "units"                     => "m",
    "positive"                  => "down",
    "axis"                      => "Z",
    "reference_datum"           => "geographical coordinates, WGS84 projection",
])

ncLAT[:] = [1.,2.,3.]
ncQC_LAT[:] = [1,1,4]

close(ds)

ds = NCDataset(fname,"r")

@test name(NCDatasets.ancillaryvariables(ds["LAT"],"status_flag")) == "QC_LAT"
@test isequal(NCDatasets.filter(ds["LAT"],:,accepted_status_flags = ["good_data","probably_good_data"]),
    [1.,2.,missing])


text:
module REPLinModule

import Base: LineEdit, REPL

## This code was adapted from the following:
##    https://github.com/Keno/Gallium.jl/blob/1ef7111880495f3c5a7989d88a5fbc60e7eeb285/src/lldbrepl.jl
## Copyright (c) 2015: Keno Fischer. Licensed under the MIT "Expat" License:
##    https://github.com/Keno/Gallium.jl/blob/b4bc668a4cbd0f2d4f63fbdb0597a1264afd7b4d/LICENSE.md

active_module = Base

"""
    REPLinModule.set(mod)

Make `mod` the active module for the REPLinModule REPL mode. Activate this mode with the `}` key.
"""
function set(mod)
    global active_module = mod
    panel.prompt = make_prompt()
    return
end

make_prompt() = string("julia $active_module> ")

function run_repl(repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl

    main_mode = mirepl.interface.modes[1]

    const keymap = Dict{Any,Any}(
        '}' => function (s,args...)
            if isempty(s)
                if !haskey(s.mode_state,panel)
                    s.mode_state[panel] = LineEdit.init_state(repl.t,panel)
                end
                LineEdit.transition(s,panel)
            else
                LineEdit.edit_insert(s,'}')
            end
        end
    )
    # Setup the repl panel
    global panel = LineEdit.Prompt(make_prompt();
        # Copy colors from the prompt object
        prompt_prefix=Base.text_colors[:green],
        prompt_suffix=main_mode.prompt_suffix,
        on_enter = Base.REPL.return_callback)

    hp = main_mode.hist
    hp.mode_mapping[:rim] = panel
    panel.hist = hp

    panel.on_done = REPL.respond(repl,panel; pass_empty = false) do line
        if !isempty(line)
            :( eval($(REPLinModule.active_module), parse($line)) )
        else
            :(  )
        end
    end

    push!(mirepl.interface.modes, panel)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]

    panel.keymap_dict = LineEdit.keymap(Dict{Any,Any}[skeymap, mk, main_mode.keymap_dict, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults])
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, keymap)
    nothing
end

function __init__()
    if isdefined(Base, :active_repl)
        run_repl(Base.active_repl)
    else
        atreplinit() do repl
            repl.interface = Base.REPL.setup_interface(repl)
            run_repl(repl)
        end
    end
end

end # module


text:
using Pkg
Pkg.add("Optim")
using Optim

# Parameters
κ = 0.5
η = 2.0
p1 = 1.0
p2 = 1.0
w = 10.0

function fprime!(storage, c1)
    storage[1] = -1*(η/(η-1))*
    ((κ*c1[1]^((η-1)/η)+(1-κ)*((w-p1*c1[1])/p2)^((η-1)/η))^((η/(η-1))-1))*
    (((η-1)/η)*κ*(c1[1])^(((η-1)/η)-1)-((η-1)/η)*(1-κ)*((w-p1*c1[1])/p2)^(((η-1)/η)-1)*(p1/p2))
end

function gprime!(storage, c1)
    storage[1] = -1*(η/(η-1))*
    ((κ*(log(exp(c1[1])))^((η-1)/η)+(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((η-1)/η))^((η/(η-1))-1))*
    (((η-1)/η)*κ*(log(exp(c1[1])))^(((η-1)/η)-1)-((η-1)/η)*(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(((η-1)/η)-1)*(p1/p2))
end

function g2prime!(storage, c1)
    storage[1] = -1*(η/(η-1))*((η/(η-1))-1)*((κ*log(exp(c1[1]))^((η-1)/η))+(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((η-1)/η))^((η/(η-1))-2)*
    (((η-1)/η)*κ*log(exp(c1[1]))^(((η-1)/η)-1)-((η-1)/η)*(1-κ)*(p1/p2)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(((η-1)/η)-1))*
    (((η-1)/η)*κ*log((exp(c1[1])))^(((η-1)/η)-1)-((η-1)/η)*(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(((η-1)/η)-1)*(p1/p2))-
    (((η-1)/η)*(((η-1)/η)-1)*κ*log(exp(c1[1]))^(((η-1)/η)-2)+((η-1)/η)*(1-κ)*(((η-1)/η)-1)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(((η-1)/η)-2)*(p1/p2)^2)*
    ((κ*log(exp(c1[1]))^((η-1)/η)+(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((η-1)/η))^((η/(η-1))-1))

end

function utility_maximizer(κ, η, p1, p2, w, initial_guess, method)
    if method == "Conjugate-Gradient"
        f(c1) = -1*(κ*c1[1]^((η-1)/η)+(1-κ)*((w-p1*c1[1])/p2)^((η-1)/η))^(η/(η-1))
        lower = [0.000001]
        upper = [w/p1]
        result = Optim.minimizer(optimize(f, fprime!, lower, upper, initial_guess, Fminbox(ConjugateGradient())))[1]
        return result
    elseif method == "Newton"
        g(c1) = -1*(κ*log(exp(c1[1]))^((η-1)/η)+(1-κ)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((η-1)/η))^(η/(η-1))
        result = Optim.minimizer(optimize(g, gprime!, g2prime!, initial_guess, Newton()))[1]
        return result
    end
end

initial_guess = [2.0]
mtd = "Newton"
c1_star = utility_maximizer(κ, η, p1, p2, w, initial_guess, mtd)
c2_star = (w-p1*c1_star)/p2


text:
using Pidoh
using Random

n = 20
problem = OneMax(n)

@testset "Simulated Annealing Correctness" begin
    algorithm = SimulatedAnnealing(FixedCooling(50.0), stop = FixedBudget(10^10))
    instance = Instance(generator = RandBitStringIP(n), problem)
    experiment = Experiment("_tests_sac", [algorithm], [instance], save = false)
    res = run(experiment)
    optx = res.jobs[1].trace.optimum
    @test sum(optx.individual) == optimum(problem)
end

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
workspace()
using Base.Test
using FactCheck
using JSON
using nmr

#addprocs(Sys.CPU_CORES)

@time cd(joinpath(Pkg.dir("nmr"), "test")) do
  function run_test(testfile)
    info("Running test file $(testfile)")
    include(testfile)
  end
  testfiles = [f for f in readdir(".") if isfile(f) && startswith(f, "test_") && endswith(f, ".jl")]
  map(run_test, testfiles)
end

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
## ACHI #####
"""
    ACHI

An ACHI (Australian Classification of Health Interventions) code.
"""
struct ACHI
  data::NTuple{7,UInt8}
end

function validachiinput(str)
  achifmt = r"^[[:digit:]]{5}(-){0,1}[[:digit:]]{2}$"
  return occursin(achifmt, str)
end

"""
    ACHI(str::String, validateinput=false)

Create an ACHI code from a string.

ACHI codes consist of 7 digits (0-9), usually with a hypen after the first 5 digits.

If `validateinput=true` then the input is checked using a regex. Note this only checks the input format; it does not
check if the code actually exists in ACHI.
"""
function ACHI(str::String, validateinput = false)
  punct = occursin("-", str) ? true : false
  if validateinput
    validachiinput(str) || throw(
      DomainError(
        str,
        "ACHI codes should have format `NNNNN-NN` where `N` is a decimal digit. The hyphen can be omitted.",
      ),
    )
  end
  ch6 = punct ? 7 : 6
  return ACHI((NTuple{5,UInt8}.(str[1:5])..., NTuple{2,UInt8}.(str[ch6:end])...))
end

## ACHI Base functions #####

function Base.print(io::IO, achi::ACHI, punct = ICDOPTS[:punct])
  if punct
    write(io, Ref(achi.data[1:5]), "-", Ref(achi.data[6:7]))
  else
    write(io, Ref(achi))
  end
end

Base.show(io::IO, ::MIME"text/plain", achi::ACHI) = print(io, achi, ICDOPTS[:punct])

Base.convert(::Type{T}, x::AbstractString) where {T<:ACHI} = T(x)

## Supporting types #####
"""
    ICD10AMAge

Age group encoding used in the ICD-10-AM electronic code lists.

  - 0-6 : 0 to 6 days
  - 11-13 : 1 to 3 weeks
  - 101-111 : 1 to 11 months
  - 201-324 : 1 to 124 years

All other values are invalid.
"""
struct ICD10AMAge
  agecode::Int16

  function ICD10AMAge(age)
    (age in 0:6 || age in 11:13 || age in 101:111 || age in 201:324) ? new(age) :
    throw(DomainError(age, "Invalid age code"))
  end
end

"""
    ICD10AMAge(agecode::String)

Convert a string representation of the age code to an ICD10AMAge.
"""
ICD10AMAge(age::String) = ICD10AMAge(parse(Int, age))

function Base.show(io::IO, age::ICD10AMAge)
  if age.agecode in 0:6
    print(io, string(age.agecode) * " days")
  elseif age.agecode in 11:13
    print(io, string(age.agecode - 10) * " weeks")
  elseif age.agecode in 101:111
    print(io, string(age.agecode - 100) * " months")
  else
    print(io, string(age.agecode - 200) * " years")
  end
end

## ACHI functions #####
"""
    isvalidcode(code:ACHI)

Takes an ACHI code and checks it is in the electronic code lists and valid.

Electronic code lists can be obtained from the [Independent Hospital Pricing
Authority](https://www.ihpa.gov.au/what-we-do/icd-10-am-achi-acs-classification),
and need to be imported into the local package before this function
can be used. See [`importachicodes`](@ref)

"Valid" codes are those that exist.
"""
function isvalidcode(achi::ACHI)
  isdefined(ICD10Utilities, :ACHIcodes) ||
    isfile(normpath(@__DIR__, "..", "data", "achicodes.jld2")) ||
    throw(
      ErrorException(
        "You need to import the ICD10AM/ACHI electronic code lists (obtained from Independent Hospital Pricing Authority at https://ar-drg.laneprint.com.au/) to use validation. Run `importachicodes(<filename>` where `<filename>` is the path to the `interven.txt` file from IHPA.",
      ),
    )
  if !isdefined(ICD10Utilities, :ACHIcodes)
    global ACHIcodes = load(normpath(@__DIR__, "..", "data", "achicodes.jld2"), "achicodes")
  end

  f(achi, codelist = ACHIcodes)::Bool = achi in codelist.achicode

  return f(achi)
end
## TO DO #####
# Need a method that takes a date and checks if the code is inactive at that date.
# Also want to check if a code is valid in a particular version


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 54640fbe-86ce-42bf-89cd-6085fa2052cb
md"""
# Data collection and usefulness

The way data is collected is immensely impactful to how useful the data is subsequently.

First, we need to differentiate between computerized data and computable data.

## PDF Files

### Scanned Documents

A scanned document stored as a PDF is computerized but it is not computable. From the computer's perspective it is simply a blob of pixels or more precisely an array of numbers between 0 and 255 that describe the brightness on the screen that we as humans read as letters, horizontal and vertical lines, letterhead, etc. As an example, consider this scanned document

![chemo worksheet](https://github.com/chapmanbe/raw-stuff/blob/main/DataSources/chemo_worksheet.png?raw=true)

I could potentially run an optical character recognition program to identiy text, but this is going to be a noisy (error prone) process that does not preserve an structure/context that may be implied in the original document. For example, running the OCR program [tesseract](https://github.com/tesseract-ocr/tesseract) on this document and then trying to copy the identified text yielded "mrieTorich" for "Antibiotics" on line 12. Recognizing the handwriting was hopeless.

### PDF generated by word processing software

## Text entry

Here is an excerpt from an actual note in my medical record:

```text
 PROBLEM:   HYPERTENSION

     1. CLINIC BP TODAY = 122/79  WITH GOAL BELOW 130/80

     2. CURRENT MEDICATION:  LISINOPRIL 10 MG DAILY

PROBLEM:  Chronic kidney disease DUE TO REMOVAL OF KIDNEY

     1. LAST CREATININE = 1.52 12/2018
           --HAVE BEEN STABLE FOR YEARS
         NORMAL MICROALBUMIN/CREATININE RATIO AND URINALYSIS
```

From a computational perspective, this string is an advance over the PDF-as-image data described previously. The computer recognizes `B` as a _B_ and `P` as a _P_ yet has no concept that `BP` represents blood pressure or that `122/79` represents the dual concepts

- _systolic blood pressure: value=122, units of measurement=mmHg_
- _diastolic blood pressure: value=79, units of measurement=mmHg_

In order to make these values computable, a natural language processing program would have to be designed and implemented to recognize this, a process that would be imperfect at best, most likely imperfect and time consuming.

By contrast, below is an excerpt from the HL7 document (converted from XML to JSON) containing the structured data entry for the blood pressure measurement taken during this visit. While hard to interpret by eye without knowing all the OID values used for HL7 documents, we can make the following observations about the observation:

1. The coding system that the data uses is clearly identified: (LOINC)
1. The concept is provided with its non-semantic ID (_"8480-6"_) as well as its text representation _"Systolic blood pressure"_
1. _"122"_ is identified as the value
1. _"mm[Hg]"_ is identified as the unit
1. Similar observations for the diastolic value

__These data are computable!__

```JSON
"component": [
                    {
                      "observation": {
                        "@classCode": "OBS",
                        "@moodCode": "EVN",
                        "templateId": [
                          {
                            "@root": "2.16.840.1.113883.10.20.22.4.27"
                          },
                          {
                            "@root": "2.16.840.1.113883.10.20.22.4.27",
                            "@extension": "2014-06-09"
                          }
                        ],
                        "id": {
                          "@root": "1.2.840.114350.1.13.90.2.7.1.2109.1",
                          "@extension": "5623842480-sysBP-22122"
                        },
                        "code": {
                          "@code": "8480-6",
                          "@codeSystem": "2.16.840.1.113883.6.1",
                          "@codeSystemName": "LOINC",
                          "originalText": "Systolic blood pressure"
                        },
                        "text": {
                          "reference": {
                            "@value": "#sysBP_5623842480"
                          }
                        },
                        "statusCode": {
                          "@code": "completed"
                        },
                        "effectiveTime": {
                          "@value": "20190318182800+0000"
                        },
                        "value": {
                          "@xsi:type": "PQ",
                          "@unit": "mm[Hg]",
                          "@value": "122",
                          "@xmlns:xsi": "http://www.w3.org/2001/XMLSchema-instance"
                        }
                      }
                    },
                    {
                      "observation": {
                        "@classCode": "OBS",
                        "@moodCode": "EVN",
                        "templateId": [
                          {
                            "@root": "2.16.840.1.113883.10.20.22.4.27"
                          },
                          {
                            "@root": "2.16.840.1.113883.10.20.22.4.27",
                            "@extension": "2014-06-09"
                          }
                        ],
                        "id": {
                          "@root": "1.2.840.114350.1.13.90.2.7.1.2109.1",
                          "@extension": "5623842480-diaBP-22122"
                        },
                        "code": {
                          "@code": "8462-4",
                          "@codeSystem": "2.16.840.1.113883.6.1",
                          "@codeSystemName": "LOINC",
                          "originalText": "Diastolic blood pressure"
                        },
                        "text": {
                          "reference": {
                            "@value": "#diaBP_5623842480"
                          }
                        },
                        "statusCode": {
                          "@code": "completed"
                        },
                        "effectiveTime": {
                          "@value": "20190318182800+0000"
                        },
                        "value": {
                          "@xsi:type": "PQ",
                          "@unit": "mm[Hg]",
                          "@value": "79",
                          "@xmlns:xsi": "http://www.w3.org/2001/XMLSchema-instance"
                        }
                      }
                    }
```

One thing that might strike you about this data is how verbose it is. To express "79 mm[Hg]" (8 characters). The number of characters used to represent this in JSON is 132 characters. This inefficiency is an important tool to reduce ambiguity.

> Redundancy—inefficient by definition—serves as the antidote to confusion. It provides second chances. Every natural language has redundancy built in; this is why people can understand text riddled with errors and why they can understand conversation in a noisy room. Most of the time, redundancy in language is just part of the background....Another specialized language provides a perfect analog: the language of aviation radio. Numbers and letters make up much of the information passed between pilots and air traffic controllers: altitudes, vectors, aircraft tail numbers, runway and taxiway identifiers, radio frequencies. This is critical communication over a notoriously noisy channel, so a specialized alphabet is employed to minimize ambiguity. The spoken letters B and V are easy to confuse; bravo and victor are safer. M and N become mike and november. In the case of numbers, five and nine , particularly prone to confusion, are spoken as *fife* and *niner*. (Gleick, James, _The Information: A History, a Theory, a Flood_)

### "Structured" Free Text Entry

The text excerpt I provided above was a small section of the SOAP note generated by my GP after the visit. Although there is structure provided by the SOAP format, it is essentially a free flowing stream of text.

Other text entry occurs in more constrained fields. For example, there might be a field for "First/Given Name" and "Surname/Family Name/Last Name". With the plethera of names in the world, this could not feasibly be entered anyway other than typing. This mgives all sorts of opportunities for typographical errors. I think that I'm `Brain` rather than `Brian` in a significant number of databases. The longer the text field, the more opportunity for typos, synonyms, ambiguous abbreviations. etc.

"""

# ╔═╡ 568ad495-ac20-40e3-9abb-220e5a0e1c78
md"""

## Problems with structured entry

Using standards for data greatly expand the usefulness of the resulting information. However, just because standards are used does not mean the standards are correct. As mentioned earlier with ICD codes, the purpose and manner in which codes are used  will introduce error and or bias into the resulting information.

Further, codes are still liable to data entry errors.

Jo-Anne Manski-Nankervis, a GP participating in this course, shared the following examples of coded data entry errors:

>One pathology company labelling both fasting and random glucose pathology tests results with the same code: 14749-6.
>
>The other example is a large general practice data set which included an ATC code for a vaccine that had been applied to other medications at random. Relying on the codes would have resulted in a underestimation of the prescription of cardiovascular medications in general practice.

### Drop down menus (and other prepopulated selection forms) for data entry

Consider this story that ran on ABC News 14 September 2021:

- [Indigenous vaccination rates in Victoria revised downwards after software error detected](https://www.abc.net.au/news/2021-09-14/indigenous-vaccination-rates-in-victoria-dashed-software-error/100459000)

Default values are problematic as they don't require conscious assertion of a fact by the person entering the data. For example, should we prepopulate all blood pressure measurements with 120 mmHg for systolic and 80 mmHg for diastolic?

Other problems with drop downs include:

- Simply selecting the first element
- Proximity errors where an adjacent menu item is selected instead of the intended/correct item.

### How can structured entry help ensure data quality?

On the positive side, coded data entry allows for checking for data integrity at the point of data entry.

- Data validation:
   - Data entry can be designed such that range checking is enforced. For example, weight must be a positive quantity or blood pressure be within physiologically possible values, with conscious overrides possible
   - Format checking:
      - Ensure that numbers are numeric
      - Phone numbers, email, etc. are in appropriate format

However, coded data entry is more tedious than simply typing free text and so a data collection strategy needs to balance ease of use with accuracy and completeness.
"""

# ╔═╡ Cell order:
# ╟─54640fbe-86ce-42bf-89cd-6085fa2052cb
# ╠═568ad495-ac20-40e3-9abb-220e5a0e1c78


text:
function submanifold_mapping(vs)
    u1 = vs[2,:] - vs[1,:]
    u2 = vs[3,:] - vs[1,:]
    u3 = cross(u1, u2)
    T = eye(4)
    T[4, 1:3] = - vs[1,:]
    M = eye(4)
    M[1:3, 1:3] = [u1 u2 u3]
    return T*M
end

function spatial_index(V::LinearAlgebraicRepresentation.Points, EV::LinearAlgebraicRepresentation.ChainOp, FE::LinearAlgebraicRepresentation.ChainOp)
    d = 3
    faces_num = size(FE, 1)
    IntervalsType = IntervalValue{Float64, Int64}
    boxes1D = Array{IntervalsType, 2}(0, d)
    for fi in 1:faces_num
        vidxs = (abs.(FE[fi:fi,:])*abs.(EV))[1,:].nzind
        intervals = map((l,u)->IntervalsType(l,u,fi), LinearAlgebraicRepresentation.bbox(V[vidxs, :])...)
        boxes1D = vcat(boxes1D, intervals)
    end
    trees = mapslices(IntervalTree{Float64, IntervalsType}, sort(boxes1D, 1), 1)

    function intersect_intervals(intervals)
        cells = Array{Int64,1}[]
        for axis in 1:d
            vs = map(i->i.value, intersect(trees[axis], intervals[axis]))
            push!(cells, vs)
        end
        mapreduce(x->x, intersect, cells)
    end

    mapping = Array{Int64,1}[]
    for fi in 1:faces_num
        cell_indexes = setdiff(intersect_intervals(boxes1D[fi, :]), [fi])
        push!(mapping, cell_indexes)
    end

    mapping
end

function face_int(V::LinearAlgebraicRepresentation.Points, EV::LinearAlgebraicRepresentation.ChainOp, face::LinearAlgebraicRepresentation.Cell)

    vs = LinearAlgebraicRepresentation.buildFV(EV, face)
    retV = LinearAlgebraicRepresentation.Points(0, 3)

    visited_verts = []
    for i in 1:length(vs)
        o = V[vs[i],:]
        j = i < length(vs) ? i+1 : 1
        d = V[vs[j],:] - o

        err = 10e-8
        if !(-err < d[3] < err)

            alpha = -o[3] / d[3]

            if -err <= alpha <= 1+err
                p = o + alpha*d

                if -err < alpha < err || 1-err < alpha < 1+err
                    if !(LinearAlgebraicRepresentation.vin(p, visited_verts))
                        push!(visited_verts, p)
                        retV = [retV; reshape(p, 1, 3)]
                    end
                else
                    retV = [retV; reshape(p, 1, 3)]
                end
            end
        end

    end

    vnum = size(retV, 1)

    if vnum == 1
        vnum = 0
        retV = LinearAlgebraicRepresentation.Points(0, 3)
    end
    enum = Int(vnum / 2)
    retEV = spzeros(Int8, enum, vnum)

    for i in 1:enum
        retEV[i, 2*i-1:2*i] = [-1, 1]
    end

    retV, retEV
end


text:
const fenjuPos =2:2
const guanlichuPos=2:10
const signalChnPos=11:11
const SignalIDPos=12:14
const SignalIDfullPos=11:14

mutable struct CsvInfos
    csvInfoSignalID2FiberPos
    csvInfoDefencePara
    csvInfoMachineInfo
    防区表fieldname
    防区表signalStringName
    防区表keyFieldNames
    防区表filename

    老版参数表fieldname
    老版参数表signalStringName
    老版参数表keyFieldNames
    老版参数表filename

    主机信息表fieldname
    主机信息表signalStringName
    主机信息表keyFieldNames
    主机信息表filename
    ParasPosNames
    function  CsvInfos()
         csvInfoSignalID2FiberPos = nothing
         csvInfoDefencePara = nothing
         csvInfoMachineInfo = nothing
         new(csvInfoSignalID2FiberPos,   csvInfoDefencePara, csvInfoMachineInfo,"","","","","","","","","","","","","")
    end
end

function findMachineRegID( ;csvInfoMachineInfo::CSVinfo ,machine友好名::String="xxxx001")::String
    if occursin( "友好名" , csvInfoMachineInfo.keyFieldNames)
        fields =split( csvInfoMachineInfo.keyFieldNames,",")
        pos友好名 =0
        pos主机编号=0
        for i = 1: length(fields )
            if fields[i] =="友好名"
                pos友好名 = i
            elseif  fields[i] =="主机编号"
                pos主机编号=i

            end
        end
        for infos in  csvInfoMachineInfo.Area2FiberLen
            if infos[pos友好名] == machine友好名
                return infos[pos主机编号]
            end
        end

    else
        # return ""
    end
    return ""
end

function getSpecMachineSlogansInfo( ; csvInfoSignalID2FiberPos1::CSVinfo ,machineRegID::String="A302042001",通道1::String="1")
    a=[]
    refobj = csvInfoSignalID2FiberPos1.Area2FiberLen
    fields =split( csvInfoSignalID2FiberPos1.keyFieldNames,",")
    pos定位界标1 = 0
    pos定位界标1光程 = 0
    pos定位界标2 = 0
    pos定位界标2光程 = 0
    for i = 1: length(fields )
        if fields[i] == "定位界标1"
            pos定位界标1 = i
        elseif  fields[i] == "定位界标1光程"
            pos定位界标1光程 = i
        elseif  fields[i] == "定位界标2光程"
            pos定位界标2光程 = i
        elseif fields[i] == "定位界标2"
            pos定位界标2 = i
        end
    end
    # "定位界标1,定位界标1光程"
    rowobj = nothing
    lastSignalID = nothing
    lastFIberPos = 0
    for i = 1:length(refobj)
        rowobj = refobj[i]
        jiebiao1 = rowobj[pos定位界标1]
        # println( jiebiao1[ guanlichuPos] ,"-",machineRegID[guanlichuPos] )
        if jiebiao1[ guanlichuPos] ==  machineRegID[guanlichuPos]  && jiebiao1[signalChnPos] == 通道1
                append!(a,[ [jiebiao1 ,parse( UInt,rowobj[pos定位界标1光程] )  ]  ] )
                lastSignalID = rowobj[pos定位界标2]
                lastFIberPos = rowobj[pos定位界标2光程]
        end
    end
    if lastSignalID  != nothing
        append!(a,[ [lastSignalID ,parse( UInt,lastFIberPos ) ]  ] )
    end

    return a
end

function getPos_Signal(位置::String,Rows=[])::String
    pos1=parse(UInt,位置)
    # println( pos1 )
    pos0 =0
    pos2 = 0
    for idx =1: length(Rows)-1
     row1 = Rows[idx]
     row2 = Rows[idx+1]
        pos0 = row1[2]
        pos2 = row2[2]
        if pos1 >= pos0 && pos1 <= pos2
            # println( row[1],"-",row[2],"-",row[3],"-",row[4] )
            if pos2 - pos1 > pos1 -pos0
                return string(row1[1][SignalIDfullPos],"+",pos1 -pos0 )  ;
                # println( 位置,":",row[1],"+",pos1 -pos0,", ",row[2], )
            else
                return string(row2[1][SignalIDfullPos],"-",pos2 -pos1 )  ;
                # println( 位置,":",row[3],"-",pos2 -pos1,", ",row[4], )
            end
        end
        if pos1 < pos0
            return string(row1[1][SignalIDfullPos],"-",pos0 -pos1)  ;
        end
    end
    row = Rows[end]
    pos2 = row[ 2 ]
    if pos1 > pos2
        return string(row[1][SignalIDfullPos],"+",pos1 -pos2)  ;
    end
    # return ""
end

function getPos_FiberLen(位置::String,Rows=[])::Int
    sign = -1
    base = nothing
    biase = nothing
    if occursin("-" ,位置)
        sign = -1
        base,biase = split(位置,"-")
    elseif occursin("+" ,位置)
        sign = 1
        base,biase = split(位置,"+")
    else
        return ""
    end

    if base == "0" && biase == "0"
        return 0
    end
    signalID= nothing
    fiberPos = -1
    for idx =1: length(Rows)
     row1 = Rows[idx]
     signalID = row1[1][SignalIDfullPos]
     if base == signalID
        fiberPos = row1[2]
        break
        end
    end

    if fiberPos == -1
        println("not find...$位置")
        return -1
    else
        # println("$fiberPos,$sign,$biase")
        return fiberPos + sign*parse(UInt,biase)
    end
end

function getSlogenCorordinats( ;通道::String,位置::String,csvInfoSignalID2FiberPos2::CSVinfo,machineRegID1::String="A302042001")::String
    Rows= getSpecMachineSlogansInfo( csvInfoSignalID2FiberPos1 = csvInfoSignalID2FiberPos2  ,
        machineRegID = machineRegID1,
        通道1 = 通道)
    posStr = getPos_Signal(位置,Rows)
    return posStr
end

function convertPara2SignalCordinate( CsvInfoObj::CsvInfos;machine友好名1="xxxx001")
    machineRegID1 = findMachineRegID( csvInfoMachineInfo=CsvInfoObj.csvInfoMachineInfo ,machine友好名=machine友好名1)
    sort!(CsvInfoObj.csvInfoDefencePara.Area2FiberLen)
    keyname = "通道,起始位置(米),结束位置(米)"

    keynames = split( keyname ,",")
    fields = split( CsvInfoObj.csvInfoDefencePara.keyFieldNames,",")
    pos通道 = 0
    pos起始位置 = 0
    pos结束位置 = 0

    for i = 1: length(fields )
        if fields[i] == keynames[1]
              pos通道 = i
        elseif  fields[i] ==  keynames[2]
              pos起始位置 = i
        elseif  fields[i] == keynames[3]
              pos结束位置 = i
        end
    end
    rows = nothing
    tongdaoID =""
    dataobjs = CsvInfoObj.csvInfoDefencePara.Area2FiberLen
    datalen = length(dataobjs  )
    for idx =1:datalen
         dataobj = dataobjs[ idx ]
         通道 = dataobj[pos通道]
         起始位置 = dataobj[pos起始位置]
         结束位置 = dataobj[pos结束位置]
         if 通道 != tongdaoID
              rows = getSpecMachineSlogansInfo( csvInfoSignalID2FiberPos1= CsvInfoObj.csvInfoSignalID2FiberPos  ,
              machineRegID=machineRegID1,
              通道1=string(通道) )
              tongdaoID= string(通道)
         end

        #  print(通道 ," ",起始位置," ",结束位置,":" )

         posStr = getPos_Signal(string(起始位置),rows)
        #  print( "  ",起始位置,":",posStr )
         dataobj[pos起始位置]= posStr
         posStr = getPos_Signal(string(结束位置),rows)
        #  println( "  ",结束位置,":",posStr )
         dataobj[pos结束位置] = posStr
        #  println( dataobj )
    end
end

function convertPara_Fiberlength( CsvInfoObj::CsvInfos;machine友好名1="xxxx001")
    machineRegID1 = findMachineRegID( csvInfoMachineInfo=CsvInfoObj.csvInfoMachineInfo ,machine友好名=machine友好名1)
    sort!(CsvInfoObj.csvInfoDefencePara.Area2FiberLen)

    keyname =  CsvInfoObj.ParasPosNames #"通道,起始位置(米),结束位置(米)"

    keynames = split( keyname ,",")
    fields = split( CsvInfoObj.csvInfoDefencePara.keyFieldNames,",")
    pos通道 = 0
    pos起始位置 = 0
    pos结束位置 = 0

    for i = 1: length(fields )
        if fields[i] == keynames[1]
              pos通道 = i
        elseif  fields[i] ==  keynames[2]
              pos起始位置 = i
        elseif  fields[i] == keynames[3]
              pos结束位置 = i
        end
    end
    rows = nothing
    tongdaoID =""
    dataobjs = CsvInfoObj.csvInfoDefencePara.Area2FiberLen
    datalen = length(dataobjs  )
    for idx =1:datalen
         dataobj = dataobjs[ idx ]
         通道 = dataobj[pos通道]
         起始位置 = dataobj[pos起始位置]
         结束位置 = dataobj[pos结束位置]
         if 通道 != tongdaoID
              rows = getSpecMachineSlogansInfo( csvInfoSignalID2FiberPos1 = CsvInfoObj.csvInfoSignalID2FiberPos  ,
              machineRegID=machineRegID1,
              通道1=string(通道) )
              tongdaoID= string(通道)
         end

        #print(通道 ," === :",起始位置," ",结束位置,":" )

         pos1 = getPos_FiberLen(string(起始位置),rows)
       #  print( "---- ",起始位置,"=>",pos1 )
         dataobj[pos起始位置]= pos1
         pos2 = getPos_FiberLen(string(结束位置),rows)
     #   println( "----  ",结束位置,"=>",pos2 )
         dataobj[pos结束位置] = pos2
        #  println( dataobj )
    end
end

function testCode()
    machineRegID1 = findMachineRegID( csvInfoMachineInfo=csvInfoMachineInfo ,machine友好名="xxxx001")
    通道="2"
    Rows = getSpecMachineSlogansInfo( csvInfoSignalID2FiberPos1=csvInfoSignalID2FiberPos  ,
         machineRegID=machineRegID1,
         通道1=通道)
    # println(machineRegID1)
    位置1="1230"
    posStr = getPos_Signal(位置1,Rows)
    println( 位置1,":",posStr )

    位置1="12310"
    posStr = getPos_Signal(位置1,Rows)
    println( 位置1,":",posStr )

    位置1="100"
    posStr = getPos_Signal(位置1,Rows)
    println( 位置1,":",posStr )
    位置1="24950"
    posStr = getPos_Signal(位置1,Rows)
    println( 位置1,":",posStr )

    位置1="24950"
    posStr = getPos_Signal(位置1,Rows)
    println( 位置1,":",posStr )

    # for row in Rows
    #      println( row[1],"-",row[2])
    # end
end

function basicInit( CsvInfoObj::CsvInfos)
    CsvInfoObj.csvInfoSignalID2FiberPos = CSVinfo(
         fieldname =CsvInfoObj.防区表fieldname,
         signalStringName=CsvInfoObj.防区表signalStringName,
         keyFieldNames =CsvInfoObj.防区表keyFieldNames
    )
    readinCSVFile!(CsvInfoObj.csvInfoSignalID2FiberPos,csvfile=CsvInfoObj.防区表filename )
    sort!(CsvInfoObj.csvInfoSignalID2FiberPos.Area2FiberLen)
    # show( csvInfoSignalID2FiberPos )

    CsvInfoObj.csvInfoMachineInfo = CSVinfo(
         fieldname = CsvInfoObj.主机信息表fieldname,
         signalStringName = CsvInfoObj.主机信息表signalStringName,
         keyFieldNames = CsvInfoObj.主机信息表keyFieldNames
    )
    readinCSVFile!(CsvInfoObj.csvInfoMachineInfo,csvfile = CsvInfoObj.主机信息表filename)
    # show(csvInfoMachineInfo )

    if CsvInfoObj.老版参数表filename !=""
        CsvInfoObj.csvInfoDefencePara= CSVinfo(
             fieldname =CsvInfoObj.老版参数表fieldname,
             signalStringName =CsvInfoObj.老版参数表signalStringName,
             keyFieldNames =CsvInfoObj.老版参数表keyFieldNames
        )
        readinCSVFile!(CsvInfoObj.csvInfoDefencePara,csvfile=CsvInfoObj.老版参数表filename)
        sort!(CsvInfoObj.csvInfoDefencePara.Area2FiberLen)
        # show(csvInfoDefencePara )
    end

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using Oceananigans
using Oceananigans.Grids: min_Δx, min_Δy, min_Δz
using JLD2
using OffsetArrays
using BenchmarkTools
using LinearAlgebra
using Plots
using Test
using Adapt

"""
This simulation is a simple 1D advection to test the
validity of the stretched WENO scheme
"""

function multiple_steps!(model)
    for i = 1:1000
        time_step!(model, 1e-6)
    end
    return nothing
end

#parameters
N    = 40
arch = CPU()
iter = 1000

# regular "stretched" grid
Freg = range(0, 1, length = N+1)

# alternating grid
Fsaw(j) = 1 / N  * (j - 1) + 0.1 * 1 / N * mod(j - 1, 2)

# center-coarsened grid
Δstr2(i, N) = i < N/4 ? 1 : ( i > N*0.75 ? 1 : ( i < N/2 ? 1.2 * (i - N/4) + 1 : 1.2 * (3*N/4 - i) + 1  )  )
Fstr2 = zeros(Float64, N+1)
for i = 2:N+1
     Fstr2[i] = Fstr2[i-1] + Δstr2(i-1, N)
end
Fstr2 ./= Fstr2[end]

# solutions
solution  = Dict()
real_sol  = Dict()
coord     = Dict()
time      = Dict()

# 1D grid constructions
grid_reg  = RectilinearGrid(arch, size = N, x = Freg,  halo = 3, topology = (Periodic, Flat, Flat))
grid_str  = RectilinearGrid(arch, size = N, x = Fsaw,  halo = 3, topology = (Periodic, Flat, Flat))
grid_str2 = RectilinearGrid(arch, size = N, x = Fstr2, halo = 3, topology = (Periodic, Flat, Flat))

# placeholder for the four different advection schemes
#  (1) WENO5(),
#  (2) WENO5(grid=grid),
#  (3) WENO5(grid=grid, stretched_smoothness=true),
#  (4) WENO5(grid=grid, stretched_smoothness=true, zweno=true)

advection = [WENO5(), WENO5(), WENO5(), WENO5()]
schemes   = [:wreg, :wstr, :wstrS, :wstrZ]

# mask for the initial condition
mask(y) = (y < 0.6 && y > 0.4) ? 1 : 0
c₀_1D(x, y, z) = @. 10 * exp(-(x - 0.5)^2 / 0.2^2)
c₀_2D(x, y, z) = mask(x) * mask(y)

# Checking the accuracy of different schemes with different settings
for (gr, grid) in enumerate([grid_reg, grid_str, grid_str2])

    U = Field{Face, Center, Center}(grid)
    parent(U) .= 1

    Δt_max   = 0.2 * min_Δx(grid)
    end_time = iter * Δt_max

    for (adv, scheme) in enumerate(advection)
        if adv == 2
            scheme = WENO5(grid = grid)
        end
        if adv == 3
            scheme = WENO5(grid = grid, stretched_smoothness = true)
        end
        if adv == 4
            scheme = WENO5(grid = grid, stretched_smoothness = true, zweno = true)
        end

        model = HydrostaticFreeSurfaceModel(architecture = arch,
                                                    grid = grid,
                                                 tracers = :c,
                                        tracer_advection = scheme,
                                              velocities = PrescribedVelocityFields(u=U),
                                                coriolis = nothing,
                                                 closure = nothing,
                                                buoyancy = nothing)

        set!(model, c=c₀_1D)
        c = model.tracers.c

        simulation = Simulation(model,
                                Δt = Δt_max,
                                stop_time = end_time)

        for i = 1:end_time/Δt_max
            csim                             = adapt(CPU(), c)
            solution[(schemes[adv], Int(i))] = csim[1:N, 1, 1]

            for j = 1:10
                time_step!(model, Δt_max)
            end
        end

        time[(schemes[adv], gr)] = @belapsed multiple_steps!($model)
    end

    x = adapt(CPU(), grid.xᶜᵃᵃ)
    x = x[1:N]

    anim = @animate for i ∈ 1:end_time/Δt_max
        plot(x, c₀_1D(x, 1, 1), seriestype=:scatter, legend = false, title ="red: U, blue: S, green: β, black: Z")
        plot!(x, solution[(schemes[1], Int(i))], linewidth = 4, linecolor =:red  , legend = false)
        plot!(x, solution[(schemes[2], Int(i))], linewidth = 3, linecolor =:blue , legend = false)
        plot!(x, solution[(schemes[3], Int(i))], linewidth = 2, linecolor =:green, legend = false)
        plot!(x, solution[(schemes[4], Int(i))], linewidth = 1, linecolor =:black, legend = false)
    end
    mp4(anim, "anim_1D_$(gr).mp4", fps = 15)

end

# """
# 2D advection test
# """

solution2D  = Dict()

grid_reg  = RectilinearGrid(size = (N, N), x = Freg,  y = Freg,  halo = (3, 3), topology = (Periodic, Periodic, Flat), architecture = arch)
grid_str  = RectilinearGrid(size = (N, N), x = Fsaw,  y = Fsaw,  halo = (3, 3), topology = (Periodic, Periodic, Flat), architecture = arch)
grid_str2 = RectilinearGrid(size = (N, N), x = Fstr2, y = Fstr2, halo = (3, 3), topology = (Periodic, Periodic, Flat), architecture = arch)

for (gr, grid) in enumerate([grid_reg, grid_str, grid_str2])

    U = Field{Face, Center, Center}(grid)
    V = Field{Center, Face, Center}(grid)

    parent(U) .= 1
    parent(V) .= 0.3

    Δt_max   = 0.1 * min_Δx(grid)  # faster in U than in V
    end_time = 2 * iter * Δt_max

    maxAdv1 = []; maxAdv2 = []; maxAdv3 = []; maxAdv4 = [];
    for (adv, scheme) in enumerate(advection)

        if adv == 2
            scheme = WENO5(grid = grid)
        end
        if adv == 3
            scheme = WENO5(grid = grid, stretched_smoothness = true)
        end
        if adv == 4
            scheme = WENO5(grid = grid, stretched_smoothness = true, zweno = true)
        end

        model = HydrostaticFreeSurfaceModel(architecture = arch,
                                                    grid = grid,
                                                tracers  = :c,
                                        tracer_advection = scheme,
                                            velocities   = PrescribedVelocityFields(u=U, v=V),
                                                coriolis = nothing,
                                                closure  = nothing,
                                                buoyancy = nothing)

        c = model.tracers.c
        set!(model, c=c₀_2D)

        simulation = Simulation(model,
                                Δt = Δt_max,
                                stop_time = end_time)

        for i = 1:end_time/Δt_max/10
            csim                               = adapt(CPU(), c)
            solution2D[(schemes[adv], Int(i))] = csim[1:N, 1:N, 1]

            for j = 1:10
                time_step!(model, Δt_max)
            end
        end
    end

    x     = adapt(CPU(), grid.xᶜᵃᵃ)[1:N]
    y     = adapt(CPU(), grid.yᵃᶜᵃ)[1:N]
    steps = ()
    anim = @animate for i ∈ 1:end_time/Δt_max/10
        plot(contourf(x, y, solution2D[(schemes[1], Int(i))], clim=(0, 1), levels = 0:0.1:1, title="Uweno"),
             contourf(x, y, solution2D[(schemes[2], Int(i))], clim=(0, 1), levels = 0:0.1:1, title="Sweno"),
             contourf(x, y, solution2D[(schemes[3], Int(i))], clim=(0, 1), levels = 0:0.1:1, title="βweno"),
             contourf(x, y, solution2D[(schemes[4], Int(i))], clim=(0, 1), levels = 0:0.1:1, title="Zweno"))
    end
    mp4(anim, "anim_2D_$gr.mp4", fps = 15)
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export id, residue_for_atom, contains, atoms

__ptr(residue::Residue) = __ptr(residue.__handle)
__const_ptr(residue::Residue) = __const_ptr(residue.__handle)

"""
Create a new residue with the given ``name`` and optional residue identifier
``id``.
"""
function Residue(name::String, id=nothing)
    if id == nothing
        ptr = @__check_ptr(lib.chfl_residue(pointer(name)))
    else
        ptr = @__check_ptr(lib.chfl_residue_with_id(pointer(name), UInt64(id)))
    end
    return Residue(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the residue at ``index`` from a ``topology``.

The residue index in the topology is not always the same as the residue
identifier.
"""
function Residue(topology::Topology, index::Integer)
    ptr = @__check_ptr(lib.chfl_residue_from_topology(__ptr(topology), UInt64(index)))
    residue = Residue(CxxPointer(ptr, is_const=true))
    copy = deepcopy(residue)
    finalize(residue)
    return copy
end

"""
Get a copy of the residue containing the atom at ``index`` in the ``topology``.

This function will return ``nothing`` if the atom is not in a residue, or if
the ``index`` is bigger than the number of atoms in the topology.
"""
function residue_for_atom(topology::Topology, index::Integer)
    ptr = lib.chfl_residue_for_atom(__ptr(topology), UInt64(index))
    if Int(ptr) == 0
        return nothing
    else
        residue = Residue(CxxPointer(ptr, is_const=true))
        copy = deepcopy(residue)
        finalize(residue)
        return copy
    end
end

"""
Get the name of a ``residue``.
"""
function name(residue::Residue)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_residue_name(
            __const_ptr(residue), buffer, size
        ))
    )
end

"""
Get the identifier of a ``residue`` in the initial topology.
"""
function id(residue::Residue)
    resid = Ref{UInt64}(0)
    __check(lib.chfl_residue_id(__const_ptr(residue), resid))
    return resid[]
end

"""
Get the number of atoms in a ``residue``.
"""
function Base.size(residue::Residue)
    count = Ref{UInt64}(0)
    __check(lib.chfl_residue_atoms_count(__const_ptr(residue), count))
    return count[]
end

"""
Add the atom at the given ``index`` in the ``residue``.
"""
function add_atom!(residue::Residue, index::Integer)
    __check(lib.chfl_residue_add_atom(__ptr(residue), UInt64(index)))
    return nothing
end

"""
Check if the atom at the given ``index`` is in the ``residue``.
"""
function contains(residue::Residue, index::Integer)
    result = Ref{UInt8}(0)
    __check(lib.chfl_residue_contains(__const_ptr(residue), UInt64(index), result))
    return convert(Bool, result[])
end

"""
Get the atoms in a ``residue``. This function returns a list of indexes.
"""
function atoms(residue::Residue)
    count = size(residue)
    result = Array{UInt64}(undef, count)
    __check(lib.chfl_residue_atoms(__const_ptr(residue), pointer(result), count))
    return result
end

"""
Set a named property for the given residue.
"""
function set_property!(residue::Residue, name::String, value)
    property = Property(value)
    __check(lib.chfl_residue_set_property(
        __ptr(residue), pointer(name), __const_ptr(property)
    ))
    return nothing
end

"""
Get a named property for the given residue.
"""
function property(residue::Residue, name::String)
    ptr = lib.chfl_residue_get_property(__const_ptr(residue), pointer(name))
    return extract(Property(CxxPointer(ptr, is_const=false)))
end

"""
Get the number of properties associated with a residue.
"""
function properties_count(residue::Residue)
    count = Ref{UInt64}(0)
    __check(lib.chfl_residue_properties_count(__const_ptr(residue), count))
    return count[]
end

"""
Get the names of all properties associated with a residue.
"""
function list_properties(residue::Residue)
    count = properties_count(residue)
    names = Array{Ptr{UInt8}}(undef, count)
    __check(lib.chfl_residue_list_properties(__const_ptr(residue), pointer(names), count))
    return map(unsafe_string, names)
end

"""
Make a deep copy of a ``residue``.
"""
function Base.deepcopy(residue::Residue)
    ptr = lib.chfl_residue_copy(__const_ptr(residue))
    return Residue(CxxPointer(ptr, is_const=false))
end


text:
using PolynomialRatios
using Polynomials
using Test

@testset "PolynomialRatios.jl" begin

    p,q = fromroots(Polynomial, [1,2,3]), fromroots(Polynomial, [2,3,4])
    r,s = ImmutablePolynomial([1,2,3], :x), ImmutablePolynomial([1,2,3], :y)
    t,u = Polynomial([1,2,pi]), ImmutablePolynomial([1.1, 2.2, 3.4], :y)

    # constructor
    @test p // q isa RationalFunction
    @test p // r isa RationalFunction
    @test_throws ArgumentError r // s

    pq = p // t # promotes to type of t
    @test pq isa RationalFunction{Float64, :x}

    # pieces
    pp,qq = p // q
    @test (pp,qq) == (p,q)

    # evaluation
    pq = p//q
    @test pq(2.5) ≈ p(2.5)/ q(2.5)
    @test pq(2) ≈ fromroots([1,3])(2) / fromroots([3,4])(2)

    # arithmetic
    rs = r // (r-1)
    x = 1.2
    @test (pq + rs)(x) ≈ (pq(x) + rs(x))
    @test (pq * rs)(x) ≈ (pq(x) * rs(x))
    @test (-pq)(x) ≈ -p(x)/q(x)

    # derivative
    pq = p // one(p)
    x = 2.3
    @test derivative(pq)(x) ≈ derivative(p)(x)
    pq = p // q
    dp,dq = derivative.((p,q))
    @test derivative(pq)(x) ≈ (dp(x)*q(x) - p(x)*dq(x))/q(x)^2

    # integral
    pq = p // (2*one(p))
    @test iszero(derivative(integrate(pq)) - pq)
    pq = p // q
    @test_throws ArgumentError integrate(pq)

    # lowest terms
    pq = p // q
    pp, qq = lowest_terms(pq)
    @test all(abs.(pp.(roots(qq))) .> 1/2)
end