text:
"""
The *LaplaceModule* contains all functionality to work with the Laplace transformation
1/x ≈ ∑ₜweight[t] * exp(-x node[t]).
The nodes and weights are contained within the *LaplacePoints* type. Hackbusch and coworkers
have optimized and pretabulated these for different x-ranges according to minmax
(Takatsuka, Ten-No, Hackbusch, JCP, 129 (2008), 044112).
Their pretabulated files can be obtained with *downloadLaplacePointsHackbusch*. Use *readLaplacePointsHackbusch* to
generate a LaplacePoints object from such a file. After downloading one can also use *findLaplacePointsHackbuschPretableLarger*/*findLaplacePointsHackbuschPretableSmaller* to obtain the best guesses for a given range. The ideal Laplace
problem considers x∈ [1,R]. For the transformation to the x∈ [A,B]-type Laplace problem utility functions are provided (*transformRangeToIdealLaplace*,*transformWeightFromIdealLaplace*,*transformNodeFromIdealLaplace*,*transformLaplacePointFromIdealLaplace*).
"""
module LaplaceModule
export LaplacePoints, downloadLaplacePointsHackbusch, findLaplacePointsHackbuschPretableLarger, findLaplacePointsHackbuschPretableSmaller, transformRangeToIdealLaplace, transformNodeFromIdealLaplace, transformWeightFromIdealLaplace, transformLaplacePointFromIdealLaplace
using ..DocumentationModule
using ..BaseModule
import Base.display
import Base.==

"""
LaplacePoints contain the collection of weights and corresponding nodes for the quadrature
according to
1/x = ∫ exp(-xt) ≈ ∑ₜweight[t] * exp(-x node[t])
"""
type LaplacePoints
  weights::Array{Float64,1}
  nodes::Array{Float64,1}
end

function display(lp::LaplacePoints)
  assert(length(lp.weights) == length(lp.nodes))
  println("  Exponential Factor         Weight")
  println("  ------------------         ------")
  for (node,weight) in zip(lp.nodes,lp.weights)
    @printf("  %.16f         %f\n",node,weight)
  end
end

"""
A dictionary mapping the number of Laplace Points to the ranges for which Hackbusch and
coworkers list pretabulated values on their webpage.
"""
const hackbuschpretables = Dict{Int64,Array{Int64,1}}(
  1=>[2:10;],
  2=>[2:10; 20:10:50],
  3=>[2:10; 20:10:100; 200],
  4=>[2:10; 20:10:100; 200:100:500],
  5=>[2:10; 20:10:100; 200:100:1000; 2000],
  6=>[2:10; 20:10:100; 200:100:1000; 2000; 3000],
  7=>[2:10; 20:10:100; 200:100:1000; 2000:1000:7000],
  8=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4; 2e4],
  9=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4],
 10=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4; 5e4; 1e5],
 11=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4; 5e4; 1e5; 2e5],
 12=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4; 5e4; 1e5:1e5:3e5],
 13=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4; 5e4; 1e5:1e5:4e5],
 14=>[      10:10:100; 200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:3e4; 5e4; 1e5:1e5:5e5; 7e5],
 15=>[      10:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6; 2e6],
 16=>[      10:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:3e6],
 17=>[      20:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:4e6],
 18=>[      40:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:5e6; 7e6],
 19=>[      50:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:5e6; 7e6; 1e7],
 20=>[      60:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:5e6; 7e6; 1e7; 2e7],
 21=>[      90:10:100; 200:100:1000; 2000:1000:1e4; 2e4:1e4:5e4; 7e4; 1e5; 2e5; 5e5; 1e6:1e6:5e6; 7e6; 1e7:1e7:3e7],
 22=>[             100:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:5e4; 7e4; 1e5:1e5:5e5; 7e5; 1e6:1e6:5e6; 7e6; 1e7:1e7:4e7],
 23=>[             200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:5e4; 7e4; 1e5:1e5:5e5; 7e5; 1e6:1e6:5e6; 7e6; 1e7:1e7:5e7; 7e7],
 24=>[             200:100:500; 700; 1000:1000:5000; 7000; 1e4:1e4:5e4; 7e4; 1e5:1e5:5e5; 7e5; 1e6:1e6:5e6; 7e6; 1e7:1e7:5e7; 7e7; 1e8],
 25=>[             200:100:500; 700; 1000:1000:10000;       2e4; 5e4;        1e5; 2e5; 5e5;    1e6; 2e6; 5e6;    1e7; 5e7;         1e8; 2e8],
 26=>[             400:100:1000;     2000:1000:10000;2e4:1e4:1e5;      2e5:1e5:5e5; 7e5;       1e6:1e6:5e6; 7e6; 1e7:1e7:5e7; 7e7; 1e8:1e8:3e8],
 27=>[                                               1e4:1e4:1e5;      2e5:1e5:1e6;            2e6:1e6:1e7; 2e7:1e7:5e7; 7e7; 1e8:1e8:4e8],
 28=>[                                                                                1e8:1e8:4e8; 7e8],
 29=>[7e8; 1e9],
 30=>[7e8; 1e9; 2e9],
 31=>[7e8; 1e9; 2e9],
 32=>[7e8; 1e9:1e9:3e9],
 33=>[7e8; 1e9:1e9:4e9],
 34=>[7e8; 1e9:1e9:5e9; 7e9],
 35=>[7e8; 1e9:1e9:3e9; 5e9; 1e10],
 36=>[1e10; 2e10],
 37=>[1e10; 2e10],
 38=>[1e10; 3e10],
 39=>[1e10; 2e10; 4e10],
 40=>[1e10; 2e10; 5e10],
 41=>[1e10; 2e10; 3e10; 7e10],
 42=>[1e10; 2e10; 5e10; 1e11],
 43=>[1e11; 2e11],
 44=>[1e11; 2e11],
 45=>[1e11; 3e11],
 46=>[1e11; 2e11; 4e11],
 47=>[1e11; 5e11],
 48=>[1e11; 2e11; 3e11; 7e11],
 49=>[1e11; 2e11; 3e11; 5e11; 1e12],
 50=>[2e8:1e8:5e8; 7e8; 1e9:1e9:7e9; 1e10:1e10:5e10; 7e10; 1e11:1e11:5e11; 7e11; 1e12; 2e12],
 51=>[1e12; 2e12],
 52=>[1e12; 3e12],
 53=>[4e11; 1e12:1e12:4e12]
)
@doc GenericCitation("http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/") hackbuschpretables

"""
Evaluates to the name of the Laplace point file which corresponds to the given number of points and limit
in Hackbusch nomenclature. See also downloadLaplacePointsHackbusch for how to obtain these files.
"""
function computeFilenameHackbusch(numberOfPoints::Integer,R::Real)
  k = numberOfPoints
  kstr = @sprintf("k%02d",k)
  Rstr = replace(@sprintf("_%.0E",R),r"\+0*(\d+)",s"\1")
  filename = "1_x$kstr$Rstr"
  return filename
end

"""
Download all files with pretabulated Laplace points to the given folder (default:"hackbusch")
"""
function downloadLaplacePointsHackbusch(targetdir::AbstractString="hackbusch")
  mkdir(targetdir)
  for (k,Rlist) in hackbuschpretables
    for R in Rlist
      filename = computeFilenameHackbusch(k,R)
      println(filename)
      download("http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/$filename","$targetdir/$filename")
    end
  end
end
@doc GenericCitation("http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/") downloadLaplacePointsHackbusch

"""
Generate a LaplacePoints object by reading in the specified file.
"""
function readLaplacePointsHackbusch(filename::AbstractString)
  k = countlines(filename)/2
  fd = open(filename)
  result = LaplacePoints(Array{Float64,1}(),Array{Float64,1}())
  for line in 1:k
    regmatch = match(Regex(" *(?P<weight>$floatregex) .*omega.*"),readline(fd))
    push!(result.weights,float(regmatch[:weight]))
  end
  for line in k+1:2k
    regmatch = match(Regex(" *(?P<node>$floatregex) .*alpha.*"),readline(fd))
    push!(result.nodes,float(regmatch[:node]))
  end
  close(fd)
  return result
end

function readLaplacePointsHackbusch(numberOfPoints::Integer,R::Real,dir::AbstractString="hackbusch")
  filename = computeFilenameHackbusch(numberOfPoints,R)
  readLaplacePointsHackbusch("$dir/$filename")
end
@doc """
Alternative to specifying the file directly one can also specify the numberOfPoints, upper limit,
and directory name (default:"hackbusch") to load the corresponding pretable file..
""" readLaplacePointsHackbusch

"""
Returns the LaplacePoints for the smallest range larger than the requested one for which
pretabulated Laplace points can be found. If the requested range is larger than any pretabulated one,
we expect the largest pretabulated one to be Rₖ, so the function issues a warning and
returns the LaplacePoints for that range.
"""
function findLaplacePointsHackbuschPretableLarger(numberOfPoints::Integer,R::Float64,dir::AbstractString="hackbusch")
  for r in hackbuschpretables[numberOfPoints]
    if R <= r
      return (readLaplacePointsHackbusch(numberOfPoints,r,dir),r)
    end
  end
  # Expect R > Rₖ, so just return largest available table for the given number of points
  warn("numberOfPoints not sufficiently large ?!")
  r = hackbuschpretables[numberOfPoints][end]
  return (readLaplacePointsHackbusch(numberOfPoints,r,dir),r)
end

"""
Returns the LaplacePoints for the largest range smaller than the requested one for which
pretabulated Laplace points can be found.
"""
function findLaplacePointsHackbuschPretableSmaller(numberOfPoints::Integer,R::Float64,dir::AbstractString="hackbusch")
  for r in reverse(hackbuschpretables[numberOfPoints])
    if R >= r
      return (readLaplacePointsHackbusch(numberOfPoints,r,dir),r)
    end
  end
  error("No pretabulated value for $numberOfPoints points and a Range of [1., $R] or smaller.")
end

"""
For the ideal Laplace problem the values x must be known to be in the range [1,R].
If 1/x can in general be within [A,B], then all nodes and weights must be scaled with 1/A,
so that x' := x/A lies within [1,B/A (=:R)].

Consider: 1/x = ∑ w exp(-α x)  =>  A/x = ∑ A w exp(-α x)
                               =>  A/x = ∑ (A w) exp(-(A α) x/A)
			       =>  1/x' = ∑ w' exp(-α' x')
"""
function transformRangeToIdealLaplace(A::Float64,B::Float64)
  return (1,B/A)
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformNodeFromIdealLaplace(α::Float64,A::Float64)
  return α/A
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformWeightFromIdealLaplace(w::Float64,A::Float64)
  return w/A
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformLaplacePointFromIdealLaplace(lp::LaplacePoints, A)
  LaplacePoints(map(x->transformWeightFromIdealLaplace(x,A),lp.weights), map(x->transformNodeFromIdealLaplace(x,A),lp.nodes))
end

"""
1/x ≈ ∑ᵢ wᵢ exp(-αᵢ x)
"""
function computeInverseByLaplaceApproximation(x,lp::LaplacePoints)
  sum=0.
  for (node,weight) in zip(lp.nodes,lp.weights)
    sum += weight*exp(-node*x)
  end
  return sum
end
@doc JournalCitation(["A. Takatsuka", "S. Ten-No", "W. Hackbusch"],"J. Chem. Phys.",129,044112,2008) computeInverseByLaplaceApproximation

"""
1/(x² + ω²) ≈ ∑ᵢ wᵢ sin(ω αᵢ)/ω exp(-αᵢ x)
"""
function computeRPADenominatorByDoubleLaplace(x::Float64,ω::Float64,lp::LaplacePoints)
  sum = 0.
  for (node,weight) in zip(lp.nodes,lp.weights)
    sum += weight * sin(ω*node)/ω * exp(-node*x)
  end
  return sum
end
@doc JournalCitation(["H. F. Schurkus","C. Ochsenfeld"],"J. Chem. Phys.",144,031101,2016) computeRPADenominatorByDoubleLaplace

==(x::LaplacePoints, y::LaplacePoints) = x.weights == y.weights && x.nodes == y.nodes

end # LaplaceModule


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

name = "p3a"
version = v"1.0.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sandialabs/p3a.git", "271cd1d77fbc850aa4ddfb015b59176926670b3b")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/p3a/
mkdir build && cd build

CMAKE_FLAGS=()

#help find mpi on mingw, subset of https://github.com/JuliaPackaging/Yggdrasil/blob/b4fdb545c3954cff218051d7520c7418991d3416/T/TauDEM/build_tarballs.jl#L28-L53
if [[ "$target" == x86_64-w64-mingw32 ]]; then
    CMAKE_FLAGS+=(
        -DMPI_HOME=${prefix}
        -DMPI_GUESS_LIBRARY_NAME=MSMPI
        -DMPI_C_LIBRARIES=msmpi64
        -DMPI_CXX_LIBRARIES=msmpi64
    )
elif [[ ${target} == *-apple-* ]]; then
    CMAKE_FLAGS+=(
        -DMPI_C_ADDITIONAL_INCLUDE_DIRS=''
        -DMPI_C_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi'
        -DMPI_CXX_ADDITIONAL_INCLUDE_DIRS=''
        -DMPI_CXX_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi'
    )
fi

cmake .. \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DBUILD_TESTING=OFF \
-DBUILD_SHARED_LIBS=ON \
-DKokkos_COMPILE_LANGUAGE=CXX \
"${CMAKE_FLAGS[@]}"

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))
#Kokkos is only available on 64 bit
filter!(p -> nbits(p) != 32, platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Kokkos_jll", uuid="c1216c3d-6bb3-5a2b-bbbf-529b35eba709"))
    Dependency(PackageSpec(name="MPICH_jll"); platforms=filter(!Sys.iswindows, platforms))
    Dependency(PackageSpec(name="MicrosoftMPI_jll"); platforms=filter(Sys.iswindows, platforms))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isfreebsd, platforms))
]

# Build the tarballs, and possibly a `build.jl` as well.
#needed a c++17 compiler, 7 and 8 kept refusing to acknowledge std::byte and a few other template things for some reason?!
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"9.1.0")

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("../server/simulation/simulation.jl")

using .Simulation

function test_axon_basic()
    println("Testing Basic Axon")

    axon = Simulation.DirectAxon()

    println("Passed")
end

function test_axon_single_stream()
    println("Testing Axon single stream")
    stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    axon = Simulation.DirectAxon()
    Simulation.add_stream!(axon, stream)

    @assert Simulation.output(axon) == 0 "Axon output not zero at start"

    expected_output = [0,0,0,0,0,1,0,0,1,0]
    for i in 1:10
        Simulation.pre!(axon)
        @assert Simulation.output(axon) == expected_output[i] "Expected output did not match: $i"
    end

    println("Passed")
end

function test_axon_two_streams()
    println("Testing Axon two streams")
    stream = Simulation.PoissonStream(UInt64(13163), 0.2)
    stream2 = Simulation.PoissonStream(UInt64(11862), 0.2)

    axon = Simulation.DirectAxon()
    Simulation.add_stream!(axon, stream)
    Simulation.add_stream!(axon, stream2)
    # for i in 1:10
    #     Simulation.step!(stream)
    #     println(Simulation.output(stream))
    # end
    # println("--------------")
    # for i in 1:10
    #     Simulation.step!(stream2)
    #     println(Simulation.output(stream2))
    # end
    # 0000010010  = stream
    # 0001001000  = stream2
    # 0001011010  = output

    expected_output = [0,0,0,1,0,1,1,0,1,0]
    for i in 1:10
        # println("i: ", i)
        Simulation.pre!(axon)
        # println(Simulation.output(axon))
        @assert Simulation.output(axon) == expected_output[i] "Expected output did not match: $i"
    end

    println("Passed")
end

function test_axon_delay_tiny()
    println("Testing Axon tiny delay")
    stream = Simulation.PoissonStream(UInt64(11862), 0.2)

    axon = Simulation.TinyDelayAxon(7)
    Simulation.add_stream!(axon, stream)

    # 0001001000  = stream

    # expected_output = [0,0,0,0,0,1,0,0,1,0]
    for i in 1:80
        Simulation.pre!(axon)
        print(bitstring(axon.shift_reg), " : $i pre -- ")
        Simulation.post!(axon)
        println(bitstring(axon.shift_reg), " : $i post")
        # @assert Simulation.output(axon) == expected_output[i] "Expected output did not match: $i"
    end

    println("Passed")
end

function test_axon_delay_small()
    println("Testing Axon small delay")
    stream = Simulation.PoissonStream(UInt64(11862), 0.2)

    axon = Simulation.SmallDelayAxon(63)
    Simulation.add_stream!(axon, stream)

    # 0001001000  = stream

    # expected_output = [0,0,0,0,0,1,0,0,1,0]
    for i in 1:80
        Simulation.pre!(axon)
        print(bitstring(axon.shift_reg), " : $i pre -- ")
        Simulation.post!(axon)
        println(bitstring(axon.shift_reg), " : $i post")
        # @assert Simulation.output(axon) == expected_output[i] "Expected output did not match: $i"
    end

    println("Passed")
end

function test_axon_delay_small2()
    println("Testing Axon small delay")
    stream = Simulation.PoissonStream(UInt64(13163), 0.2)
    stream2 = Simulation.PoissonStream(UInt64(11862), 0.2)

    axon = Simulation.MediumDelayAxon(32)
    Simulation.add_stream!(axon, stream)
    Simulation.add_stream!(axon, stream2)

    # 0001001000  = stream

    # expected_output = [0,0,0,0,0,1,0,0,1,0]
    for i in 1:80
        Simulation.pre!(axon)
        print(bitstring(axon.shift_reg), " : $i pre -- ")
        Simulation.post!(axon)
        println(bitstring(axon.shift_reg), " : $i post")
        # @assert Simulation.output(axon) == expected_output[i] "Expected output did not match: $i"
    end

    println("Passed")
end

test_axon_delay_small()


text:
using DNSS, Test

libraries = ["MinkowskiSphericalSymmetry"]

for file in libraries
    @info "Testing $file"
    include("$(file)Test.jl")
end


text:
module Highway

using JLD

import Base:
        show

export
        Vehicle, Highway1D, Highway2D,
        Empty_Cell!, Obstacle!, Change_Vehicle!, Insert_Vehicle!,  Vehicle_Behind, Vehicle_Ahead,
        Pos_next, Pos_right, Pos_left, Pos_left9,
        Merge!, Merge_Left_Right!, Merge_Right_Left!, safe_distance, Ramp!,
        AccelerationNoise, DecelerationMove, AccelerationNoiseS1, AccelerationNoiseS2,
        Measure_Fluxes!, Measure_Speeds!, add_Fluxes!, add_Speeds!, Measure_Densities!, Measure!, Measure_s!,
        Sense1, Sense2, Sections1_CR, Sections2_CR,
        S1_ramp, S1_ramp_ampli, S1_ramp_local, P_ramp_S1, P_ampli_S1, P_local_S1,
        S2_ramp, S2_ramp_ampli, S2_ramp_local, P_ramp_S2, P_ampli_S2, P_local_S2

include("types.jl")
include("aux-functions.jl")
include("AC-functions.jl")
include("merging-ramps.jl")
include("measurements.jl")
include("senses.jl")

end # module Highway


text:
@testset "dict to dataset roundtrip" begin
    rng = MersenneTwister(42)
    J = 8
    K = 6
    nchains = 4
    ndraws = 500
    vars = Dict(
        "a" => randn(rng, nchains, ndraws),
        "b" => randn(rng, nchains, ndraws, J, K)
    )
    coords = Dict("bi" => 1:J, "bj" => 1:K)
    dims = Dict("b" => ["bi", "bj"])
    attrs = Dict("mykey" => 5)

    ds = ArviZ.dict_to_dataset(vars; coords = coords, dims = dims, attrs = attrs)
    vars2, kwargs = ArviZ.dataset_to_dict(ds)
    for (k, v) in vars
        @test k ∈ keys(vars2)
        @test vars2[k] ≈ v
    end
    @test kwargs.coords == coords
    @test kwargs.dims == dims
    for (k, v) in attrs
        @test k ∈ keys(kwargs.attrs)
        @test kwargs.attrs[k] == v
    end
end


text:
function init_VariableStats()
  VarStats = VariableStats(
  zeros(Float64, 1,1), # Histogram
  [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],
  Array{Dict{Tuple{Int,Int,Int}, Float64},1}()
  )
  return VarStats
end

export init_VariableStats

function getPosNegCentroids(CCobj::CCstats, VarStats::VariableStats, QuantileScores::Array{Float32, 3}; threshigh::Float64 = 0.0, threslow::Float64 = 0.0, min_nobs_within_thres::Int = 100)
  ValueWeighting = true
  if size(QuantileScores) != size(CCobj.LabelCube) error("size  of CCobj.LabelCube and QuantileScores not matching") end
  if length(CCobj.affectedLonLatTimeScoreValue) == 0 error("CCobj.affectedLonLatTimeScoreValue is empty") end
  VarStats.High_LONcenters = zeros(Float64, CCobj.numExtremes)
  VarStats.High_LATcenters = zeros(Float64, CCobj.numExtremes)
  VarStats.High_TIMEcenters = zeros(Float64, CCobj.numExtremes)
  sweightsHigh = zeros(Float64, CCobj.numExtremes)
  valweightsHigh = zeros(Float64, CCobj.numExtremes)
  Div = zeros(Int, CCobj.numExtremes)
  VarStats.Low_LONcenters = zeros(Float64, CCobj.numExtremes)
  VarStats.Low_LATcenters = zeros(Float64, CCobj.numExtremes)
  VarStats.Low_TIMEcenters = zeros(Float64, CCobj.numExtremes)
  sweightsLow = zeros(Float64, CCobj.numExtremes)
  valweightsLow = zeros(Float64, CCobj.numExtremes)
  Div2 = zeros(Int, CCobj.numExtremes)
  for i = 1:length(CCobj.affectedLonLatTimeScoreValue)
    for lonlattimeidx = keys(CCobj.affectedLonLatTimeScoreValue[i])
      (lonidx, latidx, timeidx) = lonlattimeidx
      if ValueWeighting
        val = QuantileScores[lonidx, latidx, timeidx]
      else
        val = 1
      end
      if val > threshigh
        innerstCentroids!(sweightsHigh, valweightsHigh, VarStats.High_LATcenters, VarStats.High_LONcenters, VarStats.High_TIMEcenters, CCobj, lonlattimeidx, val, i)
        Div[i] += 1
      elseif val < threslow
        innerstCentroids!(sweightsLow, valweightsLow, VarStats.Low_LATcenters, VarStats.Low_LONcenters, VarStats.Low_TIMEcenters, CCobj, lonlattimeidx, val, i)
        Div2[i] += 1
      end
    end
  end
  broadcast!(/ ,VarStats.High_LATcenters, VarStats.High_LATcenters, sweightsHigh)
  #broadcast!(acosd ,VarStats.High_LATcenters,VarStats.High_LATcenters)
  broadcast!(/, VarStats.High_LONcenters, VarStats.High_LONcenters, valweightsHigh)
  broadcast!(/ ,VarStats.High_TIMEcenters, VarStats.High_TIMEcenters, valweightsHigh)
  broadcast!(/ ,VarStats.Low_LATcenters, VarStats.Low_LATcenters, sweightsLow)
  #broadcast!(acosd ,VarStats.Low_LATcenters,VarStats.Low_LATcenters)
  broadcast!(/, VarStats.Low_LONcenters, VarStats.Low_LONcenters, valweightsLow)
  broadcast!(/ ,VarStats.Low_TIMEcenters, VarStats.Low_TIMEcenters, valweightsLow)
  for j = 1:length(Div)
    if Div[j] < min_nobs_within_thres
      VarStats.High_LATcenters[j] = NaN; VarStats.High_LONcenters[j] = NaN; VarStats.High_TIMEcenters[j] = NaN
    end
    if Div2[j] < min_nobs_within_thres
      VarStats.Low_LATcenters[j] = NaN; VarStats.Low_LONcenters[j] = NaN; VarStats.Low_TIMEcenters[j] = NaN
    end
  end
  return Centroids = (VarStats.High_LONcenters, VarStats.High_LATcenters, VarStats.High_TIMEcenters, VarStats.Low_LONcenters, VarStats.Low_LATcenters, VarStats.Low_TIMEcenters)
end

export getPosNegCentroids

function spatDist(lon1, lon2, lat1, lat2)
    meanlat = (lat1 .+ lat2) / 2
    spatD = sqrt.((lon1 .- lon2) .* (lon1 .- lon2) + (cosd.(meanlat) .* (lat1 .- lat2)) .* (cosd.(meanlat) .* (lat1 .- lat2))) .* 111.324
    return spatD
end

export spatDist

function getDistances(VarStats::VariableStats)
  VarStats.SpatialDistance = spatDist(VarStats.Low_LONcenters, VarStats.High_LONcenters, VarStats.Low_LATcenters, VarStats.High_LATcenters)
  VarStats.TemporalDistance = (VarStats.High_TIMEcenters-VarStats.Low_TIMEcenters)
  return VarStats.SpatialDistance, VarStats.TemporalDistance # D_day negative: time of negative response > time of positive response, i.e. negative response after positive
end

export getDistances

function getPosNegSums(SumScoresDict::Array{Dict{Tuple{Int64,Int64},Float32},1})
    posres = zeros(eltype(values(SumScoresDict[1])), length(SumScoresDict))
    negres = zeros(eltype(values(SumScoresDict[1])), length(SumScoresDict))
    for lab = 1:length(SumScoresDict)
        if isassigned(SumScoresDict, lab)
            for v = values(SumScoresDict[lab])
              if !isnan(v)
                if v < 0
                    negres[lab] += v
                elseif v < 10^38
                    posres[lab] += v
                end
              end
            end
        end
    end
    bil = posres + negres
    return bil, posres, negres, posres ./ negres
end

function getBalance(VarStats::VariableStats, CCobj::CCstats, inMemVarAnomalies; multiplyByConst::Bool = true)
    SumScoresDict = get_ScoresSumAreaWeighted(CCobj, inMemVarAnomalies);
    VarStats.Balance, VarStats.PosSum, VarStats.NegSum, VarStats.Comp = getPosNegSums(SumScoresDict)
    if multiplyByConst
        kmconst = (111.324 * 1000 * abs(CCobj.Latitudes[1]-CCobj.Latitudes[2])) * (111.324 * 1000 * abs(CCobj.Longitudes[1]-CCobj.Longitudes[2])) * abs(CCobj.Times[2] - CCobj.Times[1])
        VarStats.Balance = VarStats.Balance .* kmconst
        VarStats.PosSum =  VarStats.PosSum .* kmconst
        VarStats.NegSum =  VarStats.NegSum .* kmconst
    end
    return VarStats.Balance, VarStats.NegSum, VarStats.PosSum, VarStats.Comp
end

export getBalance

function getWeightedSumOfValuesPerMask(CCobj::CCstats, landcover, xin)
    if size(landcover) != size(xin)[1:2] || size(CCobj.LabelCube) != size(xin) error("size not matching") end
    xout = Array{Dict{Int32, Float32}, 1}(undef, CCobj.numExtremes)
    concomp = CCobj.LabelCube
    if size(concomp) != size(xin) error("size of input arrays does not match") end
    if size(concomp)[1:2] != size(landcover) error("size of input arrays does not match") end
    for k = 1:size(concomp, 3)
        for j = 1:size(concomp, 2)
            for i = 1:size(concomp, 1)
                if(concomp[i, j, k] > 0)
                    label = concomp[i, j, k]
                    if !isnan(xin[i, j, k])
                        if isassigned(xout, label)
                            if(haskey(xout[label], landcover[i,j]))
                                xout[label][landcover[i,j]] = xout[label][landcover[i,j]] + xin[i, j, k] * CCobj.CosLatitudes[j]
                            else
                                xout[label][landcover[i,j]] = xin[i, j, k] * CCobj.CosLatitudes[j]
                            end
                        else
                            xout[label] = Dict(landcover[i,j] => xin[i, j, k] * CCobj.CosLatitudes[j])
                        end
                    end
                end
            end
        end
    end
    return(xout)
end

export getWeightedSumOfValuesPerMask

function getWeightedAverageOfValuesPerMask(CCobj::CCstats, landcover, xin)
    if size(landcover) != size(xin)[1:2] || size(CCobj.LabelCube) != size(xin) error("size not matching") end
    xout = Array{Dict{Int32, Float32}, 1}(undef, CCobj.numExtremes)
    xdiv = Array{Dict{Int32, Float32}, 1}(undef, CCobj.numExtremes)
    concomp = CCobj.LabelCube
    if size(concomp) != size(xin) error("size of input arrays does not match") end
    if size(concomp)[1:2] != size(landcover) error("size of input arrays does not match") end
    for k = 1:size(concomp, 3)
        for j = 1:size(concomp, 2)
            for i = 1:size(concomp, 1)
                if(concomp[i, j, k] > 0)
                    label = concomp[i, j, k]
                    if !isnan(xin[i, j, k])
                        if isassigned(xout, label)
                            if(haskey(xout[label], landcover[i,j]))
                                xout[label][landcover[i,j]] = xout[label][landcover[i,j]] + xin[i, j, k] * CCobj.CosLatitudes[j]
                                xdiv[label][landcover[i,j]] = xdiv[label][landcover[i,j]] + CCobj.CosLatitudes[j]
                            else
                                xout[label][landcover[i,j]] = xin[i, j, k] * CCobj.CosLatitudes[j]
                                xdiv[label][landcover[i,j]] = CCobj.CosLatitudes[j]
                            end
                        else
                            xout[label] = Dict(landcover[i,j] => xin[i, j, k] * CCobj.CosLatitudes[j])
                            xdiv[label] = Dict(landcover[i,j] => CCobj.CosLatitudes[j])
                        end
                    end
                end
            end
        end
    end
    for l = 1:length(xout)
        if isassigned(xout,l )
            for m = keys(xout[l])
                xout[l][m] = xout[l][m] / xdiv[l][m]
            end
        end
    end
    return(xout)
end

export getWeightedAverageOfValuesPerMask

function get_ScoresSumAreaWeighted(CCobj::CCstats, Scores)
    ScoresSum = Array{Dict{Tuple{Int,Int}, eltype(Scores)}, 1}(undef, CCobj.numExtremes)
    get_ScoresSumAreaWeighted!(ScoresSum, CCobj, Scores)
  return(ScoresSum)
end

get_ScoresSumAreaWeighted

function get_ScoresSumAreaWeighted!(xout::Array{Dict{Tuple{Int64,Int64},tp},1}, CCobj::CCstats, Scores::AbstractArray{tp, 3}) where {tp}
  if CCobj.LabelCube == zeros(Int, 1, 1, 1) error("empty CCobj.LabelCube") end
  #CCobj.affectedLonLatIdxCounts = Array(Dict{Tuple{Int,Int}, Int}, CCobj.numExtremes)
  #xout = CCobj.affectedLonLatIdxCounts
  for k = 1:size(CCobj.LabelCube, 3)
    for j = 1:size(CCobj.LabelCube, 2)
      for i = 1:size(CCobj.LabelCube, 1)
        l = CCobj.LabelCube[i,j,k]
        if l > 0
          if isassigned(xout, l)
            if haskey(xout[l], (i, j))
              xout[l][(i,j)] = xout[l][(i,j)] + Scores[i,j,k] * CCobj.CosLatitudes[j]
            else
              xout[l][(i,j)] = Scores[i,j,k] * CCobj.CosLatitudes[j]
            end
          else
            xout[l] = Dict((i, j) => Scores[i,j,k] * CCobj.CosLatitudes[j]) # lon -> i, lat -> j
          end
        end
      end
    end
  end
  return xout
end

 export get_ScoresSumAreaWeighted!


text:
# Nucleic Acid
# ============
#
# DNA and RNA types.
#
# This file is a part of the BioJulia ecosystem.
# License is MIT: https://github.com/BioJulia/NucleicAcids.jl/blob/master/LICENSE.md

# NucleicAcid Encoding
# -------------------
#
# Unambiguous nucleotides are represented in one-hot encoding as follows:
#
#   | NucleicAcid | Bits |
#   | ----------- | ---- |
#   |     A       | 0001 |
#   |     C       | 0010 |
#   |     G       | 0100 |
#   |    T/U      | 1000 |
#
# Ambiguous nucleotides are bitwise OR of these four nucleotides. For example, R
# , A or G, is represented as 0101 (= A: 0001 | G: 0100). The gap symbol is
# always 0000.  The meaningful four bits are stored in the least significant
# bits of a byte.

"""
An abstract nucleic acid type.
"""
abstract type NucleicAcid <: BioSymbol end

"""
A deoxyribonucleic acid type.
"""
primitive type DNA <: NucleicAcid 8 end

"""
A ribonucleic acid type.
"""
primitive type RNA <: NucleicAcid 8 end

prefix(::DNA) = "DNA"
prefix(::RNA) = "RNA"
type_text(x::NucleicAcid) = prefix(x)
isterm(symbol::NucleicAcid) = false

# Conversion from/to integers
# ---------------------------

Base.convert(::Type{T}, nt::UInt8) where T <: NucleicAcid = reinterpret(T, nt)
Base.convert(::Type{UInt8}, nt::T) where T <: NucleicAcid = reinterpret(UInt8, nt)
Base.convert(::Type{T}, nt::S) where {T <: Number, S <: NucleicAcid} = convert(T, convert(UInt8, nt))
Base.convert(::Type{S}, nt::T) where {T <: Number, S <: NucleicAcid} = convert(S, convert(UInt8, nt))
DNA(nt::Integer) = convert(DNA, nt)
RNA(nt::Integer) = convert(RNA, nt)

Base.convert(::Type{DNA}, nt::RNA) = DNA(convert(UInt8, nt))
DNA(nt::RNA) = convert(DNA, nt)
Base.convert(::Type{RNA}, nt::DNA) = RNA(convert(UInt8, nt))
RNA(nt::DNA) = convert(RNA, nt)

# Conversion from/to characters
# -----------------------------

function Base.convert(::Type{DNA}, c::Char)
    if c > '\uff'
        throw(InexactError(:convert, DNA, c))
    end
    @inbounds dna = char_to_dna[convert(Int, c) + 1]
    if !isvalid(DNA, dna)
        throw(InexactError(:convert, DNA, c))
    end
    return reinterpret(DNA, dna)
end
DNA(c::Char) = convert(DNA, c)

function Base.convert(::Type{RNA}, c::Char)
    if c > '\uff'
        throw(InexactError(:convert, RNA, c))
    end
    @inbounds rna = char_to_rna[convert(Int, c) + 1]
    if !isvalid(RNA, rna)
        throw(InexactError(:convert, RNA, c))
    end
    return reinterpret(RNA, rna)
end
RNA(c::Char) = convert(RNA, c)

function Base.convert(::Type{Char}, nt::DNA)
    return dna_to_char[convert(UInt8, nt) + 1]
end
Char(nt::DNA) = convert(Char, nt)

function Base.convert(::Type{Char}, nt::RNA)
    return rna_to_char[convert(UInt8, nt) + 1]
end
Char(nt::RNA) = convert(Char, nt)

# Encoding of DNA and RNA NucleicAcids
# ------------------------------------

"""
Lookup table used for converting characters to DNA symbol values

The provided `convert` method should be used rather than this table, but you can
use it if you insist and know what your are doing.

!!! note
    The array is indexed by converting a character to an integer. When indexed, it
    returns a UInt8 with the bit pattern on the corresponding nucleic acid.
    The `convert(DNA, x)` method does this for you.

!!! warning
    If you index this array with a character that is greater than '\uff', then
    you will get a bounds error. The `convert(DNA, x)` method checks such things
    to avoid this for you.

!!! warning
    If you index this array with a character that does not have a corresonding
    DNA symbol, then you get a byte with the bit pattern `0x80`, which is an
    invalid DNA symbol and will be of no use to you. The `convert(DNA, x)`
    checks such things for you and throws an exception gracefully if such a
    situation arises.
"""
const char_to_dna = [0x80 for _ in 0x00:0xff]

"""
Lookup table for converting DNA symbol values to characters

The provided `convert` method should be used rather than this table, but you can
use it if you insist and know what your are doing.

!!! note
    The array is indexed by reinterpreting a DNA symbol value as an UInt8.
    When indexed, it returns the character corresponding to the symbol.
    The `convert(Char, x::DNA)` method does this for you.

!!! warning
    If you index this array with an invalid DNA symbol, then you will hit a
    bounds error. If you construct DNA symbols properly, then this scenario
    should never occur.
"""
dna_to_char

# Derived from "The DDBJ/ENA/GenBank Feature Table Definition"
# §7.4.1 Nucleotide base code (IUPAC)
# http://www.insdc.org/documents/feature_table.html#7.4.1
const dna_to_char = let
    chararray = Vector{Char}(undef, 16)
    for (char, doc, bits) in [
        ('-', "DNA Gap",                                   0b0000),
        ('A', "DNA Adenine",                               0b0001),
        ('C', "DNA Cytosine",                              0b0010),
        ('G', "DNA Guanine",                               0b0100),
        ('T', "DNA Thymine",                               0b1000),
        ('M', "DNA Adenine or Cytosine",                   0b0011),
        ('R', "DNA Adenine or Guanine",                    0b0101),
        ('W', "DNA Adenine or Thymine",                    0b1001),
        ('S', "DNA Cytosine or Guanine",                   0b0110),
        ('Y', "DNA Cytosine or Thymine",                   0b1010),
        ('K', "DNA Guanine or Thymine",                    0b1100),
        ('V', "DNA Adenine, Cytosine or Guanine",          0b0111),
        ('H', "DNA Adenine, Cytosine or Thymine",          0b1011),
        ('D', "DNA Adenine, Guanine or Thymine",           0b1101),
        ('B', "DNA Cytosine, Guanine or Thymine",          0b1110),
        ('N', "DNA Adenine, Cytosine, Guanine or Thymine", 0b1111)]
        var = Symbol("DNA_", char != '-' ? char : "Gap")
        @eval begin
            @doc $(doc) const $(var) = reinterpret(DNA, $(bits))
            char_to_dna[$(convert(Int, char) + 1)] = char_to_dna[$(convert(Int, lowercase(char)) + 1)] = $(bits)
            $(chararray)[$(convert(Int, bits) + 1)] = $(char)
        end
    end
    Tuple(chararray)
end

@eval function alphabet(::Type{DNA})
    return $(tuple([reinterpret(DNA, x) for x in 0b0000:0b1111]...))
end

"""
    alphabet(DNA)

Get all symbols of `DNA` in sorted order.

Examples
--------

```jldoctest
julia> alphabet(DNA)
(DNA_Gap, DNA_A, DNA_C, DNA_M, DNA_G, DNA_R, DNA_S, DNA_V, DNA_T, DNA_W, DNA_Y, DNA_H, DNA_K, DNA_D, DNA_B, DNA_N)

julia> issorted(alphabet(DNA))
true

```
"""
alphabet(::Type{DNA})

"""
    ACGT

Unambiguous DNA.

Examples
--------

```jldoctest
julia> ACGT
(DNA_A, DNA_C, DNA_G, DNA_T)

```
"""
const ACGT = (DNA_A, DNA_C, DNA_G, DNA_T)

"""
    ACGTN

Unambiguous DNA and `DNA_N`.

Examples
--------

```jldoctest
julia> ACGTN
(DNA_A, DNA_C, DNA_G, DNA_T, DNA_N)

```
"""
const ACGTN = (DNA_A, DNA_C, DNA_G, DNA_T, DNA_N)

"""
Lookup table used for converting characters to RNA symbol values

The provided `convert` method should be used rather than this table, but you can
use it if you insist and know what your are doing.

!!! note
    The array is indexed by converting a character to an integer. When indexed, it
    returns a UInt8 with the bit pattern on the corresponding nucleic acid.
    The `convert(RNA, x)` method does this for you.

!!! warning
    If you index this array with a character that is greater than '\uff', then
    you will get a bounds error. The `convert(RNA, x)` method checks such things
    to avoid this for you.

!!! warning
    If you index this array with a character that does not have a corresonding
    RNA symbol, then you get a byte with the bit pattern `0x80`, which is an
    invalid RNA symbol and will be of no use to you. The `convert(RNA, x)`
    checks such things for you and throws an exception gracefully if such a
    situation arises.
"""
const char_to_rna = [0x80 for _ in 0x00:0xff]

"""
Lookup table for converting RNA symbol values to characters

The provided `convert` method should be used rather than this table, but you can
use it if you insist and know what your are doing.

!!! note
    The array is indexed by reinterpreting a RNA symbol value as an UInt8.
    When indexed, it returns the character corresponding to the symbol.
    The `convert(Char, x::RNA)` method does this for you.

!!! warning
    If you index this array with an invalid RNA symbol, then you will hit a
    bounds error. If you construct RNA symbols properly, then this scenario
    should never occur.
"""
rna_to_char

const rna_to_char = let
    chararray = Vector{Char}(undef, 16)
    for (char, doc, dna) in [
        ('-', "RNA Gap",                                  DNA_Gap),
        ('A', "RNA Adenine",                              DNA_A  ),
        ('C', "RNA Cytosine",                             DNA_C  ),
        ('G', "RNA Guanine",                              DNA_G  ),
        ('U', "RNA Uracil",                               DNA_T  ),
        ('M', "RNA Adenine or Cytosine",                  DNA_M  ),
        ('R', "RNA Adenine or Guanine",                   DNA_R  ),
        ('W', "RNA Adenine or Uracil",                    DNA_W  ),
        ('S', "RNA Cytosine or Guanine",                  DNA_S  ),
        ('Y', "RNA Cytosine or Uracil",                   DNA_Y  ),
        ('K', "RNA Guanine or Uracil",                    DNA_K  ),
        ('V', "RNA Adenine, Cytosine or Guanine",         DNA_V  ),
        ('H', "RNA Adenine, Cytosine or Uracil",          DNA_H  ),
        ('D', "RNA Adenine, Guanine or Uracil",           DNA_D  ),
        ('B', "RNA Cytosine, Guanine or Uracil",          DNA_B  ),
        ('N', "RNA Adenine, Cytosine, Guanine or Uracil", DNA_N  )]
        var = Symbol("RNA_", char != '-' ? char : "Gap")
        @eval begin
            @doc $(doc) const $(var) = reinterpret(RNA, $(dna))
            char_to_rna[$(convert(Int, char) + 1)] = char_to_rna[$(convert(Int, lowercase(char) + 1))] = reinterpret(UInt8, $(dna))
            $(chararray)[$(convert(Int, dna) + 1)] = $(char)
        end
    end
    Tuple(chararray)
end

@eval function alphabet(::Type{RNA})
    return $(tuple([reinterpret(RNA, x) for x in 0b0000:0b1111]...))
end

"""
    alphabet(RNA)

Get all symbols of `RNA` in sorted order.

Examples
--------

```jldoctest
julia> alphabet(RNA)
(RNA_Gap, RNA_A, RNA_C, RNA_M, RNA_G, RNA_R, RNA_S, RNA_V, RNA_U, RNA_W, RNA_Y, RNA_H, RNA_K, RNA_D, RNA_B, RNA_N)

julia> issorted(alphabet(RNA))
true

```
"""
alphabet(::Type{RNA})

"""
    ACGU

Unambiguous RNA.

Examples
--------

```jldoctest
julia> ACGU
(RNA_A, RNA_C, RNA_G, RNA_U)

```
"""
const ACGU = (RNA_A, RNA_C, RNA_G, RNA_U)

"""
    ACGUN

Unambiguous RNA and `RNA_N`.

Examples
--------

```jldoctest
julia> ACGUN
(RNA_A, RNA_C, RNA_G, RNA_U, RNA_N)

```
"""
const ACGUN = (RNA_A, RNA_C, RNA_G, RNA_U, RNA_N)

function Base.:~(x::N) where N <: NucleicAcid
    return reinterpret(N, ~reinterpret(UInt8, x) & 0b1111)
end

function Base.:|(x::N, y::N) where N <: NucleicAcid
    return reinterpret(N, reinterpret(UInt8, x) | reinterpret(UInt8, y))
end

function Base.:-(x::N, y::Integer) where N <: NucleicAcid
    return x + (-y)
end

function Base.:+(x::N, y::Integer) where N <: NucleicAcid
    return reinterpret(N, (convert(UInt8, x) + y % UInt8) & 0b1111)
end

"""
    gap(DNA)

Return `DNA_Gap`.
"""
gap(::Type{DNA}) = DNA_Gap

"""
    gap(RNA)

Return `RNA_Gap`.
"""
gap(::Type{RNA}) = RNA_Gap

"""
    isGC(nt::NucleicAcid)

Test if `nt` is surely either guanine or cytosine.
"""
function isGC(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b1001) == 0
end

"""
    ispurine(nt::NucleicAcid)

Test if `nt` is surely a purine.
"""
@inline function ispurine(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b1010) == 0
end

"""
    ispyrimidine(nt::NucleicAcid)

Test if `nt` is surely a pyrimidine.
"""
@inline function ispyrimidine(nt::NucleicAcid)
    bits = reinterpret(UInt8, nt)
    return bits != 0 && (bits & 0b0101) == 0
end

"""
    isambiguous(nt::NucleicAcid)

Test if `nt` is an ambiguous nucleotide.
"""
@inline function isambiguous(nt::NucleicAcid)
    return count_ones(nt) > 1
end

"""
    iscertain(nt::NucleicAcid)

Test if `nt` is a non-ambiguous nucleotide e.g. ACGT.
"""
@inline function iscertain(nt::NucleicAcid)
    return count_ones(nt) == 1
end

"""
    complement(nt::NucleicAcid)

Return the complementary nucleotide of `nt`.

This function returns the union of all possible complementary nucleotides.

Examples
--------

```jldoctest
julia> complement(DNA_A)
DNA_T

julia> complement(DNA_N)
DNA_N

julia> complement(RNA_U)
RNA_A

```
"""
function complement(nt::NucleicAcid)
    bits = compatbits(nt)
    return reinterpret(
        typeof(nt),
        (bits & 0x01) << 3 | (bits & 0x08) >> 3 |
        (bits & 0x02) << 1 | (bits & 0x04) >> 1)
end

function Base.isvalid(::Type{T}, x::Integer) where T <: NucleicAcid
    return 0 ≤ x < 16
end

function Base.isvalid(nt::NucleicAcid)
    return reinterpret(UInt8, nt) ≤ 0b1111
end

"""
    compatbits(nt::NucleicAcid)

Return the compatibility bits of `nt` as `UInt8`.

Examples
--------

```jldoctest
julia> compatbits(DNA_A)
0x01

julia> compatbits(DNA_C)
0x02

julia> compatbits(DNA_N)
0x0f

```
"""
@inline function compatbits(nt::NucleicAcid)
    return reinterpret(UInt8, nt)
end


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

name = "Arpack"
version = v"3.5.0-3"

# Collection of sources required to build ArpackMKLBuilder
sources = [
    "https://github.com/opencollab/arpack-ng.git" =>
    "b095052372aa95d4281a645ee1e367c28255c947",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd arpack*
mkdir build
cd build
# cmake -DINTERFACE64=1 -DBUILD_SHARED_LIBS=ON -DLAPACK_FOUND=true -DLAPACK_INCLUDE_DIRS=${MKLROOT}/include -DLAPACK_LIBRARIES=mkl_rt -DBLAS_FOUND=true -DBLAS_INCLUDE_DIRS=${MKLROOT}/include -DBLAS_LIBRARIES=mkl_rt ..
CMAKE_EXE_LINKER_FLAGS="${CMAKE_EXE_LINKER_FLAGS} -Wl,-rpath -Wl,$MKLROOT/lib"
# Symbols that have float32, float64, complexf32, and complexf64 support
SDCZ_SYMBOLS="axpy copy gemv geqr2 lacpy lahqr lanhs larnv lartg \
              lascl laset scal trevc trmm trsen gbmv gbtrf gbtrs \
              gttrf gttrs pttrf pttrs"
# All symbols that have float32/float64 support (including the SDCZ_SYMBOLS above)
SD_SYMBOLS="${SDCZ_SYMBOLS} dot ger labad laev2 lamch lanst lanv2 \
            lapy2 larf larfg lasr nrm2 orm2r rot steqr swap"
# All symbols that have complexf32/complexf64 support (including the SDCZ_SYMBOLS above)
CZ_SYMBOLS="${SDCZ_SYMBOLS} dotc geru unm2r"
# Add in (s|d)*_64 symbol remappings:
for sym in ${SD_SYMBOLS}; do     SYMBOL_DEFS="${SYMBOL_DEFS} -Ds${sym}=s${sym}_64 -Dd${sym}=d${sym}_64"; done
# Add in (c|z)*_64 symbol remappings:
for sym in ${CZ_SYMBOLS}; do     SYMBOL_DEFS="${SYMBOL_DEFS} -Dc${sym}=c${sym}_64 -Dz${sym}=z${sym}_64"; done
# Add one-off symbol mappings; things that don't fit into any other bucket:
for sym in scnrm2 dznrm2 csscal zdscal dgetrf dgetrs; do     SYMBOL_DEFS="${SYMBOL_DEFS} -D${sym}=${sym}_64"; done
# Set up not only lowercase symbol remappings, but uppercase as well:
SYMBOL_DEFS="${SYMBOL_DEFS} ${SYMBOL_DEFS^^}"
FFLAGS="${FFLAGS} -fdefault-integer-8 ${SYMBOL_DEFS} -ff2c"
if [[ "${target}" == *mingw* ]]; then
    FFLAGS="${FFLAGS} -L${prefix}/bin"
fi
cmake .. \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain \
    -DINTERFACE64=1 \
    -DBUILD_SHARED_LIBS=ON \
    -DBLAS_LIBRARIES="-L$MKLROOT/lib -lmkl_rt" \
    -DLAPACK_LIBRARIES="-L$MKLROOT/lib \
    -lmkl_rt" -DCMAKE_EXE_LINKER_FLAGS="${CMAKE_EXE_LINKER_FLAGS}" \
    -DCMAKE_Fortran_FLAGS="${FFLAGS}"
make clean all
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:glibc),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]
platforms = expand_gcc_versions(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://raw.githubusercontent.com/JuliaComputing/MKL.jl/76ab0ac2b671eb4e199c1b4a31872ccb22898e4f/deps/build_MKL.jl"
]

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

max_stars_repo_name: JuliaLinearAlgebra/ArpackMKLBuilder

max_issues_repo_name: JuliaLinearAlgebra/ArpackMKLBuilder

max_forks_repo_name: arnavs/ArpackMKLBuilder


text:
using Compose
function sierpinski(n)
  if n == 0
    compose(context(), polygon([(1,1), (0,1), (1/2, 0)]));
  else
    t = sierpinski(n - 1);
    compose(	context(),	(context( 1/4, 0, 1/2, 1/2), t),
			 	(context( 0, 1/2, 1/2, 1/2), t),
				(context( 1/2, 1/2, 1/2, 1/2), t));
  end
end

cx1 = compose(sierpinski(1), linewidth(0.2mm),
			fill(nothing), stroke("black"));
img = SVG("sierp1.svg", 10cm, 8.66cm); draw(img,cx1)

cx3 = compose(sierpinski(3), linewidth(0.2mm),
			fill(nothing), stroke("black"));
img = SVG("sierp3.svg", 10cm, 8.66cm); draw(img,cx3)

cx5 = compose(sierpinski(5), linewidth(0.2mm),
		fill(nothing), stroke("black"));
img = SVG("sierp5.svg", 10cm, 8.66cm); draw(img,cx5)


text:
 n=rand(1:10)
 @test diag(full(Toeplitz([-n:n]))) == zeros(n+1)


text:
logo = """

▄▄▄▄▄▄▄▄▄▄▄▄▄▄                 ▄▄▄▄▄▄▄▄▄▄▄▄▄▄
█░░░░░░░░░░░░█  ▄▄▄▄▄▄▄▄▄▄▄▄▄  █░░░░░░░░░░░░█
█▒▒▒▒▒▒▒▒▒▒▒▒█  █░░░░░░░░░░░█  █▒▒▒▒▒▒▒▒▒▒▒▒█
█░░░░░░░░░░░░█  █░░░░░░░░░░░█  █░░░░░░░░░░░░█
█▒▒▒▒▒▒▒▒▒▒▒▒█  █░░░(•••)░░░█  █▒▒▒▒▒▒▒▒▒▒▒▒█
█░░░░░░░░░░░░█━━█░(•••••••)░█━━█░░░░░░░░░░░░█
█▒▒▒▒▒▒▒▒▒▒▒▒█  █(•••••••••)█  █▒▒▒▒▒▒▒▒▒▒▒▒█
█░░░░░░░░░░░░█  █░(•••••••)░█  █░░░░░░░░░░░░█
█▒▒▒▒▒▒▒▒▒▒▒▒█  █░░░(•••)░░░█  █▒▒▒▒▒▒▒▒▒▒▒▒█
█░░░░░░░░░░░░█  ▀▀▀▀▀▀▀▀▀▀▀▀▀  █░░░░░░░░░░░░█
▀▀▀▀▀▀▀▀▀▀▀▀▀▀                 ▀▀▀▀▀▀▀▀▀▀▀▀▀▀

"""

greeting="""
                 GNSTools.jl
                   v0.1.0

                   Contact
         Sergei.Bilardi@colorado.edu
"""


text:
dsge = Model1002()

@testset "Constructing a DSGE-VAR" begin

    # Empty construction with a DSGE model
    m = DSGE.DSGEVAR(dsge, Vector{Symbol}(undef, 0), "ss0")
    @test typeof(m.dsge) == typeof(dsge)
    @test spec(m) == "dsgevar_m1002"
    @test subspec(m) == "ss0"
    @test isempty(m.observables)
    @test isempty(m.shocks)

    # Empty construction but w/some shocks
    m = DSGE.DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss0")
    @test typeof(m.dsge) == typeof(dsge)
    @test spec(m) == "dsgevar_m1002"
    @test subspec(m) == "ss0"
    @test isempty(m.observables)
    @test !isempty(m.shocks)

    # Full construction w/shocks and model
    m = DSGE.DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss1")
    @test typeof(m.dsge) == typeof(dsge)
    @test spec(m) == "dsgevar_m1002"
    @test subspec(m) == "ss1"
    @test !isempty(m.observables)
    @test !isempty(m.shocks)

    # Check errors
    @test_throws ErrorException DSGE.DSGEVAR(dsge, [:not_a_shock], "ss0")
end

@testset "Helper functions for DSGE-VAR" begin

    # Check update
    as_dsge = AnSchorfheide()
    m = DSGE.DSGEVAR(as_dsge, collect(keys(as_dsge.exogenous_shocks)), "ss0")
    DSGE.update!(m; shocks = [:rm_sh])
    @test collect(keys(m.shocks)) == [:rm_sh] && isempty(m.observables) && m.lags == 0 && m.λ == 0.
    DSGE.update!(m; observables = [:obs_gdp])
    @test collect(keys(m.shocks)) == [:rm_sh] && collect(keys(m.observables)) == [:obs_gdp] && m.lags == 0 && m.λ == 0.
    DSGE.update!(m; lags = 1)
    @test collect(keys(m.shocks)) == [:rm_sh] && collect(keys(m.observables)) == [:obs_gdp] && m.lags == 1 && m.λ == 0.
    DSGE.update!(m; λ = 0.5)
    @test collect(keys(m.shocks)) == [:rm_sh] && collect(keys(m.observables)) == [:obs_gdp] && m.lags == 1 && m.λ == 0.5

    # Check errors
    @test_throws ErrorException DSGE.update!(m, dsge; shocks = [:not_a_shock])
    @test_throws ErrorException DSGE.update!(m, dsge; observables = [:not_an_observable])
    @test_throws AssertionError DSGE.update!(m, dsge; lags = -1)

    # Check access and size functions
    m = DSGE.DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss0")
    DSGE.update!(m; observables = [:obs_spread], lags = 1)
    @test collect(keys(DSGE.get_observables(m))) == [:obs_spread]
    @test n_observables(m) == 1
    @test collect(keys(DSGE.get_shocks(m))) == collect(keys(m.dsge.exogenous_shocks))
    @test DSGE.get_λ(m) == 0. == m.λ
    @test DSGE.n_shocks(m) == length(m.dsge.exogenous_shocks)
    @test DSGE.n_lags(m) == 1
    @test DSGE.get_lags(m) == 1
    @test typeof(DSGE.get_dsge(m)) == typeof(dsge)
end

fp = dirname(@__FILE__)
@testset "Pathing for DSGEVAR" begin
    ans_base = joinpath(fp, "save", "output_data", "m1002", "ss10", "dsgevar_m1002", "ss0", "estimate")
    ans_end = ["raw", "work", "tables", "figures", "log"]
    m = DSGE.DSGEVAR(dsge, Vector{Symbol}(undef, 0), "ss0")
    m <= Setting(:saveroot, joinpath(fp, "save"))
    m <= Setting(:dataroot, joinpath(fp, "save/input_data"))
    @test rawpath(m, "estimate") == joinpath(ans_base, ans_end[1])
    @test workpath(m, "estimate") == joinpath(ans_base, ans_end[2])
    @test tablespath(m, "estimate") == joinpath(ans_base, ans_end[3])
    @test figurespath(m, "estimate") == joinpath(ans_base, ans_end[4])
    @test logpath(m, "estimate") == joinpath(ans_base, ans_end[5])
end

@testset "DSGEVAR subspecs" begin
    dsge = Model1002("ss10"; custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true),
                                       :add_NominalWageGrowth =>
                                       Setting(:add_NominalWageGrowth, true)))
    m0  = DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss0")
    m1  = DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss1")
    m10 = DSGEVAR(dsge, collect(keys(dsge.exogenous_shocks)), "ss10")

    @test m0.λ  == 0.
    @test m1.λ  == 0.5
    @test m10.λ == 0.5
    @test m0.lags  == 0
    @test m1.lags  == 4
    @test m10.lags == 4
    @test isempty(m0.observables)
    @test collect(keys(m1.observables))  == [:obs_hours, :obs_gdpdeflator]
    @test collect(keys(m10.observables)) == [:obs_hours, :obs_gdpdeflator,
                                            :laborshare_t, :NominalWageGrowth]
    @test collect(values(m1.observables))  == collect(1:2)
    @test collect(values(m10.observables)) == collect(1:4)
end

nothing


text:
# Generate functions for wrapping and setting a callback
macro callback(ex)
	transform = ex.head == :->
	def = transform ? ex.args[1] : ex                  # Foo(x::T1, y::T2, etc)
	name = string(def.args[1], "Callback")             # FooCallback
	args = def.args[2:end]                             # (x::T1, y::T2, etc)
	values = transform ? ex.args[2].args[2].args : map(argname, args)  # (x, y, etc) or whatever came after the ->

	setter = symbol("Set", name)                       # SetFooCallback
	libsetter = Expr(:quote, symbol("glfw", setter))   # glfwSetFooCallback
	wrapper = gensym(string(name, "Wrapper"))          # FooCallbackWrapper

	window_arg = filter(iswindow, args)                # :(window::Window)
	window_value = map(argname, window_arg)            # :window
	handle_type = map(_ -> :WindowHandle, window_arg)  # :WindowHandle

	wrapper_args = map(win2handle, args)
	wrapper_types = Expr(:tuple, map(argtype, wrapper_args)...)
	window_lookup = map(x -> :($x = Base.cconvert(Window, $(win2handle(x)));), window_value)

	if isempty(window_arg)
		callback_ref = gensym(name)
		declare_callback_ref = :($callback_ref = Ref{Function}())
		callback_ref = :($callback_ref[])
	else
		global _window_callbacks_len += 1
		idx = _window_callbacks_len
		callback_ref = :($(window_value[1]).callbacks[$idx])
		declare_callback_ref = nothing
	end

	ex = quote
		$declare_callback_ref

		# Set the callback function
		function $setter($(window_arg...), callback::Function)
			$callback_ref = callback # Prevent Julia function from being garbage-collected
			cfunptr = cfunction($wrapper, Void, $wrapper_types)
			ccall( ($libsetter, lib), Void, ($(handle_type...), Ptr{Void}), $(window_value...), cfunptr)
		end

		# Unset the callback function
		function $setter($(window_arg...), ::Void)
			ccall( ($libsetter, lib), Void, ($(handle_type...), Ptr{Void}), $(window_value...), C_NULL)
			$callback_ref = undef # Allow former callback function to be garbage-collected
			return nothing
		end

		# Julia callback wrapper that can be passed to `cfunction`
		$wrapper($(wrapper_args...)) = ($(window_lookup...); $callback_ref($(values...)); return nothing)
	end
	esc(ex)
end

# Helper functions for the callback macro
argname(ex) = ex.args[1]
argtype(ex) = ex.args[2]
iswindow(ex) = argtype(ex) == :Window
win2handle(name::Symbol) = symbol(name, "_handle")
win2handle(ex) = iswindow(ex) ? :($(win2handle(argname(ex)))::WindowHandle) : ex
undef(any...) = throw(UndefRefError())

# Size of vector to create during Window construction
_window_callbacks_len = 0


text:
# This section defines functions that need to be overloaded so custom Algorithm works alongside other Algorithm with Baytes.
############################################################################################

#########################################
#!NOTE: New algorithm needs to be a subtype of AbstractAlgorithm, and functions for propose and propose! need to be provided.
"""
$(TYPEDEF)
Abstract super type for Algorithms. New algorithm needs to be subtype of `AbstractAlgorithm`.
"""
abstract type AbstractAlgorithm end

"""
    $(FUNCTIONNAME)
Propose new parameter with given algorithm, keeping objective fixed.

# Examples
```julia
propose(_rng::AbstractRNG, algorithm::AbstractAlgorithm, objective::AbstractObjective)
```

"""
function propose end

"""
    $(FUNCTIONNAME)
Inplace version of [`propose`](@ref).

# Examples
```julia
propose!(_rng::AbstractRNG, algorithm::AbstractAlgorithm, model::AbstractModelWrapper, data, temperature, update)
```

"""
function propose! end

#########################################
#NOTE: These functions are needed so sampler can communicate with eacher other in SMC. Not needed outside SMC.
"""
    $(FUNCTIONNAME)
Substitute result position in algorithm with new result.

# Examples
```julia
result!(algorithm::AbstractAlgorithm, result::AbstractResult)
```

"""
function result! end

"""
    $(FUNCTIONNAME)
Show log density result of algorithm.

# Examples
```julia
get_result(algorithm::AbstractAlgorithm)
```

"""
function get_result end

"""
    $(FUNCTIONNAME)
Get tagged parameter of AbstractAlgorithm..
"""
function get_tagged(algorithm::AbstractAlgorithm)
    return algorithm.tune.tagged
end

#########################################
#!NOTE: New Algorithm needs a constructor so given rng, model, data, chains, temperature, algorithm can be initiated.
"""
$(TYPEDEF)
Abstract type to construct Algorithms. New algorithm needs to define a functor to initiate corresponding `AbstractAlgorithm`.
"""
abstract type AbstractConstructor end

"""
    $(FUNCTIONNAME)
Return tagged symbol of AbstractConstructor.
"""
function get_sym(constructor::AbstractConstructor)
    return constructor.sym
end

#########################################
#!NOTE: New Algorithm needs to return a new model parameter and diagnostics that are a subtype of AbstractDiagnostics.
"""
$(TYPEDEF)
Abstract super type for Algorithm Diagnostics. Once `AbstractAlgorithm` has run, return new model parameter and `AbstractDiagnostics` of `AbstractAlgorithm`.
"""
abstract type AbstractDiagnostics end

"""
    $(FUNCTIONNAME)
Return prediction of AbstractAlgorithm diagnostics.

"""
function get_prediction(diagnostics::AbstractDiagnostics)
    return get_prediction(diagnostics.base)
end

"""
    $(FUNCTIONNAME)
Show results of input type, in this case a vector of `AbstractDiagnostics` of corresponding `AbstractAlgorithm`.

# Examples
```julia
results(diagnosticsᵛ::AbstractVector{AbstractDiagnostics}, algorithm, Ndigits, quantiles)
```

"""
function results end

"""
    $(FUNCTIONNAME)
Show values of Algorihm diagnostics.

# Examples
```julia
generate_showvalues(diagnostics::AbstractDiagnostics)
```

"""
function generate_showvalues end

"""
    $(FUNCTIONNAME)
Infer output type from input. Infer `AbstractDiagnostics` from corresponding `AbstractAlgorithm`

# Examples
```julia
infer(_rng, diagnostics::Type{AbstractDiagnostics}, algorithm, model, data)
```

"""
function infer end

############################################################################################
# Export
export
    AbstractAlgorithm,
    propose,
    propose!,

    AbstractConstructor,

    AbstractDiagnostics,
    infer,
    results


text:
## test/extractions.jl --- check if important properties are preserved

@testset "Model extractions (PDB); residue unsorted" begin
    Trajectory(joinpath(dataloc, fname)) do traj
        fr = read(traj)

        extracted = extractframe(fr)
        top = Topology(fr)

        ## equal number of atoms
        @test length(fr) == length(extracted.at_list)
        ## TRIVIAL: step number for trajectories
        @test step(fr) == extracted.step
        ## all positions are the same and in order
        @test let Q = reinterpret(reshape, SVector{3,Float64}, positions(fr))
            flag = true
            for (ch_, ext_) in zip(Q, extracted.at_pos)
                flag &= (ch_ == ext_)
                flag || break
            end
            flag
        end

        ## check residue base contents and order
        @test let Q = [Residue(top, i-1) for i in 1:count_residues(top)]
            flag = true
            for (ch_, ext_) in zip(Q, extracted.res_list)
                flag &= (Symbol(name(ch_)) == ext_.name)
                flag &= all(∈(atoms(ch_).+1), values(ext_.at_dict))
                flag &= begin
                    l = "is_standard_pdb" in list_properties(ch_)
                    l || return true
                    l &= property(ch_, "is_standard_pdb") == ext_.standard_pdb
                end
                flag || break
            end
            flag
        end

        ## check atom base contents and order
        @test let Q = [Atom(fr, i-1) for i in 1:size(fr)]
            flag = true
            for (ch_, ext_) in zip(Q, extracted.at_list)
                flag &= (Symbol(name(ch_)) == ext_.name)
                flag &= (type(ch_) == ext_.type)
                flag &= (mass(ch_) == ext_.mass)
                flag &= (charge(ch_) == ext_.charge)
                flag || break
            end
            flag
        end
    end
end

@testset "Model extractions (PDB); residue sorted" begin
    unsorted = Trajectory(joinpath(dataloc, fname)) do traj
        fr = read(traj)
        extractframe(fr; sort=false)
    end

    sorted = Trajectory(joinpath(dataloc, fname)) do traj
        fr = read(traj)
        extractframe(fr; sort=true)
    end

    ## check if the res list for the "sorted" extraction is actually sorted
    @test let (; res_list) = sorted
        issorted(res_list;
                 lt=(x,y)->isless(x.standard_pdb, y.standard_pdb), rev=true)
    end

    ## TRIVIAL: check if the order of the atoms and positions are maintained
    @test let (pl1, pl2) = (unsorted.at_pos, sorted.at_pos)
        all(zip(pl1, pl2)) do (p1, p2)
            isapprox(p1, p2; rtol=1e-9)
        end
    end

    @test let (al1, al2) = (unsorted.at_list, sorted.at_list)
        all(zip(al1, al2)) do (a1, a2)
            ## just check the names ig
            a1.name == a2.name
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ ce8179c0-d85e-11ec-1eb8-775c4b10457c
using PlutoPlotly

# ╔═╡ ab3f0662-6bef-4470-8ea4-fe8aad0aba26
using PlutoUI

# ╔═╡ a287d938-99ab-4f0d-a13d-2e5eafd3cc1e
using CSV, DataFrames

# ╔═╡ 4dca7d00-888f-471b-9067-b3358d5f9392
using RDatasets # This is needed to use the `dataset` function below

# ╔═╡ 8c3477ca-8043-41e5-bc0b-6bbbad44e8e1
begin
	using Random
	Random.seed!(42)

	N = 100
	random_x = range(0, stop=1, length=N)
	random_y0 = randn(N) .+ 5
	random_y1 = randn(N)
	random_y2 = randn(N) .- 5

	plot([
	    scatter(x=random_x, y=random_y0, mode="markers", name="markers"),
	    scatter(x=random_x, y=random_y1, mode="lines", name="lines"),
	    scatter(x=random_x, y=random_y2, mode="markers+lines", name="markers+lines")
	])
end

# ╔═╡ f5499780-449c-4e04-a944-385b6f2a4e2f
TableOfContents()

# ╔═╡ bf5c0b14-ff0b-4cfd-a847-3ff586c78ae1
md"""
# Preamble

This is a simple adaptation of the example codes at [https://plotly.com/julia/line-and-scatter/](https://plotly.com/julia/line-and-scatter/) to use `PlutoPlotly.jl` rather than `PlotlyJS.jl`
"""

# ╔═╡ 8a4aa83f-99ca-4fe2-b8fe-1437d741deed
md"""
# Scatter Plots in Julia
"""

# ╔═╡ d732d5e0-4a82-46a7-acfb-0299ab767777
md"""
## Scatter Trace Type
The scatter trace type can be used to represent scatter charts (one point or marker per observation), line charts (a line drawn between each point), or bubble charts (points with size proportional to a dimension of the observation).

To draw a scatter chart, use the scatter trace type and set the mode parameter to markers.
"""

# ╔═╡ 083502fe-7997-46f8-8cb7-a8a213d448c8
# x and y given as arrays
plot(scatter(x=1:10, y=rand(10), mode="markers"))

# ╔═╡ 0daf552b-d434-41ab-88d2-6cb7ff864c61
# x and y given as DataFrame columns
df = dataset("datasets", "iris") # Synthax of `dataset` seems changed from plotly.com example

# ╔═╡ f5bad65a-62c4-4852-b151-8a62b225c30c
plot(scatter(df; x=:SepalWidth, y=:SepalLength, mode="markers")) # Column Names are now CamelCase rather than snake_case

# ╔═╡ a1834f35-9550-4b25-b480-c7a706a20ed8
md"""
### Set size and color with column names
Note that you can set `marker_size` via column name and generate multiple traces using group.
"""

# ╔═╡ 506ecf90-d909-4c33-b535-727f6befc73c
plot(scatter(
    df; x=:SepalWidth, y=:SepalLength, group=:Species,
    marker=attr(size=:PetalLength, sizeref=maximum(df.PetalLength) / (20^2), sizemode="area"),
    mode="markers"
)) # scatter call had to be added here compared to the web documentation, code was also slightly wrong (had `color` instead of `group`)

# ╔═╡ 5cd83e71-fd95-4ce6-8ec9-a766d0bddfd1
md"""
## Line plots
By setting mode to lines, you can draw a line chart.
"""

# ╔═╡ 7f7f1c58-01ef-4928-a2a5-53660b29af14
t = 0:0.01:2π

# ╔═╡ fcc6e511-9c6f-4067-ba98-eeb4f2dc9f37
plot(
    scatter(x=t, y=cos.(t), mode="lines"),
    Layout(yaxis_title="cos(t)", xaxis_title="t")
)

# ╔═╡ e7a0a93e-667a-4d69-b7c2-f0aae0895955
html"""
You can also plot functions directly:

This functionality from the plotly.com page does not seem to actually work
"""

# ╔═╡ 62d6a403-139b-436c-920e-860452ce56ee
# Doesn't do anything, just errors on the javascript side
plot(cos, 0, 2π, mode="lines", Layout(title="cos(t)", template="none"))

# ╔═╡ 03e5a1b6-d690-4321-8dc4-612adf6986ba
md"""
As well as dataframes
"""

# ╔═╡ 3930f946-548a-4688-b026-db77c7f8c61e
let # Put in let to avoid clash with existing df
	df = dataset("datasets", "gapminder") # Can't seem to find this dataset
	df_ocean = df[df.continent .== "Oceania", :]
	plot(
	    df_ocean, x=:year, y=:lifeExp, color=:country, mode="lines"
	)
end

# ╔═╡ f6802ed8-e7bf-4faa-bc30-f4bbc6a97f50
md"""
## Line and Scatter Plots
Use mode argument to choose between markers, lines, or a combination of both.
"""

# ╔═╡ 8141fb07-97ef-4f79-94fb-945f579e35de
md"""
## Bubble Scatter Plots
In [bubble charts](https://en.wikipedia.org/wiki/Bubble_chart), a third dimension of the data is shown through the size of markers. For more examples, see the bubble chart docs
"""

# ╔═╡ b8483605-1ec9-4ddd-87a5-a890c4ae8999
plot(scatter(
    x=1:4, y=10:13, mode="markers", marker=attr(size=40:20:100, color=0:3)
))

# ╔═╡ 50d85c2a-cff1-4be7-8be5-a015120f5628
md"""
## Style Scatter Plots
There are many properties of the scatter trace type that control differetn aspects of the appearance of the trace. Here are a few examples
"""

# ╔═╡ 83226702-c2d9-483c-a05d-3ff6f58ddf05
let
	p = plot(
	    [sin, cos], 0, 10, mode="markers", marker=attr(size=10, line_width=2),
	    Layout(title="Styled Scatter", yaxis_zeroline=false, xaxis_zeroline=false)
	)
	restyle!(p.Plot, 1, marker_color="rgba(152, 0, 0, 0.8)")
	restyle!(p.Plot, 2, marker_color="rgba(255, 182, 193, 0.9)")
	p
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
PlutoPlotly = "8e989ff0-3d88-8e9f-f020-2b208a939ff0"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
CSV = "~0.10.4"
DataFrames = "~1.3.4"
PlutoPlotly = "~0.3.4"
PlutoUI = "~0.7.39"
RDatasets = "~0.7.7"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "873fb188a4b9d76549b81465b1f75c82aaf59238"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.4"

[[deps.CategoricalArrays]]
deps = ["DataAPI", "Future", "Missings", "Printf", "Requires", "Statistics", "Unicode"]
git-tree-sha1 = "109664d3a6f2202b1225478335ea8fea3cd8706b"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.10.5"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.3"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "7297381ccb5df764549818d9a7d57e45f1057d30"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.18.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "a985dc37e357a3b22b260a5def99f3530fb415d3"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.2"

[[deps.ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.8"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.10.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "daa21eb85147f72e41f6352a57fccea377e310a9"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.4"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "cc1a8e22627f33c789ab60b36a9132ac050bbf75"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.12"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.ExprTools]]
git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.8"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "9267e5f50b0e12fdfd5a2455534345c4cf2c7f7a"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.14.0"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "129b104185df66e408edd6625d480b7f9e9823a0"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.18"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
deps = ["Tricks"]
git-tree-sha1 = "c47c5fa4c5308f27ccaac35504858d8914e102f9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.4"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "336cc738f03e069ef2cac55a104eb823455dca75"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.4"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.15"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlotlyBase]]
deps = ["ColorSchemes", "Dates", "DelimitedFiles", "DocStringExtensions", "JSON", "LaTeXStrings", "Logging", "Parameters", "Pkg", "REPL", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "180d744848ba316a3d0fdf4dbd34b77c7242963a"
uuid = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
version = "0.8.18"

[[deps.PlutoPlotly]]
deps = ["AbstractPlutoDingetjes", "Dates", "HypertextLiteral", "InteractiveUtils", "LaTeXStrings", "Markdown", "PlotlyBase", "PlutoUI", "Reexport"]
git-tree-sha1 = "b470931aa2a8112c8b08e66ea096c6c62c60571e"
uuid = "8e989ff0-3d88-8e9f-f020-2b208a939ff0"
version = "0.3.4"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "8d1f54886b9037091edf146b517989fc4a09efec"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.39"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.2"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.RData]]
deps = ["CategoricalArrays", "CodecZlib", "DataFrames", "Dates", "FileIO", "Requires", "TimeZones", "Unicode"]
git-tree-sha1 = "19e47a495dfb7240eb44dc6971d660f7e4244a72"
uuid = "df47a6cb-8c03-5eed-afd8-b6050d6c41da"
version = "0.8.3"

[[deps.RDatasets]]
deps = ["CSV", "CodecZlib", "DataFrames", "FileIO", "Printf", "RData", "Reexport"]
git-tree-sha1 = "2720e6f6afb3e562ccb70a6b62f8f308ff810333"
uuid = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
version = "0.7.7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "bc40f042cfcc56230f781d92db71f0e21496dffd"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.5"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TimeZones]]
deps = ["Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "0a359b0ee27e4fbc90d9b3da1f48ddc6f98a0c9e"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.7.3"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.Tricks]]
git-tree-sha1 = "6bac775f2d42a611cdfcd1fb217ee719630c4175"
uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775"
version = "0.1.6"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═ce8179c0-d85e-11ec-1eb8-775c4b10457c
# ╠═ab3f0662-6bef-4470-8ea4-fe8aad0aba26
# ╠═a287d938-99ab-4f0d-a13d-2e5eafd3cc1e
# ╠═4dca7d00-888f-471b-9067-b3358d5f9392
# ╠═f5499780-449c-4e04-a944-385b6f2a4e2f
# ╟─bf5c0b14-ff0b-4cfd-a847-3ff586c78ae1
# ╠═8a4aa83f-99ca-4fe2-b8fe-1437d741deed
# ╟─d732d5e0-4a82-46a7-acfb-0299ab767777
# ╠═083502fe-7997-46f8-8cb7-a8a213d448c8
# ╠═0daf552b-d434-41ab-88d2-6cb7ff864c61
# ╠═f5bad65a-62c4-4852-b151-8a62b225c30c
# ╟─a1834f35-9550-4b25-b480-c7a706a20ed8
# ╠═506ecf90-d909-4c33-b535-727f6befc73c
# ╟─5cd83e71-fd95-4ce6-8ec9-a766d0bddfd1
# ╠═7f7f1c58-01ef-4928-a2a5-53660b29af14
# ╠═fcc6e511-9c6f-4067-ba98-eeb4f2dc9f37
# ╠═e7a0a93e-667a-4d69-b7c2-f0aae0895955
# ╠═62d6a403-139b-436c-920e-860452ce56ee
# ╟─03e5a1b6-d690-4321-8dc4-612adf6986ba
# ╠═3930f946-548a-4688-b026-db77c7f8c61e
# ╟─f6802ed8-e7bf-4faa-bc30-f4bbc6a97f50
# ╠═8c3477ca-8043-41e5-bc0b-6bbbad44e8e1
# ╟─8141fb07-97ef-4f79-94fb-945f579e35de
# ╠═b8483605-1ec9-4ddd-87a5-a890c4ae8999
# ╟─50d85c2a-cff1-4be7-8be5-a015120f5628
# ╠═83226702-c2d9-483c-a05d-3ff6f58ddf05
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
Module for providing MPI functionality for `Rimu`. This module is unexported. To use it, run

```
using Rimu.RMPI
```
"""
module RMPI

import MPI
using Rimu
using LinearAlgebra
using Rimu.ConsistentRNG
using Random
using StaticArrays

import Rimu: sort_into_targets!

export MPIData
export mpi_rank, is_mpi_root, @mpi_root, mpi_barrier
export mpi_comm, mpi_root, mpi_size, mpi_seed_CRNGs!, mpi_allprintln

function __init__()
    # Initialise the MPI library once at runtime.
    MPI.Initialized() || MPI.Init()
    # make sure that MPI ranks have independent random numbers
    mpi_seed_CRNGs!()
end

const mpi_registry = Dict{Int,Any}()
# The registry keeps references to distributed structures that must not be
# garbage collected at random
abstract type DistributeStrategy end

include("mpidata.jl")
include("helpers.jl")
include("noexchange.jl")
include("pointtopoint.jl")
include("onesided.jl")
include("alltoall.jl")

end # module RMPI


text:
using Plots; pyplot()

years = ["pf-5", "pf-6", "pf-7", "pf-8", "pf-9"]
gdp = [8014, 7957, 7929, 8041, 8048]

# Plot the graph and create a line chart, years on x-axis, gdp on y-axis
plot(years, gdp, marker='o', linestyle=:solid, linecolor=( :green))
# add a title
plot!(title = "Algorithms")
# add a label to the y-axis
plot!(ylabel = "Counts \$")


text:
# This file was generated, do not modify it. # hide
gdf = groupby(iris, :Species);


text:
using Test
using SimpleLife

A = zeros(Int, 5, 5)
@test life_step(A) == A

B = ones(Int,5,5)
@test sum(life_step(B)) == 4


text:
using Oceananigans.Grids: AbstractGrid

#####
##### Momentum flux operators
#####

@inline momentum_flux_huu(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_uu(i, j, k, grid, advection, solution.uh, solution.uh) / solution.h[i, j, k]

@inline momentum_flux_huv(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_uv(i, j, k, grid, advection, solution.uh, solution.vh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_hvu(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_vu(i, j, k, grid, advection, solution.uh, solution.vh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)

@inline momentum_flux_hvv(i, j, k, grid, advection, solution) =
    @inbounds momentum_flux_vv(i, j, k, grid, advection, solution.vh, solution.vh) / solution.h[i, j, k]

#####
##### Momentum flux divergence operators
#####

@inline div_hUu(i, j, k, grid, advection, solution) =
    1 / Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, momentum_flux_huu, advection, solution) +
                               δyᵃᶜᵃ(i, j, k, grid, momentum_flux_huv, advection, solution))

@inline div_hUv(i, j, k, grid, advection, solution) =
    1 / Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_hvu, advection, solution) +
                               δyᵃᶠᵃ(i, j, k, grid, momentum_flux_hvv, advection, solution))

# Support for no advection
@inline div_hUu(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)
@inline div_hUv(i, j, k, grid::AbstractGrid{FT}, ::Nothing, solution) where FT = zero(FT)

#####
##### Mass flux and flux divergence operators
#####

@inline mass_flux_x(i, j, k, grid, uh) = @inbounds Ax_ψᵃᵃᶠ(i, j, k, grid, uh)
@inline mass_flux_y(i, j, k, grid, vh) = @inbounds Ay_ψᵃᵃᶠ(i, j, k, grid, vh)

div_UV(i, j, k, grid, advection, solution) =
    1 / Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, mass_flux_x, advection, solution.uh) +
                               δyᵃᶜᵃ(i, j, k, grid, mass_flux_y, advection, solution.vh))


text:
using Mimi
using OptiMimi
using Test

# Basic model to optimize:
# Each period, you have 1 unit to divide between consuming and staving off disaster.
# If disaster happens, you get 0 consumption for the remainder of the game.
# Utility is defined as sqrt(c).  The probability of disaster is 1-c in each period.
# Assume 5% discounting and 10 periods.

## Analytical solution

## V_t = max_c u(c) + gamma (1 - c) V_t+1
## u'(c) - gamma V_t+1 = 0
## Let u(c) = sqrt(c), u'(c) = .5 c^(-.5)
## c = 1/(2 gamma V_t+1)^2
gamma = exp(-.05) # 5% discount rate

# consume everything in last period
VV_reverse = [1.]
consumption_reverse = [1.]

for tt in 10:-1:1
    consumption = 1 / (2 * gamma * VV_reverse[end])^2
    VV = sqrt(consumption) + gamma * (1 - consumption) * VV_reverse[end]
    push!(consumption_reverse, consumption)
    push!(VV_reverse, VV)
end

println(reverse(consumption_reverse))
println(reverse(VV_reverse))

# OptiMimi solution

# Create simple model
@defcomp Bellmano begin
    # The x-value to evaluate the quadratic
    consumption = Parameter(index=[time])

    # The y-value of the quadratic at the x-value
    utility = Variable(index=[time])
    disaster::Bool = Variable(index=[time])
    bonus = Variable()

    function run_timestep(p, v, d, t)
        v.bonus = 0
        if is_first(tt) || !v.disaster[tt-1]
            v.utility[tt] = sqrt(p.consumption[tt])
            v.disaster[tt] = rand() < p.consumption[tt]
	    if tt == 10 && !v.disaster[tt]
	        v.bonus = 1
	    end
        else
            v.utility[tt] = 0
            v.disaster[tt] = true
        end
    end
end

m = Model()
set_dimension!(m, :time, collect(1:10))

bellmano = add_comp!(m, Bellmano)
bellmano[:consumption] = repmat([.25], 10)

run(m)
m[:Bellmano, :utility]

import OptiMimi.uncertainproblem
using DataFrames

df = DataFrame(algorithm=[], mcperlife=[], cons1=[], cons2=[], cons3=[], cons4=[], cons5=[], cons6=[], cons7=[], cons8=[], cons9=[], cons10=[], cons1serr=[], cons2serr=[], cons3serr=[], cons4serr=[], cons5serr=[], cons6serr=[], cons7serr=[], cons8serr=[], cons9serr=[], cons10serr=[], time=[])
push!(df, [:exact; Inf; reverse(consumption_reverse)[1:10]; repmat([0], 10); 0])

for samplemult in [1, 2, 4]
    for mcperlife in [5, 20, 50]
        println(mcperlife)

        prob = uncertainproblem(m, [:Bellmano], [:consumption], [0.], [1.], m -> sum(m[:Bellmano, :utility] .* exp(-(0:9) * .05)) + m[:Bellmano, :bonus] * exp(-10 * .05), (model) -> nothing)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :social, mcperlife, samplemult)
        time = toc()
        push!(df, [:social; mcperlife; soln.xmean; soln.xserr; time])

        println(df)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :biological, mcperlife, 110000*samplemult)
        time = toc()
        push!(df, [:biological; mcperlife; soln.xmean; soln.xserr; time])

        println(df)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :sampled, mcperlife, 300*samplemult)
        time = toc()
        push!(df, [:sampled; mcperlife; soln.xmean; soln.xserr; time])

        println(df)
    end
    writetable("attempts-compare.csv", df)
end


text:
"""

  FitMethods(f :: Func)

  Generates the methods for each fit function which allow calls using only
  upper or lower bounds to be defined, in any order.

"""
macro FitMethods(f)
  @eval begin
    T = AbstractArray{<:Real}
    # No options provided
    $f(X :: T, Y :: T, l :: lower; s :: Float64, options :: Options) = $f(X,Y,l=l,s=s,options=options)
    $f(X :: T, Y :: T, u :: upper; kargs)= $f(X,Y,u=u,kargs...)
    $f(X :: T, Y :: T, l :: lower, u :: upper; kargs) = $f(X,Y,l=l,u=u,kargs...)
    $f(X :: T, Y :: T, u :: upper, l :: lower; kargs) = $f(X,Y,l=l,u=u,kargs...)
  end
end

macro FitMethodsExponential(f)
  @eval begin
    # No options provided
    $f(X :: T, Y :: T, l :: lower; n :: Int = 1) = $f(X,Y,l=l,n=n)
    $f(X :: T, Y :: T, u :: upper; n :: Int = 1) = $f(X,Y,u=u,n=n)
    $f(X :: T, Y :: T, l :: lower, u :: upper; n :: Int = 1) = $f(X,Y,l=l,u=u,n=n)
    $f(X :: T, Y :: T, u :: upper, l :: lower; n :: Int = 1) = $f(X,Y,l=l,u=u,n=n)
    # Options provided
    $f(X :: T, Y :: T, options :: Options; n :: Int = 1) = $f(X,Y,options=options,n=n)
    $f(X :: T, Y :: T, l :: lower, options :: Options; n :: Int = 1) = $f(X,Y,l=l,options=options,n=n)
    $f(X :: T, Y :: T, u :: upper, options :: Options; n :: Int = 1) = $f(X,Y,u=u,options=options,n=n)
    $f(X :: T, Y :: T, l :: lower, u :: upper, options :: Options; n :: Int = 1) = $f(X,Y,l=l,u=u,options=options,n=n)
    $f(X :: T, Y :: T, u :: upper, l :: lower, options :: Options; n :: Int = 1) = $f(X,Y,l=l,u=u,options=options,n=n)
  end
end


text:
using Distributed
addprocs(8)
@everywhere include("reproducibility.jl")
include("reproducibility.jl")
@everywhere using LightGraphs
@everywhere using LinearAlgebra
@everywhere using QSWalk
@everywhere using SparseArrays
@everywhere using NPZ

data_folder_name = "convergence_qsw"

##

a = [0 0 0 0 1 1 0; 0 0 0 0 1 0 0; 0 0 0 0 1 1 0; 0 0 0 0 1 1 0; 1 1 1 1 0 1 0; 1 0 1 1 1 0 1; 0 0 0 0 0 1 0]
l, vset = nm_lind(a)
h = nm_glob_ham(a)
hrot = nm_loc_ham(vset)
s = sparse(evolve_generator(h, [l], hrot))

init1 = nm_init(vset[[1]], vset)
init2 = nm_init(vlist(vset), vset)
times = 0:10:1000

println("result 1")
result1 = hcat(pmap(t->nm_measurement(evolve(s, init1, t), vset), times)...)
npzwrite("$data_folder_name/special_nmgqsw_1.npz", result1)
println("result 2")
result2 = hcat(pmap(t->nm_measurement(evolve(s, init2, t), vset), times)...)
npzwrite("$data_folder_name/special_nmgqsw_2.npz", result2)


text:
include("./training_set_generator.jl")
using TrainingSetGenerator

#Perceptron realization class
type Perceptron
    inputs_number::Int8
    weights::Array{Float64}
    learning_rate::Float64
end

#Predicting class based on the observation vector
function predict(observation, perceptron)
    sum = dot(observation[1], perceptron.weights)
    return return_sign(sum)
end

#Signum
function return_sign(number)
    if number > 0
        return 1
    else
        return -1
    end
end

#Method to recalculate weights according to the error
function recalculate_weight(observation, perceptron)
        label = observation[2]
        prediction = predict(observation, perceptron)
        error = label - prediction
        #Recalculate weights
        perceptron.weights += error * observation[1] * perceptron.learning_rate
        return perceptron, prediction, label
end

#Calculating error rate of the current epoch
function error_rate(predictions, targets)
    total_positives = 0
    for i = 1:length(targets)
        if targets[i] == predictions[i]
            total_positives +=1
        end
    end
    accuracy = total_positives/length(targets)
    return accuracy
end

#Method to train perceptron
function train(inputs, perceptron)
    predictions = []
    true_labels = []
    for i=1:length(inputs)
        perceptron, prediction, label = recalculate_weight(inputs[i], perceptron)
        push!(predictions, prediction)
        push!(true_labels, label)
    end
    accuracy = error_rate(predictions, true_labels)
    return perceptron, accuracy
end

#Method to test perceptron
function test(inputs, perceptron)
    testing_predictions =  []
    true_labels = []
    for i=1:length(inputs)
        prediction = predict(inputs[i], perceptron)
        push!(testing_predictions, prediction)
        push!(true_labels, inputs[i][2])
    end
    accuracy = error_rate(testing_predictions, true_labels)
    println(string("Accuracy on the testing data is: ", accuracy * 100, "%"))
    return perceptron
end

#Metod to run perceptron
function run_perceptron(epochs)

    #Initializing perceptron with, three inputs(two coordinates and a bias) with random weights
    inputs_number = 3
    bias = 1
    learning_rate = 0.2
    perceptron = Perceptron(inputs_number, [rand(-1:1), rand(-1:1), bias], learning_rate)
    accuracy = 0

    sample_training_set = TrainingSetGenerator.generate_data(200)

    for i= 1:epochs
        perceptron, accuracy = train(sample_training_set, perceptron)
    end

    println(string("Accuracy on the training data is: ", accuracy * 100, "%"))
    testing_data = TrainingSetGenerator.generate_data(20)
    test(testing_data, perceptron)

end

max_stars_repo_name: AlexAdamenko/Perceptron

max_issues_repo_name: AlexAdamenko/Perceptron

max_forks_repo_name: AlexAdamenko/Perceptron