text:
# @file runtests.jl
# @author Ryan Curtin
#
# Tests for the Julia bindings.
using Pkg
Pkg.activate(".")
using Test
using mlpack: test_julia_binding, GaussianKernel, serialize_bin, deserialize_bin
using Serialization
import Base.Filesystem

# The return order for the binding is this:
#
# col_out, double_out, int_out, matrix_and_info_out, matrix_out, model_bw_out,
# model_out, row_out, str_vector_out, string_out, ucol_out, umatrix_out,
# urow_out, vector_out
#
# That's a lot of parameters!  But this is an atypical binding...

# Test that when we run the binding correctly (with correct input parameters),
# we get the expected output.
@testset "TestRunBindingCorrectly" begin
  _, dblOut, intOut, _, _, _, _, _, _, strOut, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         flag1=true)

  @test strOut == "hello2"
  @test intOut == 13
  @test dblOut == 5.0
end

# If we forget the mandatory flag, we should get wrong results.
@testset "TestRunBindingNoFlag" begin
  _, dblOut, intOut, _, _, _, _, _, _, strOut, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello")

  @test strOut != "hello2"
  @test intOut != 13
  @test dblOut != 5.0
end

# If we give the wrong string, we should get wrong results.
@testset "TestRunBindingWrongString" begin
  _, _, _, _, _, _, _, _, _, strOut, _, _, _, _ =
      test_julia_binding(4.0, 12, "goodbye",
                         flag1=true)

  @test strOut != "hello2"
end

# If we give the wrong int, we should get wrong results.
@testset "TestRunBindingWrongInt" begin
  _, _, intOut, _, _, _, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 15, "hello",
                         flag1=true)

  @test intOut != 13
end

# If we give the wrong double, we should get wrong results.
@testset "TestRunBindingWrongDouble" begin
  _, dblOut, _, _, _, _, _, _, _, _, _, _, _, _ =
      test_julia_binding(2.0, 12, "hello",
                         flag1=true)

  @test dblOut != 5.0
end

# If we give the second flag, this should fail.
@testset "TestRunBadFlag" begin
  _, dblOut, intOut, _, _, _, _, _, _, strOut, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         flag1=true,
                         flag2=true)

  @test strOut != "hello2"
  @test intOut != 13
  @test dblOut != 5.0
end

# The matrix we pass in, we should get back with the third dimension doubled and
# the fifth forgotten.
@testset "TestMatrix" begin
  x = rand(100, 5)

  _, _, _, _, matOut, _, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         matrix_in=x,
                         points_are_rows=true)

  @test size(matOut, 1) == 100
  @test size(matOut, 2) == 4
  @test typeof(matOut[1, 1]) == Float64
  for i in [0, 1, 3]
    for j in 1:100
      @test matOut[j, i + 1] == x[j, i + 1]
    end
  end

  for j in 1:100
    @test matOut[j, 3] == 2 * x[j, 3]
  end
end

# The matrix we pass in, we should get back with the third dimension doubled and
# the fifth forgotten.  This test is column major.
@testset "TestMatrixColMajor" begin
  x = rand(5, 100)

  _, _, _, _, matOut, _, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         matrix_in=x,
                         points_are_rows=false)

  @test size(matOut, 1) == 4
  @test size(matOut, 2) == 100
  @test typeof(matOut[1, 1]) == Float64
  for i in 1:100
    for j in [0, 1, 3]
      @test matOut[j + 1, i] == x[j + 1, i]
    end
  end

  for j in 1:100
    @test matOut[3, j] == 2 * x[3, j]
  end
end

# Same as TestMatrix but with an unsigned matrix.
@testset "TestUMatrix" begin
  # Generate a random matrix of integers.
  x = convert(Array{Int, 2}, rand(1:500, (100, 5)))

  _, _, _, _, _, _, _, _, _, _, _, umatOut, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         umatrix_in=x,
                         points_are_rows=true)

  @test size(umatOut, 1) == 100
  @test size(umatOut, 2) == 4
  @test typeof(umatOut[1, 1]) == Int
  for i in [0, 1, 3]
    for j in 1:100
      @test umatOut[j, i + 1] == x[j, i + 1]
    end
  end

  for j in 1;100
    # Since we subtract one when we convert to C++, and then add one when we
    # convert back, we get a slightly different result here.
    @test umatOut[j, 3] == 2 * x[j, 3] - 1
  end
end

# Same as TestMatrix but with an unsigned column major matrix.
@testset "TestUMatrixColMajor" begin
  # Generate a random matrix of integers.
  x = convert(Array{Int, 2}, rand(1:500, (5, 100)))

  _, _, _, _, _, _, _, _, _, _, _, umatOut, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         umatrix_in=x,
                         points_are_rows=false)

  @test size(umatOut, 1) == 4
  @test size(umatOut, 2) == 100
  @test typeof(umatOut[1, 1]) == Int
  for i in 1:100
    for j in [0, 1, 3]
      @test umatOut[j + 1, i] == x[j + 1, i]
    end
  end

  for j in 1;100
    # Since we subtract one when we convert to C++, and then add one when we
    # convert back, we get a slightly different result here.
    @test umatOut[3, j] == 2 * x[3, j] - 1
  end
end

# Test a column vector input parameter.
@testset "TestCol" begin
  x = rand(100)

  colOut, _, _, _, _, _, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         col_in=x)

  @test size(colOut, 1) == 100
  @test typeof(colOut) == Array{Float64, 1}

  for i in 1:100
    @test colOut[i] == 2 * x[i]
  end
end

# Test an unsigned column vector input parameter.
@testset "TestUCol" begin
  x = convert(Array{Int, 1}, rand(1:500, 100))

  _, _, _, _, _, _, _, _, _, _, ucolOut, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         ucol_in=x)

  @test size(ucolOut, 1) == 100
  @test typeof(ucolOut) == Array{Int, 1}
  for i in 1:100
    # Since we subtract one when we convert to C++, and then add one when we
    # convert back, we get a slightly different result here.
    @test ucolOut[i] == 2 * x[i] - 1
  end
end

# Test a row vector input parameter.
@testset "TestRow" begin
  x = rand(100)

  _, _, _, _, _, _, _, rowOut, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         row_in=x)

  @test size(rowOut, 1) == 100
  @test typeof(rowOut) == Array{Float64, 1}
  for i in 1:100
    @test rowOut[i] == 2 * x[i]
  end
end

# Test an unsigned row vector input parameter.
@testset "TestURow" begin
  x = convert(Array{Int, 1}, rand(1:500, 100))

  _, _, _, _, _, _, _, _, _, _, _, _, urowOut, _ =
      test_julia_binding(4.0, 12, "hello",
                         urow_in=x)

  @test size(urowOut, 1) == 100
  @test typeof(urowOut) == Array{Int, 1}
  for i in 1:100
    # Since we subtract one when we convert to C++, and then add one when we
    # convert back, we get a slightly different result here.
    @test urowOut[i] == 2 * x[i] - 1
  end
end

# Test that we can pass a matrix with all numeric features.
@testset "TestMatrixAndInfo" begin
  x = rand(Float64, (10, 100))
  # Dimension information.
  dims = [false, false, false, false, false, false, false, false, false, false]
  z = x

  _, _, _, matrix_and_info_out, _, _, _,  _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         matrix_and_info_in=(dims, z),
                         points_are_rows=false)

  @test size(matrix_and_info_out, 1) == 10
  @test size(matrix_and_info_out, 2) == 100

  for i in 1:100
    for j in 1:10
      @test matrix_and_info_out[j, i] == 2.0 * z[j, i]
    end
  end
end

# Test that we can pass a matrix with all numeric features.
@testset "TestMatrixAndInfoRowMajor" begin
  x = rand(Float64, (100, 10))
  # Dimension information.
  dims = [false, false, false, false, false, false, false, false, false, false]
  z = x

  _, _, _, matrix_and_info_out, _, _, _,  _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         matrix_and_info_in=(dims, z),
                         points_are_rows=true)

  @test size(matrix_and_info_out, 1) == 100
  @test size(matrix_and_info_out, 2) == 10

  for i in 1:100
    for j in 1:10
      @test matrix_and_info_out[i, j] == 2.0 * z[i, j]
    end
  end
end

# Test that we can pass a vector of ints and get back that same vector but with
# the last element removed.
@testset "TestIntVector" begin
  x = [1, 2, 3, 4, 5]

  _, _, _, _, _, _, _, _, _, _, _, _, _, vecOut =
      test_julia_binding(4.0, 12, "hello",
                         vector_in=x)

  @test vecOut == [1, 2, 3, 4]
end

# Test that we can pass a vector of strings and get back that same vector but
# with the last element removed.
@testset "TestStringVector" begin
  x = ["one", "two", "three", "four", "five"]

  _, _, _, _, _, _, _, _, strVecOut, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         str_vector_in=x)

  @test strVecOut == ["one", "two", "three", "four"]
end

# First create a GaussianKernel object, then send it back and make sure we get
# the right double value.
@testset "TestModel" begin
  _, _, _, _, _, _, modelOut, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         build_model=true)

  _, _, _, _, _, bwOut, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         model_in=modelOut)

  @test bwOut == 20.0
end

# Test that we can serialize a model and then use it again.
@testset "TestStreamSerialization" begin
  _, _, _, _, _, _, modelOut, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         build_model=true)

  stream = IOBuffer()
  serialize_bin(stream, modelOut)

  newStream = IOBuffer(copy(stream.data))
  newModel = deserialize_bin(newStream, GaussianKernel)

  _, _, _, _, _, bwOut, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         model_in=newModel)
end

# Test that we can serialize a model as part of a larger tuple.
@testset "TestStreamTupleSerialization" begin
  _, _, _, _, _, _, modelOut, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         build_model=true)

  stream = IOBuffer()
  serialize(stream, (modelOut, 3, 4, 5))

  newStream = IOBuffer(copy(stream.data))
  (newModel, a, b, c) = deserialize(newStream)

  _, _, _, _, _, bwOut, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         model_in=newModel)

  @test a == 3
  @test b == 4
  @test c == 5
end

@testset "TestFileSerialization" begin
  _, _, _, _, _, _, modelOut, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         build_model=true)

  open("model.bin", "w") do io
    serialize_bin(io, modelOut)
  end

  local newModel
  open("model.bin", "r") do io
    newModel = deserialize_bin(io, GaussianKernel)
  end

  _, _, _, _, _, bwOut, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         model_in=newModel)

  Filesystem.rm("model.bin")
end

@testset "TestBaseSerialization" begin
  _, _, _, _, _, _, modelOut, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         build_model=true)

  serialize("model.bin", modelOut)
  newModel = deserialize("model.bin")

  _, _, _, _, _, bwOut, _, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello",
                         model_in=newModel)

  Filesystem.rm("model.bin")
end

# Ensure that we don't accidentally free a model multiple times.
@testset "TestMultipleModelDealloc" begin
  _, _, _, _, _, _, model, _, _, _, _, _, _, _ =
      test_julia_binding(4.0, 12, "hello", build_model=true)

  begin
    for i = 1:100
      out = test_julia_binding(4.0, 12, "hello", model_in=model,
          duplicate_model=true)
    end
  end

  # This should free the other models.  It's likely to crash if a model might be
  # freed multiple times.
  GC.gc()
end

max_stars_repo_name: oblanchet/mlpack

max_issues_repo_name: oblanchet/mlpack

max_forks_repo_name: oblanchet/mlpack


text:
using SparseDiffTools
using LightGraphs
using Test

using Random
Random.seed!(123)

#= Test data =#
test_graphs = Array{SimpleGraph, 1}(undef, 0)

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

#=
 Coloring needs to satisfy two conditions:

1. every pair of adjacent vertices receives distinct colors
(a distance-1 coloring)

2. For any vertex v, any color that leads to a two-colored path
involving v and three other vertices  is  impermissible  for  v.
In other words, every path on four vertices uses at least three
colors.
=#

#Sample graph from Gebremedhin AH, Manne F, Pothen A. **What color is your Jacobian? Graph coloring for computing derivatives.**

#=
     (2)
    /  \
   /    \
 (1)----(3)----(4)

=#

gx = SimpleGraph(4)

add_edge!(gx,1,2)
add_edge!(gx,1,3)
add_edge!(gx,2,3)
add_edge!(gx,3,4)

push!(test_graphs, gx)

#begin testing
for i in 1:6
    g = test_graphs[i]

    out_colors1 = SparseDiffTools.color_graph(g,SparseDiffTools.GreedyStar1Color())
    out_colors2 = SparseDiffTools.color_graph(g,SparseDiffTools.GreedyStar2Color())

    #test condition 1
    for v = vertices(g)
        color = out_colors1[v]
        for j in inneighbors(g, v)
            @test out_colors1[j] != color
        end
    end

    #test condition 2
    for j = vertices(g)
        walk = LightGraphs.self_avoiding_walk(g, j, 4)
        walk_colors = zeros(Int64, 0)
        if length(walk) >= 4
            for t in walk
                push!(walk_colors, out_colors1[t])
            end
            @test length(unique(walk_colors)) >= 3
        end
    end

    #test condition 1
    for v = vertices(g)
        color = out_colors2[v]
        for j in inneighbors(g, v)
            @test out_colors2[j] != color
        end
    end

    #test condition 2
    for j = vertices(g)
        walk = LightGraphs.self_avoiding_walk(g, j, 4)
        walk_colors = zeros(Int64, 0)
        if length(walk) >= 4
            for t in walk
                push!(walk_colors, out_colors2[t])
            end
            @test length(unique(walk_colors)) >= 3
        end
    end

end

ext: jl


text:
# Autogenerated wrapper script for Taskflow_jll for any
JLLWrappers.@generate_wrapper_header("Taskflow")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Observables
import Makie

export Sphere
export RGBSphere
export update

"""
    Represents a sphere.

fields: q, radius, segments, color and observable.
"""
mutable struct Sphere <: Sprite
    q::Biquaternion
    radius::Float64
    segments::Int
    color::Observables.Observable{Array{Makie.ColorTypes.RGBA{Float32},2}}
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

mutable struct RGBSphere <: Sprite
    q::Biquaternion
    radius::Float64
    segments::Int
    color::Any
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Sphere(q, scene)

Construct a sphere with the given configuration `q` and `scene`, and the optional arguments:
`radius`, `segments`, `color` and `transparency`.
"""
function Sphere(q::Biquaternion,
                scene::Makie.Scene;
                radius::Float64 = 1.0,
                segments::Int = 36,
                color::Makie.RGBAf = Makie.RGBAf(1.0, 0.2705, 0.0, 0.5),
                transparency::Bool = false)
    sphere = constructsphere(q, radius, segments = segments)
    colorarray = Observables.Observable(fill(color, segments, segments))
    observable = buildsurface(scene, sphere, colorarray, transparency = transparency)
    Sphere(q, radius, segments, colorarray, observable)
end

function RGBSphere(q::Biquaternion,
                scene::Makie.Scene,
                color::Any;
                radius::Float64 = 1.0,
                segments::Int = 36,
                transparency::Bool = false)
    sphere = constructsphere(q, radius, segments = segments)
    observable = buildsurface(scene, sphere, color, transparency = transparency)
    RGBSphere(q, radius, segments, color, observable)
end

"""
    update(sphere, q)

Update a Sphere by changing its observable with the given `sphere` and configuration `q`.
"""
function update(sphere::Sphere, q::Biquaternion)
    sphere.q = q
    value = constructsphere(sphere.q, sphere.radius, segments = sphere.segments)
    updatesurface(value, sphere.observable)
end

"""
    update(sphere, q)

Update a Sphere by changing its observable with the given `sphere` and configuration `q`.
"""
function update(sphere::RGBSphere, q::Biquaternion)
    sphere.q = q
    value = constructsphere(sphere.q, sphere.radius, segments = sphere.segments)
    updatesurface(value, sphere.observable)
end

"""
    update(sphere, color)

Update a Sphere by changing its observable with the given `sphere` and `color`.
"""
function update(sphere::Sphere, color::Makie.RGBAf)
    sphere.color[] = fill(color, sphere.segments, sphere.segments)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
function _test_interface(::Type{T}) where T<:AbstractGraph
    @testset "$(T)" begin
        g = T(5)
        @testset "basics" begin
            @test countedges(g) == 0
            @test countvertices(g) == 5
            edges = [(1, 2, 3), (1, 5, 4), (2, 4, 2), (4, 3, 7), (5, 2, 1)]
            map(edges) do (i, j, w)
                setedge!(g, i, j, w)
            end
            @test countedges(g) == 5
            map(edges) do (i, j, w)
                @test isedge(g, i, j)
                @test getweight(g, i, j) == w
            end
            @test !isedge(g, 3, 4)
            @test first(g, 1) == 2
            @test first(g, 2) == 4
            @test first(g, 3) == 6
            @test first(g, 4) == 3
            @test first(g, 5) == 2
            @test next(g, 1, first(g, 1)) == 5
            @test next(g, 2, first(g, 2)) == 6
            @test next(g, 4, first(g, 4)) == 6
            @test next(g, 5, first(g, 5)) == 6
            @test_throws ArgumentError setedge!(g, 1, 4, -1)
            map(edges) do (i, j, w)
                @test getweight(g, i, j) == w
            end
            map(1:5) do v
                @test getmark(g, v) == 0
            end
            setmark!(g, 1, 1)
            @test getmark(g, 1) == 1
        end
        @testset "remove!" begin
            remove!(g, 1, 5)
            @test !isedge(g, 1, 5)
            remove!(g, 1, 2)
            @test !isedge(g, 1, 2)
            @test getweight(g, 1, 5) == getweight(g, 1, 2) == 0
            @test countedges(g) == 3
        end
    end
end


text:
"""
MPS is a series of tensors, and also a vector in Hilbert space.
"""
mutable struct MPS{FloatType<:Number,ArrayType<:AbstractArray{FloatType,3}} <: AbstractMPS
    L::Int #length
    S::Int #physical index
    bdim::Bdims #bond dimension
    tensor::Array{ArrayType,1}
end

function MPS(L::Int,S::Int,D::Int;FloatType=Float64,ArrayType=Array)
    bdim = Bdims(L,append!([D for i=1:L-1],1))
    tensor = [zeros(FloatType,bdim[i-1],S,bdim[i]) for i=1:L]
    return MPS{FloatType,ArrayType{FloatType,3}}(L,S,bdim,tensor)
end

"""
overlap of two mps: mps1 and mps2.
the same as transpose(mps2)*mps1

Automatically using CUDA
"""
overlap(mps1::AbstractMPS,mps2::AbstractMPS;ArraType=CuArray) = _overlap(mps1,mps2,ArraType)

"""
    Compress a MPS to bond dimension: Dcut
"""
function compress!(mps::AbstractMPS,Dcut::Int,epsilon=1E-13;use_cuda=false)
    normalization!(mps,LeftNormalization())
    res = cutoff!(mps,Dcut,LeftNormalization())
    return res
end

function dag!(mps::MPS{T,V}) where T<:Complex where V
    for i=1:1:mps.L
        mps[i] = conj.(mps[i])
    end
    return mps
end

function dag(mps::MPS{T,V}) where T<:Complex where V
    temp = deepcopy(mps)
    dag!(temp)
    return temp
end

dag(mps::MPS{T,V}) where T<: Real where V = mps
dag!(mps::MPS{T,V}) where T<: Real where V = mps

"""
    display mps(array...)
    Just think that it display an element for a vector in Hilbert space.
"""
function disp(mps::AbstractMPS,array::Array{Int,1})
    @assert mps.L == size(array)[1]
    t = transpose(mps[1][1,array[1],:])
    for i = 2:mps.L
        t = t*mps[i][:,array[i],:]
    end
    return t[1,1]
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
module SeasonalTrendLoess

using Statistics

export stl

include("loess.jl")
include("stl.jl")

end # module


text:
#!/usr/bin/env julia

# Note!  This requires GuassianMixtures, which is not in REQURE.

using FITSIO
using GaussianMixtures
using Distributions

function mag_to_nanomaggies(mag::Float64)
	10 ^ (( mag - 22.5 ) / -2.5)
end

@vectorize_1arg Float64 mag_to_nanomaggies

function read_r_colors(prior_file)
	fits = FITS(prior_file)
  cat = fits[2]
	num_rows, = read_key(cat, "NAXIS2")
	table = Array(Float64, num_rows, 12)
	for i in 1:12
    col_name, = read_key(cat, "TTYPE$i")
		table[:, i] = read(cat, col_name)
	end
	close(fits)

	t2 = table[:,:]
	for err_col in 8:12
		t2 = t2[t2[:,err_col] .< median(table[:,err_col]) / 2, :]
	end

	for b in 2:7
		t2 = t2[t2[:, b] .> -50, :]
	end

	t3 = mag_to_nanomaggies(t2[:, 3:7])
	colors = Array(Float64, size(t2)[1], 4)
	for i in 1:4
		colors[:, i] = log(t3[:, i + 1] ./ t3[:, i])
	end

	#TODO: convert 'r' band to nanomaggies
	colors, mag_to_nanomaggies(table[table[:, 5] .> 1, 5])
end

function vecmat_to_tensor(vecmat::Vector{Matrix{Float64}})
    ret = Array(Float64, size(vecmat[1], 1), size(vecmat[1], 2), length(vecmat))
    for i in 1:length(vecmat)
        ret[:, :, i] = vecmat[i]
    end
    ret
end

if length(ARGS) != 2
    println("usage: gen_priors.jl [catalog.fits] [out_file.dat]")
else
    c0, r0 = read_r_colors(ARGS[1])

    fit_r = fit_mle(LogNormal, r0)

    D = 64
    c0_train = c0
    #c0_test = c0[120001:end, :]
    fit_gmm = GMM(D, c0_train, kind=:full, method=:split)
    println("train avll:", GaussianMixtures.avll(fit_gmm, c0_train))
    #println("test avll:", GaussianMixtures.avll(fit_gmm, c0_test))

    out_file = open(ARGS[2], "w+")
    serialize(out_file, (
        params(fit_r),
        weights(fit_gmm),
        means(fit_gmm)',
        vecmat_to_tensor(covars(fit_gmm))))
    close(out_file)
end


text:
# -------------------- #
# Dispatcher functions #
# -------------------- #

const BASE_TYPES = [:spli, :cheb, :lin]
const ABSR_MAP = Dict(
    :none => Direct(),
    :direct => Direct(),
    :tensor => Tensor(),
    :expanded => Expanded(),
)
get_bformat(b::T) where T<: BasisMatrix{Direct} = :direct
get_bformat(b::T) where T<:BasisMatrix{Expanded} = :expanded
get_bformat(b::T) where T<:BasisMatrix{Tensor} = :tensor

function to_dict(bm::BasisMatrix)
    B = Dict{Symbol, Any}()
    B[:order] = bm.order
    B[:format] = get_bformat(bm)
    B[:vals] = bm.vals
    B
end

function bm_from_dict(B::Dict)
    arr_type = eltype(B[:vals])
    bm = BasisMatrix{typeof(ABSR_MAP[B[:format]]),arr_type}(B[:order], B[:vals])
    bm
end

base_exists(s::Symbol) = s in BASE_TYPES

basedef(s::Symbol, args...) =
    s == :spli ? splidef(args...) :
    s == :cheb ? chebdef(args...) :
    s == :lin  ? lindef(args...)  :
    error("somehow you snuck through here you 👺")

basenode(s::Symbol, args...) =
    s == :spli ? splinode(args...) :
    s == :cheb ? chebnode(args...) :
    s == :lin  ? linnode(args...)  :
    error("somehow you snuck through here you 👺")

BasisMatrices.evalbase(s::Symbol, args...) =
    s == :spli ? splibase(args...) :
    s == :cheb ? chebbase(args...) :
    s == :lin  ? linbase(args...)  :
    error("somehow you snuck through here you 👺")

# Helper function
function squeeze_trail(x::AbstractArray)
    sz = size(x)
    squeezers = Int[]
    n = length(sz)
    for i=n:-1:1
        if sz[i] == 1
            push!(squeezers, i)
        else
            break
        end
    end
    squeeze(x, tuple(squeezers...))
end

# ---------------------------- #
# Generic translated functions #
# ---------------------------- #

# from fundef.m  -- DONE
function fundef(foo...)
    d = length(foo)  # 89
    n = zeros(Int, d)  # 93
    b = zeros(d)  # 94
    a = zeros(d)  # 95
    p = Array{Any}(undef, d)  # 96
    _params = Array{BasisMatrices.BasisParams}(undef, d)

    basetype = Array{Symbol}(undef, d)
    for j=1:d
        basetype[j] = foo[j][1]  # 99
        !(base_exists(basetype[j])) && error("Unknown basis $(foo[j][1])")
        n[j], a[j], b[j], p[j], _params[j] = basedef(basetype[j], foo[j][2:end]...)  # 124
    end

    # package output. Lines 143-150
    g = Dict{Symbol, Any}()
    g[:d] = d
    g[:n] = n
    g[:a] = a
    g[:b] = b
    g[:basetype] = basetype
    g[:params] = p
    g[:_basis_params] = _params
    g[:_basis] = Basis(_params...)
    g
end

# fundefn.m
function fundefn(basistype::Symbol, n, a, b, order=3)
    d = length(n)
    length(a) != d && error("a must be same dimension as n")
    length(b) != d && error("b must be same dimension as n")
    any(a .> b) && error("left endpoints must be less than right endpoints")
    any(n .< 2) && error("n(i) must be greater than 1")

    params = Array{Any}(undef, 1, d)
    if basistype == :cheb
        for i=1:d params[i] = Any[:cheb, n[i], a[i], b[i]] end
    elseif basistype == :spli
        for i=1:d params[i] = Any[:spli, [a[i], b[i]], n[i]-order+1, order] end
    elseif basistype == :lin
        for i=1:d params[i] = Any[:lin, [a[i], b[i]], n[i]] end
    end

    fundef(params...)
end

# funnode.m -- DONE
funnode(basis::Dict) = nodes(basis[:_basis])

# funbase.m -- DONE
function funbase(basis::Dict, x=funnode(basis)[1], order=fill(0, 1, basis[:d]))
    BasisMatrix(basis[:_basis], Expanded(), x, order).vals[1]
end

# funbasex.m -- DONE
function funbasex(basis::Dict{Symbol}, x=funnode(basis)[1], order=0,
                  bformat::Symbol=:none)
    to_dict(BasisMatrix(basis[:_basis], ABSR_MAP[bformat], x, order))
end

# funfitf.m -- DONE
funfitf(basis, f::Function, args...) = funfitf(basis[:_basis], f, args...)

# funfitxy.m -- DONE
function funfitxy(basis, x, y)
    c, bm = funfitxy(basis[:_basis], x, y)
    return c, to_dict(bm)
end

# funeval.m -- DONE
function funeval(c, basis::Dict, B, _order=0)
    isempty(c) && error("missing basis coefficients")
    order = BasisMatrices._check_order(basis[:d], _order)

    if isa(B, Dict)  # B is a basis structure
        bm = bm_from_dict(B)
        y = funeval(c, bm, order)
        return y, B
    else
        bm = BasisMatrix(basis[:_basis], B, order)
        y = funeval(c, bm, bm.order)
        return y, to_dict(bm)
    end
end

# fund.m
function fund(c, basis, x, hess_opt)
    # TODO: come back when I need this. I think I should probably do something
    #       like what optim does and write functions `f`, `fg!` and `fgh!` to
    #       replicate this for the type of basis instead of this function
    nothing
end

# funbconv.m  -- DONE
function funbconv(b::Dict, order=fill(0, 1, size(b[:order], 2)),
                  format::Symbol=:expanded)
    bm = bm_from_dict(b)
    new_bm = convert(typeof(ABSR_MAP[format]), bm, order)
    to_dict(new_bm)
end


text:
using SparseArrays
using LinearAlgebra
using MatrixDepot
using BenchmarkTools
using UnicodePlots
using PrettyTables
using ChainPartitioners
using Cthulhu
using Profile

for mtx in [
            #"DIMACS10/chesapeake",
            #"HB/can_292",
            "Boeing/ct20stif",
            #"Schmid/thermal1",
            #"Rothberg/3dtube",
           ]
    A = permutedims(1.0 * sparse(mdopen(mtx).A))
    (m, n) = size(A)

    y = zeros(m)
    x = rand(n)
    ref_time = @belapsed(mul!($y, $A, $x))

    println(mtx)
    rows = []
    for K = [4, 8, 16]
        mdl = ConstrainedCost(AffineConnectivityModel{Int64}(0, 0, 0, 1), VertexCount(), ceil(Int, n/K * 1.5))

        (j_lo, j′_hi) = ChainPartitioners.column_constraints(A, K, mdl.w, mdl.w_max)
        show(histogram(j′_hi .- j_lo, nbins = 15, closed = :left))
        println()

        #println()
        #println()
        #println(spy(A, maxwidth=50, maxheight=50, title="$mtx"))
        for (key, method) in [
            ("dynamic", DynamicTotalSplitter(mdl)),
            ("dynamic′", DynamicTotalChunker(mdl)),
            ("quadrangle", ConvexTotalSplitter(mdl)),

        ]
            setup_time = @belapsed(partition_stripe($A, $K, $method))/ref_time
            Φ = nothing
            #@profile begin
                Φ = partition_stripe(A, K, method)
            #end
            #Profile.print()
            comm = total_value(A, Φ, mdl)
            @assert issorted(Φ.spl)
            @assert Φ.spl[1] == 1
            @assert Φ.spl[end] == n + 1

            push!(rows, [key K setup_time comm])
        end
    end
    pretty_table(vcat(rows...), ["method", "K", "setuptime", "communication"])
end


text:
export
    moon_apogee_perigee_time_low,
    moon_illuminated_fraction_high,
    moon_illuminated_fraction_low,
    lunation,
    moon_mean_ascending_node_longitude,
    moon_mean_perigee_longitude,
    moon_altitude,
    moon_constants,
    moon_dimension,
    moon_dimension3,
    moon_horizontal_parallax,
    moon_latitude,
    moon_longitude,
    moon_nearest_phase,
    moon_radius,
    moon_riseset,
    moon_age_location,
    moon_node,
    moon_true_ascending_node_longitude,
    moon_rst_altitude

#=
Lunar position model ELP2000-82 of Chapront.

The resulting values are for the equinox of date and have been adjusted
for light-time.

This is the simplified version of Jean Meeus, _Astronomical Algorithms_,
second edition, 1998, Willmann-Bell, Inc.

=#

# [Meeus-1998: table 47.A]
#
#    D, M, M1, F, l, r

moon_table_a = Array[
    [0,  0,  1,  0, 6288774, -20905355],
    [2,  0, -1,  0, 1274027,  -3699111],
    [2,  0,  0,  0,  658314,  -2955968],
    [0,  0,  2,  0,  213618,   -569925],
    [0,  1,  0,  0, -185116,     48888],
    [0,  0,  0,  2, -114332,     -3149],
    [2,  0, -2,  0,   58793,    246158],
    [2, -1, -1,  0,   57066,   -152138],
    [2,  0,  1,  0,   53322,   -170733],
    [2, -1,  0,  0,   45758,   -204586],
    [0,  1, -1,  0,  -40923,   -129620],
    [1,  0,  0,  0,  -34720,    108743],
    [0,  1,  1,  0,  -30383,    104755],
    [2,  0,  0, -2,   15327,     10321],
    [0,  0,  1,  2,  -12528,        0],
    [0,  0,  1, -2,   10980,     79661],
    [4,  0, -1,  0,   10675,    -34782],
    [0,  0,  3,  0,   10034,    -23210],
    [4,  0, -2,  0,    8548,    -21636],
    [2,  1, -1,  0,   -7888,     24208],
    [2,  1,  0,  0,   -6766,     30824],
    [1,  0, -1,  0,   -5163,     -8379],
    [1,  1,  0,  0,    4987,    -16675],
    [2, -1,  1,  0,    4036,    -12831],
    [2,  0,  2,  0,    3994,    -10445],
    [4,  0,  0,  0,    3861,    -11650],
    [2,  0, -3,  0,    3665,     14403],
    [0,  1, -2,  0,   -2689,     -7003],
    [2,  0, -1,  2,   -2602,        0],
    [2, -1, -2,  0,    2390,     10056],
    [1,  0,  1,  0,   -2348,      6322],
    [2, -2,  0,  0,    2236,     -9884],
    [0,  1,  2,  0,   -2120,      5751],
    [0,  2,  0,  0,   -2069,        0],
    [2, -2, -1,  0,    2048,     -4950],
    [2,  0,  1, -2,   -1773,      4130],
    [2,  0,  0,  2,   -1595,        0],
    [4, -1, -1,  0,    1215,     -3958],
    [0,  0,  2,  2,   -1110,        0],
    [3,  0, -1,  0,    -892,      3258],
    [2,  1,  1,  0,    -810,      2616],
    [4, -1, -2,  0,     759,     -1897],
    [0,  2, -1,  0,    -713,     -2117],
    [2,  2, -1,  0,    -700,      2354],
    [2,  1, -2,  0,     691,        0],
    [2, -1,  0, -2,     596,        0],
    [4,  0,  1,  0,     549,     -1423],
    [0,  0,  4,  0,     537,     -1117],
    [4, -1,  0,  0,     520,     -1571],
    [1,  0, -2,  0,    -487,     -1739],
    [2,  1,  0, -2,    -399,        0],
    [0,  0,  2, -2,    -381,     -4421],
    [1,  1,  1,  0,     351,        0],
    [3,  0, -2,  0,    -340,        0],
    [4,  0, -3,  0,     330,        0],
    [2, -1,  2,  0,     327,        0],
    [0,  2,  1,  0,    -323,      1165],
    [1,  1, -1,  0,     299,        0],
    [2,  0,  3,  0,     294,        0],
    [2,  0, -1, -2,       0,      8752]]

# [Meeus-1998: table 47.B]
#
#    D, M, M1, F, b

moon_table_b = Array[
    [0,  0,  0,  1, 5128122],
    [0,  0,  1,  1,  280602],
    [0,  0,  1, -1,  277693],
    [2,  0,  0, -1,  173237],
    [2,  0, -1,  1,   55413],
    [2,  0, -1, -1,   46271],
    [2,  0,  0,  1,   32573],
    [0,  0,  2,  1,   17198],
    [2,  0,  1, -1,    9266],
    [0,  0,  2, -1,    8822],
    [2, -1,  0, -1,    8216],
    [2,  0, -2, -1,    4324],
    [2,  0,  1,  1,    4200],
    [2,  1,  0, -1,   -3359],
    [2, -1, -1,  1,    2463],
    [2, -1,  0,  1,    2211],
    [2, -1, -1, -1,    2065],
    [0,  1, -1, -1,   -1870],
    [4,  0, -1, -1,    1828],
    [0,  1,  0,  1,   -1794],
    [0,  0,  0,  3,   -1749],
    [0,  1, -1,  1,   -1565],
    [1,  0,  0,  1,   -1491],
    [0,  1,  1,  1,   -1475],
    [0,  1,  1, -1,   -1410],
    [0,  1,  0, -1,   -1344],
    [1,  0,  0, -1,   -1335],
    [0,  0,  3,  1,    1107],
    [4,  0,  0, -1,    1021],
    [4,  0, -1,  1,     833],
    [0,  0,  1, -3,     777],
    [4,  0, -2,  1,     671],
    [2,  0,  0, -3,     607],
    [2,  0,  2, -1,     596],
    [2, -1,  1, -1,     491],
    [2,  0, -2,  1,    -451],
    [0,  0,  3, -1,     439],
    [2,  0,  2,  1,     422],
    [2,  0, -3, -1,     421],
    [2,  1, -1,  1,    -366],
    [2,  1,  0,  1,    -351],
    [4,  0,  0,  1,     331],
    [2, -1,  1,  1,     315],
    [2, -2,  0, -1,     302],
    [0,  0,  1,  3,    -283],
    [2,  1,  1, -1,    -229],
    [1,  1,  0, -1,     223],
    [1,  1,  0,  1,     223],
    [0,  1, -2, -1,    -220],
    [2,  1, -1, -1,    -220],
    [1,  0,  1,  1,    -185],
    [2, -1, -2, -1,     181],
    [0,  1,  2,  1,    -177],
    [4,  0, -2, -1,     176],
    [4, -1, -1, -1,     166],
    [1,  0,  1, -1,    -164],
    [4,  0,  1, -1,     132],
    [1,  0, -1, -1,    -119],
    [4, -1,  0, -1,     115],
    [2, -2,  0,  1,     107]]

# Meeus 47.6

function eccent_E(T)
    E = polynomial([1.0, -0.002516, -0.0000074], T)
    E2 = E*E
    return (E, E2)
end

# Calculate values required by several other functions

function moon_constants(T)
    # Constant terms.
    _kL1 = [deg2rad(218.3164477), deg2rad(481267.88123421), deg2rad(-0.0015786), deg2rad( 1.0/  538841), deg2rad(-1.0/ 65194000)]
    _kD  = [deg2rad(297.8501921), deg2rad(445267.1114034),  deg2rad(-0.0018819), deg2rad( 1.0/  545868), deg2rad(-1.0/113065000)]
    _kM  = [deg2rad(357.5291092), deg2rad( 35999.0502909),  deg2rad(-0.0001536), deg2rad( 1.0/24490000)]
    _kM1 = [deg2rad(134.9633964), deg2rad(477198.8675055),  deg2rad( 0.0087414), deg2rad( 1.0/   69699), deg2rad(-1.0/ 14712000)]
    _kF  = [deg2rad( 93.2720950), deg2rad(483202.0175233),  deg2rad(-0.0036539), deg2rad(-1.0/ 3526000), deg2rad( 1.0/863310000)]

    _kA1 = [deg2rad(119.75), deg2rad(   131.849)]
    _kA2 = [deg2rad( 53.09), deg2rad(479264.290)]
    _kA3 = [deg2rad(313.45), deg2rad(481266.484)]

    L1 = mod2pi(polynomial(_kL1, T))
    D  = mod2pi(polynomial(_kD,  T))
    M  = mod2pi(polynomial(_kM,  T))
    M1 = mod2pi(polynomial(_kM1, T))
    F  = mod2pi(polynomial(_kF,  T))

    A1 = mod2pi(polynomial(_kA1, T))
    A2 = mod2pi(polynomial(_kA2, T))
    A3 = mod2pi(polynomial(_kA3, T))

    (E, E2) = eccent_E(T)

    return (L1, D, M, M1, F, A1, A2, A3, E, E2)
end

# ELP2000 lunar position calculations

"""
    moon_dimension3(jd)

Return geocentric ecliptic longitude, latitude and radius.

        When we need all three dimensions it is more efficient to combine the
        calculations in one routine.

        Parameters:
            jd : Julian Day in dynamical time

        Returns:
            geocentric longitude in radians
            geocentric latitude in radians
            radius in km, Earth's center to Moon's center

"""
function moon_dimension3(jd)
    T = jd_to_jcent(jd)
    L1, D, M, M1, F, A1, A2, A3, E, E2 = moon_constants(T)

    # longitude and radius

    lsum = 0.0
    rsum = 0.0
    for v in moon_table_a
        (tD, tM, tM1, tF, tl, tr) = v
        arg = tD * D + tM * M + tM1 * M1 + tF * F
        if abs(tM) == 1
            tl = tl*E
            tr = tr*E
        elseif abs(tM) == 2
            tl = tl*E2
            tr = tr*E2
        end
        lsum = lsum + tl * sin(arg)
        rsum = rsum + tr * cos(arg)
    end
    #
    # latitude
    #
    bsum = 0.0
    for v in moon_table_b
        (tD, tM, tM1, tF, tb) = v
        arg = tD * D + tM * M + tM1 * M1 + tF * F
        if abs(tM) == 1
         tb = tb*E
        elseif abs(tM) == 2
         tb = tb*E2
        end
        bsum = bsum + tb * sin(arg)
    end

    lsum = lsum + 3958 * sin(A1) +
            1962 * sin(L1 - F) +
             318 * sin(A2)

    bsum = bsum - 2235 * sin(L1) +
            382 * sin(A3) +
            175 * sin(A1 - F) +
            175 * sin(A1 + F) +
            127 * sin(L1 - M1) -
            115 * sin(L1 + M1)

    long = L1 + deg2rad(lsum / 1000000)
    lat = deg2rad(bsum / 1000000)
    dist = 385000.56 + rsum / 1000
    return (long, lat, dist)
end

"""
    moon_longitude(jd)

Return the geocentric ecliptic longitude in radians.

A subset of the logic in moon_dimension3()

"""
function moon_longitude(jd)
    T = jd_to_jcent(jd)
    L1, D, M, M1, F, A1, A2, A3, E, E2 = moon_constants(T)
    lsum = 0.0
    for v in moon_table_a
        (tD, tM, tM1, tF, tl, tr) = v
        arg = tD * D + tM * M +tM1 * M1 + tF * F
        if abs(tM) == 1
            tl = tl*E
        elseif abs(tM) == 2
            tl = tl*E2
        end
        lsum = lsum + tl * sin(arg)
    end

    lsum = lsum + 3958 * sin(A1) +
                  1962 * sin(L1 - F) +
                   318 * sin(A2)

    longitude = L1 + deg2rad(lsum / 1000000)
    return longitude
end

"""
    moon_latitude(jd)

Return the geocentric ecliptic latitude in radians.
 A subset of the logic in moon_dimension3()

"""
function moon_latitude(jd)
    T = jd_to_jcent(jd)
    L1, D, M, M1, F, A1, A2, A3, E, E2 = moon_constants(T)

    bsum = 0.0
    for v in moon_table_b
        tD, tM, tM1, tF, tb = v
        arg = tD * D + tM * M + tM1 * M1 + tF * F
        if abs(tM) == 1
            tb = tb*E
        elseif abs(tM) == 2
            tb = tb*E2
        end
        bsum = bsum + tb * sin(arg)
    end

    bsum = bsum - 2235 * sin(L1) +
                   382 * sin(A3) +
                   175 * sin(A1 - F) +
                   175 * sin(A1 + F) +
                   127 * sin(L1 - M1) -
                   115 * sin(L1 + M1)

    lat = deg2rad(bsum / 1000000)
    return lat
end

"""
    moon_radius(jd)

Return the geocentric radius in km.

A subset of the logic in moon_dimension3()
"""
function moon_radius(jd)
    T = jd_to_jcent(jd)
    (L1, D, M, M1, F, A1, A2, A3, E, E2) = moon_constants(T)

    rsum = 0.0
    for v in moon_table_a
        (tD, tM, tM1,tF, tl, tr) = v
        arg = tD * D + tM * M + tM1 * M1 + tF * F
        if abs(tM) == 1
          tr = tr*E
        elseif abs(tM) == 2
          tr = tr*E2
        end
        rsum = rsum + tr * cos(arg)
    end
    dist = 385000.56 + rsum / 1000
    return dist
end

"""
    moon_dimension(jd, dim)

Return one of geocentric ecliptic longitude, latitude and radius.

Parameters:
    jd : Julian Day in dynamical time
    dim : "L" (longitude") or "B" (latitude) or "R" (radius)

Returns:
    longitude in radians or
    latitude in radians or
    radius in km, Earth's center to Moon's center
"""
function moon_dimension(jd, dim)
    if dim == "L"
        return moon_longitude(jd)
    end
    if dim == "B"
        return moon_latitude(jd)
    end
    if dim == "R"
        return moon_radius(jd)
    end
    error("unknown dimension = $dim")
end

"""
    moon_age_location(jd)

Simple moon phase calculator adapted from the basic program found at
http://www.skyandtelescope.com/resources/software/3304911.html.
This function helps anyone who needs to know the Moon's
phase (age), distance, and position along the ecliptic on
any date within several thousand years in the past or future.
To illustrate its application, Bradley Schaefer applied it
to a number of famous events influenced by the Moon in
World War II.  His article appeared in Sky & Telescope for
April 1994, page 86.

Parameters:
    jd : Julian Day
Returns :
    moon's age from new in days
    distance from  anomalistic phase
    ecliptic latitude
    ecliptic longitude
"""
function moon_age_location(jd)
    global earth_equ_radius
    v  = (jd - 2451550.1) / 29.530588853
    ip = v- floor(Integer, v)
    ag = ip * 29.530588853   # Moon's age from new moon in days
    ip = ip * (pi * 2)       # Converts to radian

    # Calculate distance from anomalistic phase
    v= (jd - 2451562.2) / 27.55454988
    dp = v- floor(Integer, v)
    dp = dp * (pi * 2)
    di = 60.4-3.3 * cos(dp)-.6 * cos(2 * ip - dp) - .5 * cos(2 * ip)

    # Calculate ecliptic latitude from nodal (draconic) phase
    v = (jd - 2451565.2) /27.212220817
    np = v- floor(Integer, v)
    np = np * (pi * 2)
    la = 5.1*sin(np)

    # Calculate ecliptic longitude from sidereal motion
    v=(jd - 2451555.8) / 27.321582241
    rp = v -  floor(Integer, v)
    lo = 360 * rp + 6.3 * sin(dp) + 1.3 * sin(2 * ip - dp) + .7 * sin(2 * ip)

    return (ag, di * earth_equ_radius, deg2rad(la), deg2rad(lo))
end

"""
    moon_nearest_phase(jd, phase=0)

Find the Julian Date of the nearest moon phase (new, 1st quarter, full, last quarter) to the given Julian date.

Parameters:
    jd : the required Julian Day in dynamical time
    phase: integer between 0 and 3: 0-> new moon, 1-> first quarter, 2-> full moon, 3-> second quarter
      defaulting to 0, new moon
Returns :
    the Julian Date for the type of moon phase asked

[Meeus - chapter 49]
"""
function moon_nearest_phase(jd, phase=0)
     _kMean =   [2451550.09766,  29.530588861*1236.85,  0.00015437, -0.000000150, 0.00000000073]
     _kM    =   [      2.55340,  29.105356700*1236.85, -0.0000014,  -0.00000011]
     _kM1   =   [    201.56430, 385.816935280*1236.85,  0.0107582,   0.00001238, -0.000000058]
     _kF    =   [    160.7108,  390.670502848*1236.85, -0.0016118,  -0.00000227, -0.000000011]
     _kO    =   [    124.7746,   -1.563755880*1236.85,  0.0020672,  -0.00000215]
     _cor   =   Dict(
               0.000325 => Array([299.77,  0.107408*1236.85, -0.009173]),
               0.000165 => Array([251.88,  0.016321*1236.85]),
               0.000164 => Array([251.83, 26.651886*1236.85]),
               0.000126 => Array([349.42, 36.412478*1236.85]),
               0.000110 => Array([ 86.66, 18.206239*1236.85]),
               0.000062 => Array([141.74, 53.303771*1236.85]),
               0.000060 => Array([207.14,  2.453732*1236.85]),
               0.000056 => Array([154.84,  7.306860*1236.85]),
               0.000047 => Array([ 34.52, 27.261239*1236.85]),
               0.000042 => Array([207.19,  0.121824*1236.85]),
               0.000040 => Array([291.34,  1.844379*1236.85]),
               0.000037 => Array([161.72, 24.198154*1236.85]),
               0.000035 => Array([239.56, 25.513099*1236.85]),
               0.000023 => Array([331.55,  3.592518*1236.85])
               )
    (yr, mo, d) = jd_to_cal(jd)
    n = cal_to_day_of_year(yr, mo, d)

    if is_leap_year(yr)
        f = n/366 else f = n/365
    end

    k = (yr+f-2000)*12.3685
    # Round to the nearest multiple of 0.25 according to the needed phase
    k = round(k-0.25*phase)+0.25*phase
    T = k/1236.85
    mean_phase = polynomial(_kMean, T)

    (E, E2) = eccent_E(T)
    # M, M1, F, O, new moon, full moon
    _nf_cor = Array[
               [ 0, 1, 0, 0, -0.40720,    -0.40614],
               [ 1, 0, 0, 0,  0.17241*E,   0.17302 * E],
               [ 0, 2, 0, 0,  0.01608,     0.01614],
               [ 0, 0, 2, 0,  0.01039,     0.01043],
               [-1, 1, 0, 0,  0.00739*E,   0.00734 * E],
               [ 1, 1, 0, 0, -0.00514*E,  -0.00515 * E],
               [ 2, 0, 0, 0,  0.00208*E2,  0.00209 * E2],
               [ 0, 1,-2, 0, -0.00111,    -0.00111],
               [ 0, 1, 2, 0, -0.00057,    -0.00057],
               [ 1, 2, 0, 0,  0.00056*E,   0.00056 * E],
               [ 0, 3, 0, 0, -0.00042,    -0.00042],
               [ 1, 0, 2, 0,  0.00042*E,   0.00042 * E],
               [ 1, 0,-2, 0,  0.00038*E,   0.00038 * E],
               [-1, 2, 0, 0, -0.00024*E,  -0.00024 * E],
               [ 0, 0, 0, 1, -0.00017,    -0.00017],
               [ 2, 1, 0, 0, -0.00007,    -0.00007],
               [ 0, 2,-2, 0,  0.00004,     0.00004],
               [ 3, 0, 0, 0,  0.00004,     0.00004],
               [ 1, 1,-2, 0,  0.00003,     0.00003],
               [ 0, 2, 2, 0,  0.00003,     0.00003],
               [ 1, 1, 2, 0, -0.00003,    -0.00003],
               [-1, 1, 2, 0,  0.00003,     0.00003],
               [-1, 1,-2, 0, -0.00002,    -0.00002],
               [ 1, 3, 0, 0, -0.00002,    -0.00002],
               [ 0, 4, 0, 0,  0.00002,     0.00002]
           ]

    # M, M1, F, O, first and last quarters
    _qq_cor = Array[
               [ 0, 1, 0, 0, -0.62801],
               [ 1, 0, 0, 0,  0.17172 * E],
               [ 1, 1, 0, 0, -0.01183 * E],
               [ 0, 2, 0, 0,  0.00862],
               [ 0, 0, 2, 0,  0.00804],
               [-1, 1, 0, 0,  0.00454 * E],
               [ 2, 0, 0, 0,  0.00204 * E2],
               [ 0, 1,-2, 0, -0.00180],
               [ 0, 1, 2, 0, -0.00070],
               [ 0, 3, 0, 0, -0.00040],
               [-1, 2, 0, 0, -0.00034 * E],
               [ 1, 0, 2, 0,  0.00032 * E],
               [ 1, 0,-2, 0,  0.00032 * E],
               [ 2, 1, 0, 0, -0.00028 * E2],
               [ 1, 2, 0, 0,  0.00027 * E],
               [ 0, 0, 0, 1, -0.00017],
               [-1, 1,-2, 0, -0.00005],
               [ 0, 2, 2, 0,  0.00004],
               [ 1, 1, 2, 0, -0.00004],
               [-2, 1, 0, 0,  0.00004],
               [ 1, 1,-2, 0,  0.00003],
               [ 3, 0, 0, 0,  0.00003],
               [ 0, 2,-2, 0,  0.00002],
               [-1, 1, 2, 0,  0.00002],
               [ 1, 3, 0, 0, -0.00002]
            ]

    M  = mod2pi(deg2rad(polynomial(_kM,  T)))
    M1 = mod2pi(deg2rad(polynomial(_kM1, T)))
    F  = mod2pi(deg2rad(polynomial(_kF, T)))
    O  = mod2pi(deg2rad(polynomial(_kO, T)))

    if phase == 0
    # New moon
        for (kM, kM1, kF, kO, coefn) in _nf_cor
           s = sin(kM*M+kM1*M1+kF*F+kO*O)
           mean_phase = mean_phase + s*coefn
      end
    end

    if phase == 2
    # Full moon
        for (kM, kM1, kF, kO, coefn, coeff) in _nf_cor
           s = sin(kM*M+kM1*M1+kF*F+kO*O)
           mean_phase = mean_phase + s*coeff
        end
    end

    if phase == 1 || phase == 3
    # first or last quarters
        for (kM, kM1, kF, kO, coef) in _qq_cor
           s = sin(kM*M+kM1*M1+kF*F+kO*O)
           mean_phase = mean_phase + s*coef
        end

        W =  0.00306
                - 0.00038*E*cos(M)
                + 0.00026*cos(M1)
                - 0.00002*cos(M1-M)
                + 0.00002*cos(M1+M)
                + 0.00002*cos(2*F)
        if phase == 3
            W = -W
        end
        mean_phase = mean_phase +W
    end

    for (k, poly) in _cor
        mean_phase = mean_phase + k * sin(mod2pi(deg2rad(polynomial(poly, T))))
    end
    return mean_phase
end

"""
    moon_illuminated_fraction_high(jd)

Compute the moon illuminated fraction and position's angle of the moon's bright limb

Parameters:
            jd : Julian Day in dynamical time

Returns :
            f, a (fraction, angle)

[Meeus - chapter 48 or 49 depending on edition]
"""
function moon_illuminated_fraction_high(jd)
    global km_per_au
    (lambda, bbeta, delta) = moon_dimension3(jd)  # Moon's geocentric longitude, latitude and radius
    (lambda0, bbeta0, R)   = sun_dimension3(jd)   # Sun's geocentric longitude, latitude and radius

    R = R * km_per_au
    psi = acos(cos(bbeta) * cos(lambda - lambda0))
    i = atan(R * sin(psi), delta - R * cos(psi))
    k = (1+cos(i))/2

    o = true_obliquity(jd)
    (asc,  decl)  = ecl_to_equ(lambda, bbeta, o)
    (asc0, decl0) = ecl_to_equ(lambda0, bbeta0, o)
    khi = atan(cos(decl0) * sin(asc0 - asc), sin(decl0) * cos(decl) - cos(decl0) * sin(decl) * cos(asc0 - asc))
    return (k, mod2pi(khi))
end

"""
    moon_illuminated_fraction_low(jd)

Same as above - lower precision. Simpler algorithm. Does not return the angle of the bright limb
"""
function moon_illuminated_fraction_low(jd)
    T = jd_to_jcent(jd)
    (L1, D, M, M1) = moon_constants(T)

    i = (180 - rad2deg(D)
             - 6.289*sin(M1)
             + 2.100 * sin(M)
             - 1.274 * sin(2 * D - M1)
             - 0.658 * sin(2 * D)
             - 0.214 * sin(2 * M1)
             - 0.110 * sin(D))
    k = (1 + cos(deg2rad(i)))/2
    return (round(k, digits=2))
end
"""
    moon_mean_ascending_node_longitude(jd)

Returns the longitude of the mean ascending node and of the mean perigee of the moon

Parameters:
    jd: Julian day (dynamical time)

Returns
    mean ascending node or mean perigee longitude in radians
"""
function moon_mean_ascending_node_longitude(jd)
    _O = [125.0445479, -1934.1362891, 0.0020754, 1/467441, -1/60616000]
    T = jd_to_jcent(jd)
    return mod2pi(deg2rad(polynomial(_O, T)))
end

"""
    moon_true_ascending_node_longitude(jd)

"""
function moon_true_ascending_node_longitude(jd)
    m = moon_mean_ascending_node_longitude(jd)
    T = jd_to_jcent(jd)
    (L1, D, M, M1, F) = moon_constants(T)
    return mod2pi(
          - deg2rad(1.4979)*sin(2*(D-F))
          - deg2rad(0.1500)*sin(M)
          - deg2rad(0.1226)*sin(2*D)
          + deg2rad(0.1176)*sin(2*F)
          - deg2rad(0.0801)*sin(2*(M1-F))
          + m)
end

"""
    moon_mean_perigee_longitude(jd)

"""
function moon_mean_perigee_longitude(jd)
    _P = [83.3532465, 4096.0137287, -0.01032, -1/80053, -1/189999000]
    T = jd_to_jcent(jd)
    return mod2pi(deg2rad(polynomial(_P, T)))
end

"""
    moon_apogee_perigee_time_low(jd, apo_nperi)

Returns the jd for apogee of perigee, and the corresponding equatorial horizontal parallax

Parameters:
    jd: Julian day (dynamical time)
    apo_nperi : 1 for an apogee, 0 for a perigee

Returns
    jd for apogee or perigee
    parallax in seconds
"""
function moon_apogee_perigee_time_low(jd, apo_nperi)
    (yr, mo, d) = jd_to_cal(jd)
    n = cal_to_day_of_year(yr, mo, d)
    if is_leap_year(yr)
        f = n/366
    else
        f = n/365
    end
    k = (yr + f - 1999.97) * 13.2555
    # Round to the nearest multiple of 0.5 according to apogee or perigee
    k = round(k - 0.5 * apo_nperi) + 0.5 * apo_nperi

    T = k/1325.55
    # k = T*1325.55
    _mT = [ 2451534.6698,  27.55454989*1325.55, -0.0006691, -0.0000001098, 0.0000000052]
    _mD = [ 171.9179,     335.91060460*1325.55, -0.0100383, -0.00001156,   0.000000055]
    _mM = [ 347.3477,      27.15777210*1325.55, -0.0008130, -0.0000010 ]
    _mF = [ 316.6109,     364.52879110*1325.55, -0.0125053, -0.0000148 ]
    # D, M, F, coef - for apogee
    _mA = Array[
         [ 2, 0, 0,  0.4392 ],
         [ 4, 0, 0,  0.0684 ],
         [ 0, 1, 0,  0.0456-0.00011*T ],
         [ 2,-1, 0,  0.0426-0.00011*T ],
         [ 0, 0, 2,  0.0212 ],
         [ 1, 0, 0, -0.0189 ],
         [ 6, 0, 0,  0.0144 ],
         [ 4,-1, 0,  0.0113 ],
         [ 2, 0, 2,  0.0047 ],
         [ 1, 1, 0,  0.0036 ],
         [ 8, 0, 0,  0.0035 ],
         [ 6,-1, 0,  0.0034 ],
         [ 2, 0,-2, -0.0034 ],
         [ 2,-2, 0,  0.0022 ],
         [ 3, 0, 0, -0.0017 ],
         [ 4, 0, 2,  0.0013 ],
         [ 8,-1, 0,  0.0011 ],
         [ 4,-2, 0,  0.0010 ],
         [10, 0, 0,  0.0009 ],
         [ 3, 1, 0,  0.0007 ],
         [ 0, 2, 0,  0.0006 ],
         [ 2, 1, 0,  0.0005 ],
         [ 2, 2, 0,  0.0005 ],
         [ 6, 0, 2,  0.0004 ],
         [ 6,-2, 0,  0.0004 ],
         [10,-1, 0,  0.0004 ],
         [ 5, 0, 0, -0.0004 ],
         [ 4, 0,-2, -0.0004 ],
         [ 0, 1, 2,  0.0003 ],
         [12, 0, 0,  0.0003 ],
         [ 2,-1, 2,  0.0003 ],
         [ 1,-1, 0, -0.0003 ]
         ]
    # D, M, F, coef - for perigee
    _mP = Array[
           [ 2, 0, 0, -1.6769 ],
           [ 4, 0, 0,  0.4589 ],
           [ 6, 0, 0, -0.1856 ],
           [ 8, 0, 0,  0.0883 ],
           [ 2,-1, 0, -0.0773 + 0.00019*T ],
           [ 0, 1, 0,  0.0502 - 0.00013*T ],
           [10, 0, 0, -0.0460 ],
           [ 4,-1, 0,  0.0422 - 0.00011*T ],
           [ 6,-1, 0, -0.0256 ],
           [12, 0, 0,  0.0253 ],
           [ 1, 0, 0,  0.0237 ],
           [ 8,-1, 0,  0.0162 ],
           [14, 0, 0, -0.0145 ],
           [ 0, 0, 2,  0.0129 ],
           [ 3, 0, 0, -0.0112 ],
           [10,-1, 0, -0.0104 ],
           [16, 0, 0,  0.0086 ],
           [12,-1, 0,  0.0069 ],
           [ 5, 0, 0,  0.0066 ],
           [ 2, 0, 2, -0.0053 ],
           [18, 0, 0, -0.0052 ],
           [14,-1, 0, -0.0046 ],
           [ 7, 0, 0, -0.0041 ],
           [ 2, 1, 0,  0.0040 ],
           [20, 0, 0,  0.0032 ],
           [ 1, 1, 0, -0.0032 ],
           [16,-1, 0,  0.0031 ],
           [ 4, 1, 0, -0.0029 ],
           [ 9, 0, 0,  0.0027 ],
           [ 4, 0, 2,  0.0027 ],
           [ 2,-2, 0, -0.0027 ],
           [ 4,-2, 0,  0.0024 ],
           [ 6,-2, 0, -0.0021 ],
           [22, 0, 0, -0.0021 ],
           [18,-1, 0, -0.0021 ],
           [ 6, 1, 0,  0.0019 ],
           [11, 0, 0, -0.0018 ],
           [ 8, 1, 0, -0.0014 ],
           [ 4, 0,-2, -0.0014 ],
           [ 6, 0 ,2, -0.0014 ],
           [ 3, 1, 0,  0.0014 ],
           [ 5, 1, 0, -0.0014 ],
           [13, 0, 0,  0.0013 ],
           [20,-1, 0,  0.0013 ],
           [ 3, 2, 0,  0.0011 ],
           [ 4,-2, 2, -0.0011 ],
           [ 1, 2, 0, -0.0010 ],
           [22,-1, 0, -0.0009 ],
           [ 0, 0, 4, -0.0008 ],
           [ 6, 0,-2,  0.0008 ],
           [ 2, 1,-2,  0.0008 ],
           [ 0, 2, 0,  0.0007 ],
           [ 0,-1, 2,  0.0007 ],
           [ 2, 0, 4,  0.0007 ],
           [ 0,-2, 2, -0.0006 ],
           [ 2, 2,-2, -0.0006 ],
           [24, 0, 0,  0.0006 ],
           [ 4, 0,-4,  0.0005 ],
           [ 2, 2, 0,  0.0005 ],
           [ 1,-1, 0, -0.0004 ]
        ]

    _pA =  Array[
           [ 2, 0, 0, -9.147 ],
           [ 1, 0, 0, -0.841 ],
           [ 0, 0, 2,  0.697 ],
           [ 0, 1, 0, -0.656 + 0.0016*T ],
           [ 4, 0, 0,  0.355 ],
           [ 2,-1, 0,  0.159 ],
           [ 1, 1, 0,  0.127 ],
           [ 4,-1, 0,  0.065 ],
           [ 6, 0, 0,  0.052 ],
           [ 2, 1, 0,  0.043 ],
           [ 2, 0, 2,  0.031 ],
           [ 2, 0,-2, -0.023 ],
           [ 2,-2, 0,  0.022 ],
           [ 2, 2, 0,  0.019 ],
           [ 0, 2, 0, -0.016 ],
           [ 6,-1, 0,  0.014 ],
           [ 8, 0, 0,  0.010 ]]

    _pP =  Array[
           [ 2, 0, 0, 63.224 ],
           [ 4, 0, 0, -6.990 ],
           [ 2,-1, 0,  2.834-0.0071*T ],
           [ 6, 0, 0,  1.927 ],
           [ 1, 0, 0, -1.263 ],
           [ 8, 0, 0, -0.702 ],
           [ 0, 1, 0,  0.696-0.0017*T ],
           [ 0, 0, 2, -0,690 ],
           [ 4,-1, 0, -0.629+0.0016*T ],
           [ 2, 0, 2, -0.392 ],
           [10, 0, 0,  0.297 ],
           [ 6,-1, 0,  0.260 ],
           [ 3, 0, 0,  0.201 ],
           [ 2, 1, 0, -0.161 ],
           [ 1, 1, 0,  0.157 ],
           [12, 0, 0, -0.138 ],
           [ 8,-1, 0, -0.127 ],
           [ 2, 0,-2,  0.104 ],
           [ 5, 0, 0, -0.079 ],
           [14, 0, 0,  0.068 ],
           [10,-1, 0,  0.067 ],
           [ 4, 1, 0,  0.054 ],
           [12,-1, 0, -0.038 ],
           [ 4,-2, 0, -0.038 ],
           [ 7, 0, 0,  0.037 ],
           [ 4, 0, 2, -0.037 ],
           [16, 0, 0, -0.035 ],
           [ 3, 1, 0, -0.030 ],
           [ 1,-1, 0,  0.029 ],
           [ 6, 1, 0, -0.025 ],
           [ 0, 2, 0,  0.023 ],
           [14,-1, 0,  0.023 ],
           [ 2, 2, 0, -0.023 ],
           [ 6,-2, 0,  0.022 ],
           [ 2,-1,-2, -0.021 ],
           [ 9, 0, 0, -0.020 ],
           [18, 0, 0,  0.019 ],
           [ 6, 0, 2,  0.017 ],
           [ 0,-1, 2,  0.014 ],
           [16,-1, 0, -0.014 ],
           [ 4, 0,-2,  0.013 ],
           [ 8, 1, 0,  0.012 ],
           [11, 0, 0,  0.011 ],
           [ 5, 1, 0,  0.010 ],
           [20, 0, 0, -0.010 ]]

    mean_pa = polynomial(_mT, T)
    D = polynomial(_mD, T)
    M = polynomial(_mM, T)
    F = polynomial(_mF, T)

    if apo_nperi == 1
        _m = _mA
        _p = _pA
        parallax = 3245.251 # In arcseconds: 3245".251 according to Meeus table 50.B
    else
        _m = _mP
        _p = _pP
        parallax = 3629.215  # 3629".215 according to Meeus table 50.B
    end

    c = 0
    for (aD, aM, aF, coef) in _m
        c = c+coef*sin(deg2rad(aD*D+aM*M+aF*F))
    end
    mean_pa = mean_pa + c

    for (pD, pM, pF, coef) in _p
        parallax = parallax + coef * cos(deg2rad(pD * D + pM * M + pF * F))
    end

    return (round(mean_pa, digits=4), round(parallax, digits=3))
end

"""
    moon_horizontal_parallax(jd)

"""
function moon_horizontal_parallax(jd)
    global earth_equ_radius
    R = moon_dimension(jd, "R")
    return asin(deg2rad(earth_equ_radius/R))
end

"""
    moon_node(jd, desc_not_asc)

Passage of the moon through the nodes

Parameter
    - jd: julian date
    - desc_not_asc: 0 for the ascending node, 1 for the descending node
Returns
    - julian day of the closer passage through the node
"""
function moon_node(jd, desc_not_asc)
    (yr, mo, d) = jd_to_cal(jd)
    n = cal_to_day_of_year(yr, mo, d)
    if is_leap_year(yr)
        f = n/366
    else
        f = n/365
    end
    k = (yr + f - 2000.05) * 13.4223
    # Round to the nearest multiple of 0.5 according to ascending or descending node
    k = round(k - 0.5 * desc_not_asc) + 0.5 * desc_not_asc
    T = k/1342.33
    # k = T*1342.33 ?
    _mD  = [183.6380, 331.73735682*1342.33, 0.0014852, 0.00000209, -0.00000001  ]
    _mM  = [ 17.4006,  26.82037250*1342.33, 0.0001186, 0.00000006              ]
    _mM1 = [ 38.3776, 355.52747313*1342.33, 0.0123499, 0.000014627,-0.000000069 ]
    _mO  = [123.9767,  -1.44098956*1342.33, 0.0020608, 0.00000214, -0.000000016 ]
    _mV  = [299.7500, 132.85,             -0.009173                         ]
    D =  polynomial(_mD, T)
    M =  polynomial(_mM, T)
    M1 = polynomial(_mM1, T)
    O =  polynomial(_mO, T)
    V =  polynomial(_mV, T)
    P =  O + 272.75 - 2.3 * T

    E, E2 = eccent_E(T)

    _mN = Array[
        [ 0, 0, 1, -0.4721  ],
        [ 2, 0, 0, -0.1649  ],
        [ 2, 0,-1, -0.0868  ],
        [ 2, 0, 1,  0.0084  ],
        [ 2,-1, 0, -0.0083*E],
        [ 2,-1,-1, -0.0039*E],
        [ 0, 0 ,2,  0.0034  ],
        [ 2, 0,-2, -0.0031  ],
        [ 2, 1, 0,  0.0030*E],
        [ 0, 1,-1,  0.0028*E],
        [ 0, 1, 0,  0.0026*E],
        [ 4, 0, 0,  0.0025  ],
        [ 1, 0, 0,  0.0024  ],
        [ 0, 1, 1,  0.0022*E],
        [ 4, 0,-1,  0.0014  ],
        [ 2, 1,-1,  0.0005*E],
        [ 2,-1, 1,  0.0004*E],
        [ 2,-2, 0, -0.0003*E],
        [ 4,-1, 0,  0.0003*E]
    ]

    node = polynomial([2451565.1619, 27.212220817 * 1342.33, 0.0002762, 0.000000021, -0.000000000088], T )
    node = node + 0.0017 * sin(deg2rad(O)) + 0.0003 * sin(deg2rad(V)) + 0.0003 * sin(deg2rad(P))

    for (nD, nM, nM1, coef) in _mN
        c = coef * sin(deg2rad(nD * D + nM * M+ nM1 * M1))
        node = node + c
    end
    return node
end

"""
    lunation(jd, system=false)

Returns the lunation number

Parameters
   - JD julian date (DT or UTC)
   - System can be "brown" (from the first 1923 new moon), "meeus" ( from the first 2000 new moon)
      "islamic" (from the beginning of islamic calendar) or "goldstein" (Goldstein system)
Returns
   - The lunation number
"""
function lunation(jd, system=false)
    if system == false
        system = "brown"
    end
    # moon_nearest_phase will not always return the closest full moon if we are close to a new moon.
    # So we take the next and previous full moon, as well as the new moon in between, and determine
    # the lunation from these three values
    jdf1 = moon_nearest_phase(jd, 2) # Full moon preceding or following jd
    if jdf1 <= jd
        # jdf1 is the previous full moon
        jdf2 = moon_nearest_phase(jdf1 + 29, 2) # jdf2 is the next full moon
    else
        jdf2 = moon_nearest_phase(jdf1 - 29, 2)
        (jdf1, jdf2) = (jdf2, jdf1)
    end
    # now jdf1 is the previous full moon, and jdf2 the next. Find the new moon between jdf1 and jdf2
    jdn = moon_nearest_phase(jdf1 + 14, 0)
    # Compute n in the brown system
    if jd < jdn
        n = ceil((jdf1 - 2423436.6120014)/29.53058853)
    else
        n = ceil((jdf2 - 2423436.6120014)/29.53058853)
    end

    if system == "brown"
        return n
    end
    if system == "meeus"
        return n - 953
    end
    if system == "islamic"
        return n + 17038
    end
    if system == "goldstein"
        return n + 37105
    end
    error("unknown lunation system")
end

"""
    moon_altitude(jd)

Returns the moon's altitude at time JD

Parameters
   - JD julian date

Returns
   - altitude in radians
"""
function moon_altitude(jd)
    global longitude
    (l, b, r)   = moon_dimension3(jd) # geocentric longitude in radians, geocentric latitude in radians, radius in km
    o           = true_obliquity(jd)
    (ra, de)    = ecl_to_equ(l, b, o)
    H           = apparent_sidereal_time_greenwich(jd) - longitude - ra
    (A, h)      = equ_to_horiz(H, de) # azi/alt
    # h = h + 0.7275 * moon_horizontal_parallax(jd) - deg2rad(34/60) -- Semi-diameter of the moon
    # TODO moon_altitude() - add effect of atmospheric refraction
    return h
end

"""
    moon_riseset(jd)

Rising and setting time of the moon for julian day jd

TODO Doesn't work.

# TODO moon_riseset() get this function working
# TODO moon_riseset() why does this return k

"""
function moon_riseset(jd)
    (moonrise, moonset, k, nsteps, uroot) = (false, false, 0, 24, false)
    (jd_list, v_list) = ([0, 1, 2], [0.0, 0.0, 0.0])

    for (i, v) in enumerate(jd_list)
        v_list[i] = moon_altitude(jd + v/nsteps)
    end

    for i in 1:nsteps
        (a, b, c) = quadratic_interpolation(jd_list, v_list)
        s1, s2 = quadratic_roots(a, b, c)

        if s2 == false || s1 == false
            uroot = true
        end # Very unlikely to ever happen due to the precision used

        if s1 >= 0 && s1 <= 2
            slope = a * s1 + b
            j = jd + (i - 1 + s1) / nsteps

            # Convert Julian Day from dynamical to terrestrial universal time
            ut = dt_to_ut(j)

            #println("j $j $(jd_to_date(j)) ut $ut \n slope $slope $(jd_to_date(ut))")

            if slope > 0
                moonrise = ut
            else
                moonset = ut
            end
            if uroot
                k = i
                break
            end
        end

        if (s2 >= 0) && (s2 <= 2)
            slope = a * s2 + b
            j = jd + (i - 1 + s2) / nsteps
            ut = dt_to_ut(j)
            if slope > 0
                moonrise = ut
            else
                moonset = ut
            end
        end

        if moonrise != false && moonset != false  # TODO is this right?
            k = i
            break
        end

        v_list[1], v_list[2] = v_list[2], v_list[3]
        v_list[3] = moon_altitude(jd + (i + 2) / nsteps)
    end
    if uroot
        if moonrise != false
            moonset = moonrise
        else
            moonrise = moonset
        end
    end
    return (moonrise, moonset, k)
end

"""
    moon_rst_altitude(r)

Return the standard altitude of the Moon.

Parameters:
    r : Distance between the centers of the Earth and Moon, in km.

Returns:
    Standard altitude in radians.
"""
function moon_rst_altitude(r)
    global earth_equ_radius, standard_rst_altitude
    # horizontal parallax
    parallax = asin(earth_equ_radius / r)
    return 0.7275 * parallax + standard_rst_altitude
end


text:
mutable struct Catcher{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Ball}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    ball_reward::T
    ball_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(Catcher)
@generate_setters(Catcher)

function Catcher(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, BALL)
    world = GridWorldBase(objects, height, width)

    ball_pos = CartesianIndex(1, 1)
    world[BALL, ball_pos] = true

    agent_pos = CartesianIndex(height, 1)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    terminal_reward = -one(T)
    ball_reward = one(T)
    done = false

    env = Catcher(world, agent_pos, reward, rng, terminal_reward, ball_reward, ball_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.StateStyle(env::Catcher) = RLBase.InternalState{Any}()

RLBase.state_space(env::Catcher, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const CATCHER_LAYERS = SA.SVector(2)
RLBase.state(env::Catcher, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), CATCHER_LAYERS)

RLBase.state_space(env::Catcher, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::Catcher, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::Catcher, player::RLBase.DefaultPlayer) = (MOVE_LEFT, MOVE_RIGHT, NO_MOVE)
RLBase.reward(env::Catcher, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::Catcher) = get_done(env)

function RLBase.reset!(env::Catcher{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)
    height = get_height(env)
    width = get_width(env)

    world[AGENT, get_agent_pos(env)] = false
    world[BALL, get_ball_pos(env)] = false

    new_ball_pos = CartesianIndex(1, rand(rng, 1:width))
    set_ball_pos!(env, new_ball_pos)
    world[BALL, new_ball_pos] = true

    agent_start_pos = CartesianIndex(height, rand(rng, 1:width))
    world[AGENT, agent_start_pos] = true
    set_agent_pos!(env, agent_start_pos)

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::Catcher{T})(action::Union{MoveLeft, MoveRight, NoMove}) where {T}
    world = get_world(env)
    height = get_height(env)

    old_agent_pos = get_agent_pos(env)
    world[AGENT, old_agent_pos] = false
    new_agent_pos = wrap_agent(env, move(action, old_agent_pos))
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    old_ball_pos = get_ball_pos(env)
    world[BALL, old_ball_pos] = false
    new_ball_pos = wrap_ball(env, move(DOWN, old_ball_pos))
    set_ball_pos!(env, new_ball_pos)
    world[BALL, new_ball_pos] = true

    if new_ball_pos[1] == height
        if new_ball_pos[2] == new_agent_pos[2]
            set_done!(env, false)
            set_reward!(env, env.ball_reward)
        else
            set_done!(env, true)
            set_reward!(env, get_terminal_reward(env))
        end
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function wrap_ball(env::Catcher, pos::CartesianIndex{2})
    height = get_height(env)
    width = get_width(env)
    rng = get_rng(env)
    i = pos[1]

    if i > height
        return CartesianIndex(1, rand(rng, 1:width))
    else
        return pos
    end
end

function wrap_agent(env::Catcher, pos::CartesianIndex{2})
    height = get_height(env)
    width = get_width(env)
    i = pos[1]
    j = pos[2]

    if j < 1
        return CartesianIndex(i, width)
    elseif j > width
        return CartesianIndex(i, 1)
    else
        return pos
    end
end


text:
function HMMbackward(v,phghm,pvgh;UseLogArray=true)
    #%HMMBACKWARD HMM Backward Pass
    #% logbeta=HMMbackward(v,phghm,pvgh)
    #%
    #% Inputs:
    #% v : visible (observation) sequence being a vector v=[2 1 3 3 1 ...]
    #% phghm : homogeneous transition distribution phghm(i,j)=p(h(t)=i|h(t-1)=j)
    #% pvgh : homogeneous emission disrtribution pvgh(i,j)=p(v(t)=i|h(t)=j)
    #%
    #% Outputs:
    #% logbeta: log beta messages: log p(v(t+1:T)|h(t))
    #           p(v(t:T)|h(t)) = sum_{h(t)} {p(v(t)| h(t)) * p(h(t)|h(t-1)) * p(v(t+1:T)|h(t))}
    #           beta(t-1) = p(v(t:T)|h(t)
    #           beta(t) = p(v(t+1:T)|h(t)
    #% See also HMMbackward.m, HMMviterbi.m, demoHMMinference.m

    T=length(v); H=size(phghm,1)

    if !UseLogArray # beta recursion (not recommended due to numerical underflow)
        beta=zeros(H,T)
        beta[:,T]=ones(H,1)
        for t=T:-1:2
            beta[:,t-1]=phghm'*(beta[:,t].*pvgh[v[t],:]')
        end
        return log(beta)
    end

    logpvgh=LogArray(log(pvgh))
    logphghm=LogArray(log(phghm))
    logbetaoutput=zeros(H,T)
    logbeta = LogArray(zeros(H,1))
    for t=2:T
        logbeta=logphghm.'*(logbeta.*logpvgh[v[t],:].')
        logbetaoutput[:,t]=(logbeta).content
    end

    return logbetaoutput

end


text:
function test_sqrt_all()
    print(rpad("Math.Series.[1/]sqrt(x)...", 50, ' '))

    sqrt = Brainstorm._Math._Series._Sqrt.sqrt
    rsqrt = Brainstorm._Math._Series._Sqrt.rsqrt

    setprecision(4_000_000) do
        @test all([sqrt(d, x)  - trunc(BigInt, Base.sqrt(  big(x))*(big(10)^d)) for x in (2, 3, 5, 17, 1_000_003, 1_000_033, 1_000_037), d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] .== 0)
        @test all([rsqrt(d, x) - trunc(BigInt, Base.sqrt(1/big(x))*(big(10)^d)) for x in (2, 3, 5, 17, 1_000_003, 1_000_033, 1_000_037), d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] .== 0)
    end

    println("PASS")
end


text:
using HDF5

abstract type Solution{dType, tType, N} end
abstract type DeterministicSolution{dType, tType, N} <: Solution{dType, tType, N} end
abstract type StochasticSolution{dType, tType, NQ, NW} <: Solution{dType, tType, NQ} end

hdf5(sol::Solution)   = error("hdf5() not implemented for ", typeof(sol))
time(sol::Solution)   = error("time() not implemented for ", typeof(sol))
ntime(sol::Solution)  = error("ntime() not implemented for ", typeof(sol))
nsave(sol::Solution)  = error("nsave() not implemented for ", typeof(sol))
offset(sol::Solution) = error("offset() not implemented for ", typeof(sol))

create_hdf5(sol::Solution, file) = error("create_hdf5() not implemented for ", typeof(sol))
CommonFunctions.write_to_hdf5(sol::Solution, h5::HDF5File, offset=0) = error("write_to_hdf5() not implemented for ", typeof(sol))

"Create solution for ODE and split ODE."
function Solution(equation::Union{ODE,SODE}, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SolutionODE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for partitioned ODE."
function Solution(equation::PODE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SolutionPODE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for variational ODE."
function Solution(equation::VODE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SolutionPODE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for implicit ODE."
function Solution(equation::IODE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SSolutionPDAE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for DAE."
function Solution(equation::DAE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SolutionDAE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for partitioned DAE."
function Solution(equation::PDAE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SSolutionPDAE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for implicit DAE."
function Solution(equation::IDAE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE)
    SSolutionPDAE(equation, Δt, ntime, nsave, nwrite)
end

"Create solution for SDE."
function Solution(equation::SDE, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE; K::Int=0, conv::String="strong")
    SolutionSDE(equation, Δt, ntime, nsave, K=K, conv=conv)
end

"Create solution for PSDE."
function Solution(equation::Union{PSDE,SPSDE}, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE; K::Int=0, conv::String="strong")
    SolutionPSDE(equation, Δt, ntime, nsave, K=K, conv=conv)
end

"Print error for solutions of equations not implemented, yet."
function Solution(equation::Equation, Δt, ntime::Int, nsave::Int=DEFAULT_NSAVE)
    error("No solution found for equation ", equation)
end

"createHDF5: Creates or opens HDF5 file."
function createHDF5(sol::Solution, file::AbstractString; overwrite=true)
    if overwrite
        flag = "w"
        get_config(:verbosity) > 1 ? @info("Creating HDF5 file ", file) : nothing
        isfile(file) ? @warn("Overwriting existing HDF5 file.") : nothing
    else
        flag = "r+"
        get_config(:verbosity) > 1 ? @info("Opening HDF5 file ", file) : nothing
    end

    # create or open HDF5 file
    h5 = h5open(file, flag)

    return h5
end

"save_attributes: Saves common attributes of Solution to HDF5 file."
function save_attributes(sol::Solution)
    save_attributes(sol, hdf5(sol))
end

"save_attributes: Saves attributes of Deterministic Solutions to HDF5 file."
function save_attributes(sol::DeterministicSolution, h5::HDF5File)
    # save attributes
    attrs(h5)["ntime"] = ntime(sol)
    attrs(h5)["nsave"] = nsave(sol)
end

"save_attributes: Saves attributes of Stochastic Solutions to HDF5 file."
function save_attributes(sol::StochasticSolution, h5::HDF5File)
    attrs(h5)["nsave"] = solution.nsave
    attrs(h5)["conv"] = solution.conv
    attrs(h5)["nd"] = solution.nd
    attrs(h5)["nm"] = solution.nm
    attrs(h5)["ns"] = solution.ns
    attrs(h5)["ni"] = solution.ni
    attrs(h5)["K"] = solution.K
end

"write_to_hdf5: Wrapper for saving Solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::Solution)
    write_to_hdf5(solution, hdf5(solution), offset(solution))
end

"Creates HDF5 file, writes solution to file, and closes file."
function CommonFunctions.write_to_hdf5(solution::Solution, file::AbstractString)
    h5 = create_hdf5(solution, file)
    write_to_hdf5(solution, h5)
    close(h5)
end


text:
x = 42
while x > 0
	println(x)
	x = x - 1
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
#**************************************************************************************
# test_Input_Console.jl
# ======================= part of the GeoEfficiency.jl package.
#
# test all the input to the package from the console.
#
#**************************************************************************************

using .MathConstants

@debug("GeoEfficiency.input")
@testset "GeoEfficiency.input" begin
	@test H.exec_console_unattended(G.input, []) == ""
    @test H.exec_console_unattended(G.input, "") == ""
    @test H.exec_console_unattended(G.input, [1]) == "1"
	@test H.exec_console_unattended(G.input, "1") == "1"
	@test H.exec_console_unattended(G.input,  "anbfyiQERFC") == "anbfyiQERFC"
	@test readavailable(stdin.buffer) |> String == ""   # test that no thing is left in the stdin
end #testset_input

@debug("GeoEfficiency.getfloat")
@testset "GeoEfficiency.getfloat" begin

	@debug("getfloat - different ways to input numbers")
	@test   0.0     ==  G.getfloat("\njust press return: ",value="")
	@test   0.0     ==  G.getfloat("\ninput '0.0', then press return: ",value="0.0")
	@test   0.0     ==  G.getfloat("\ninput '0', then press return: ",value="0")
	@test   0.0     ==  G.getfloat("\ninput '+0', then press return: ",value="+0")
	@test   0.0     ==  G.getfloat("\ninput '-0', then press return: ",value="-0")
	@test   0.0     ==  G.getfloat("\ninput '1*0im', then press return: ",value="1*0im")

	@test   1.0     ==  G.getfloat("\njust press return: ",value="1.0")
	@test   1.0     ==  G.getfloat("\ninput 1, then press return: ",value="+1.0")
	@test   -1.0    ==  G.getfloat("\ninput 1.0, then press return: ",value="-1.0")
	@test   1.0     ==  G.getfloat("\ninput '1*0im', then press return: ",value="1+1*0im")

	@test   2000.0  ==  G.getfloat("\ninput '2e3', then press return: ",value="2e3")
	@test   2000.0  ==  G.getfloat("\ninput '2e3', then press return: ",value="2E3")
	@test   -2000.0 ==  G.getfloat("\ninput '-2e3', then press return: ",value="-2e3")
	@test   -2000.0 ==  G.getfloat("\ninput '-2e3', then press return: ",value="-2E3")
	@test   0.034   ==  G.getfloat("\ninput '3.4e-2', then press return: ",value="3.4e-2")
	@test   -0.034  ==  G.getfloat("\ninput '-3.4e-2', then press return: ",value="-3.4e-2")

	@debug("getfloat - mathematical expressions")
	@test   0.5           		==  G.getfloat("\ninput 1/2, then press return: ",value="1/2")
	@test   0.75          		==  G.getfloat("\ninput 3//4, then press return: ",value="3//4")
	@test   pi/2 		        ≈   G.getfloat("\ninput 'pi/2', then press return: ",value="pi/2")
	@test   e                   ≈   G.getfloat("\ninput 'e', then press return: ",value="e")
	@test   e^3                 ≈   G.getfloat("\ninput 'e^3', then press return: ",value="e^3")
	@test   Base.sin(0.1) 		≈   G.getfloat("\ninput 'sin(0.1)', then press return: ",value="sin(0.1)")
	@test   Base.sin(0.1) + e^3 ≈   G.getfloat("\ninput 'sin(0.1)', then press return: ",value="sin(0.1)+e^3")

	@debug("getfloat - invalide console input")
	readavailable(stdin.buffer) == UInt8[] || @warn "buffer not empty, see the pervious process to `stdin`"  buffer = String(bffr)
	write(stdin.buffer,"\n"); 				@test 0.0 == G.getfloat("\njust press return: ")	# valide input but for completness
	write(stdin.buffer,"\n" * "3\n"); 		@test 3.0 == G.getfloat("\nthe first time just press return, then input `3` : ", 0.1, 4.0)
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput 5, then press return: ", 0.0, 4.0)
	write(stdin.buffer,"-1\n" * "3\n"); 	@test 3.0 == G.getfloat("\ninput -1, then press return: ", 0.0, 4.0)

	write(stdin.buffer,"1.2f\n" * "3\n");	 @test 3.0 == G.getfloat("\nthe first time input '1.2f', then input `3` : ")
	write(stdin.buffer,"one\n" * "3\n");	 @test 3.0 == G.getfloat("\nthe first time input 'one', then input `3` : ") # trying to input any string, only valid number should accepted.
	write(stdin.buffer,"3.4e-2.0\n" * "3\n");@test 3.0  == G.getfloat("\nthe first time input '3.4e-2.0', then input `3` : ")

	for i = 0:5
		write(stdin.buffer,"1.2+2im\n"^i * "3\n")
		@test G.getfloat("\nthe first time input '1.2+2im': ") == 3.0
	end # for

	@debug("getfloat - interval boundaries")
	# '1', '5' are to represent the lower, upper limits; while '3' represents an intermediate value.

	@test 1.0 == G.getfloat("\ninput '1', then press Return: ", 1.0, 5.0, value="1")					# by default lower limit is valide
	@test 1.0 == G.getfloat("\ninput '1', then press Return: ", 1.0, 5.0, value="1", lower=true)		# now explicit lower
	@test 1.0 == G.getfloat("\ninput '1', then press Return: ", 1.0, 5.0, value="1", lower=true, upper=false)	# now explicit lower - explicit upper
	@test 1.0 == G.getfloat("\ninput '1', then press Return: ", 1.0, 5.0, value="1", lower=true, upper=true)	# now explicit lower- modified upper

	write(stdin.buffer,"1\n" * "3\n"); 		@test 3.0 == G.getfloat("\ninput '1' & Return, then '3' & Return: ", 1.0, 5.0, lower=false)					# now lower limit `1` is not valide
	write(stdin.buffer,"1\n" * "3\n"); 		@test 3.0 == G.getfloat("\ninput '1' & Return, then '3' & Return: ", 1.0, 5.0, lower=false, upper=false)	# now lower limit `1` is not valide - explicit upper
	write(stdin.buffer,"1\n" * "3\n"); 		@test 3.0 == G.getfloat("\ninput '1' & Return, then '3' & Return: ", 1.0, 5.0, lower=false, upper=true)	# now lower limit `1` is not valide - modified upper

	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0)		# by default upper limit is not valide
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0, lower=true)		# by default upper limit is not valide - explicit lower
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0, lower=false)		# by default upper limit is not valide - modified lower
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0, 			  upper=false)		# explicit upper limit is not valide
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0, lower=true, upper=false,)		# explicit upper limit is not valide - explicit lower
	write(stdin.buffer,"5\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5' & Return, then '3' & Return:", 1.0, 5.0, lower=false, upper=false, )	# explicit upper limit is not valide - modified lower

	@test 5.0 == G.getfloat("\ninput 5, then press Return: ", 1.0, 5.0, value="5", upper=true)		# now upper limit `5` is valide
	@test 5.0 == G.getfloat("\ninput 5, then press Return: ", 1.0, 5.0, value="5", upper=true, lower=true)		# now input `5` is valide	 - explicit lower
	@test 5.0 == G.getfloat("\ninput 5, then press Return: ", 1.0, 5.0, value="5", upper=true, lower=false)		# now input `5` is valide	 - modified lower

	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3")		# by default intermediate value accepted
	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=true)
	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=true, upper=false)
	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=true, upper=true)

	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=false)
	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=false, upper=false)
	@test 3.0 == G.getfloat("\ninput '3', then press Return: ", 1.0, 5.0, value="3", lower=false, upper=true)

	write(stdin.buffer,"5\n" * "1\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5', then '1', then '3': ", 3.0, 3.0)								# when only one value (ex. 3.0) accepted
	write(stdin.buffer,"5\n" * "1\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5', then '1', then '3': ", 3.0, 3.0, 			 upper =true)	# when only one value (ex. 3.0) accepted
	write(stdin.buffer,"5\n" * "1\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5', then '1', then '3': ", 3.0, 3.0, lower =true, upper =true)	# when only one value (ex. 3.0) accepted
	write(stdin.buffer,"5\n" * "1\n" * "3\n");  	@test 3.0 == G.getfloat("\ninput '5', then '1', then '3': ", 3.0, 3.0, lower =false, upper =true)	# when only one value (ex. 3.0) accepted

	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 3.0, 3.0, lower =false, upper =false, value="3")
	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 5.0, 1.0, value="3")
	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 5.0, 1.0, lower =true, value="3")
	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 5.0, 1.0, lower =true, upper =false, value="3")
	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 5.0, 1.0, lower =false, upper =false, value="3")
	@test_throws 	ArgumentError	3.0 == G.getfloat("\ninput '3', then press Return: ", 5.0, 1.0, lower =false, upper =true, value="3")
end #testset_getfloat


text:
using Test

using Beamforming
using Beamforming.Synth
import Seis

@testset "Synth" begin
    @testset "synthetic_arrival" begin
        let t = synthetic_arrival(1, 2)
            @test length(t) == 20
            @test (-2.5, -2.5) <= extrema(t.sta.lon) <= (2.5, 2.5)
            @test t[1].delta == 0.01
            @test minimum(minimum, trace.(t)) == 0
            @test all(x -> Seis.nsamples(x) == 10001, t)
            @test all(x -> x === missing, t.evt.lon)
        end
        let t = synthetic_arrival(1, 2, delta=0.02, nsta=10, lon=rand(10),
                lat=rand(10), distance=80, noise=0.1)
            @test length(t) == 10
            @test (0, 0) <= extrema(t.sta.lon) <= (1, 1)
            @test t[1].delta == 0.02
            @test minimum(minimum, trace.(t)) < 0
            @test all(x -> Seis.nsamples(x) == 10001, t)
            @test all(x -> x !== missing, t.evt.lon)
            @test Beamforming.delta(t[1].evt.lon, t[1].evt.lat,
                t[1].sta.lon, t[1].sta.lon, true) ≈ 80 atol=2
            @test all(isapprox.(Seis.backazimuth.(t), atand(1, 2), atol=2))
        end
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
function quarter_date_to_number(date::Date)
    y = Dates.year(date)
    q = Dates.quarterofyear(date)
    return y + (q-1)*0.25
end

function quarter_number_to_date(datenum::Real)
    if datenum % 0.25 != 0
        throw(DomainError())
    end

    y = convert(Int, floor(datenum))
    q = convert(Int, (datenum % 1) / 0.25) + 1
    return quartertodate("$y-Q$q")
end

function plot_extension()
    be = typeof(Plots.backend())
    if be == Plots.GRBackend
        :pdf
    elseif be in [Plots.PlotlyBackend, Plots.PlotlyJSBackend]
        :html
    else
        :pdf
    end
end

function describe_series(m::AbstractDSGEModel, var::Symbol, class::Symbol;
                         detexify::Bool = false)
    res = if class in [:obs, :pseudo]
        dict = if class == :obs
            m.observable_mappings
        elseif class == :pseudo
            m.pseudo_observable_mappings
        end
        dict[var].name
    elseif class == :states
        string(var)
    elseif class in [:shocks, :stdshocks]
        replace(string(var), r"_sh$" => "")
    else
        error("Invalid class: " * string(class))
    end

    detexify ? DSGE.detexify(res) : res
end

function describe_series(m::AbstractDSGEVARModel, var::Symbol, class::Symbol;
                         detexify::Bool = false)
    return describe_series(get_dsge(m), var, class; detexify = detexify)
end

function series_ylabel(m::AbstractDSGEModel, var::Symbol, class::Symbol;
                       untrans::Bool = false,
                       fourquarter::Bool = false)
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    end

    if class in [:obs, :pseudo]
        dict = if class == :obs
            m.observable_mappings
        elseif class == :pseudo
            m.pseudo_observable_mappings
        end
        transform = dict[var].rev_transform

        if transform in [loggrowthtopct_annualized_percapita, loggrowthtopct_annualized, loggrowthtopct, loggrowthtopct_percapita]
            if untrans
                return "Q/Q Log Growth Rate"
            elseif fourquarter
                return "Percent 4Q Growth"
            else
                return "Percent Q/Q Annualized"
            end
        elseif transform in [logleveltopct_annualized_percapita, logleveltopct_annualized]
            if untrans
                return "Log Level"
            elseif fourquarter
                return "Percent 4Q Growth"
            else
                return "Percent Q/Q Annualized"
            end
        elseif transform in [quartertoannual]
            if untrans
                return "Percent Q/Q"
            else
                return "Percent Annualized"
            end
        elseif transform == identity
            ""
        else
            error("series_ylabel not implemented for transform: $transform")
        end
    elseif class == :stdshocks
        return "Standard Deviations"
    elseif class in [:states, :shocks]
        return ""
    else
        error("Invalid class: " * string(class))
    end
end

function save_plot(output_file::String = ""; verbose::Symbol = :low)
    if !isempty(output_file)
        output_dir = dirname(output_file)
        !isdir(output_dir) && mkpath(output_dir)
        @show output_file
        Plots.savefig(output_file)

        println(verbose, :low, "Saved $output_file")
    end
end

function save_plot(p::Plots.Plot, output_file::String = ""; verbose::Symbol = :low)
    if !isempty(output_file)
        output_dir = dirname(output_file)
        !isdir(output_dir) && mkpath(output_dir)
        @show output_file
        Plots.savefig(p, output_file)

        println(verbose, :low, "Saved $output_file")
    end
end

function rescale_annotations_fontsize!(p::Plots.Plot{Plots.GRBackend}, fontsize::Int64)
    for i in 1:length(p.series_list)
        if typeof(p.series_list[i].plotattributes[:series_annotations]) == Nothing
            continue
        else
            p.series_list[i].plotattributes[:series_annotations].font.pointsize = fontsize
        end
    end
end

# horiz_position can be :hcenter, :left, :right
# vert_position can be :vcenter, :top, :bottom
# reverse_position Bool to indicate if you want left and right to actually correspond to left and right..
function relocate_annotations!(p::Plots.Plot{Plots.GRBackend},
                               horiz_position::Symbol, vert_position::Symbol,
                               reverse_position::Bool = true)
    if reverse_position
        horiz_position =
        if horiz_position == :right
            :left
        elseif horiz_position == :left
            :right
        else
            :hcenter
        end
        vert_position =
        if vert_position == :top
            :bottom
        elseif vert_position == :bottom
            :top
        else
            :vcenter
        end
    end

    for i in 1:length(p.series_list)
        if typeof(p.series_list[i].plotattributes[:series_annotations]) == Nothing
            continue
        else
            p.series_list[i].plotattributes[:series_annotations].font.halign = horiz_position
            p.series_list[i].plotattributes[:series_annotations].font.valign = vert_position
        end
    end
end

function date_to_float(datevec::Vector{Date}, reps::Int64 = 1)
    date_floats = zeros(reps, length(datevec))
    tofloats = Dict{String,Float64}("03" => .25, "06" => .5, "09" => .75, "12" => 1.)
    for j = 1:length(datevec)
        tmp = string(datevec[j])
        s = tmp[6:7]
        num = parse(Float64,tmp[1:4]) + tofloats[s]
        date_floats[:,j] .= num
    end
    return date_floats
end


text:
println("	MODERN")
@test GMT.helper_sub_F("1/2") == "1/2"
println("    SUBPLOT1")
subplot(grid=(1,1), limits="0/5/0/5", frame="west", F=([1 2]), name=:png)
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", frame="west", F=([1 2]), Vd=dbg2), "-F1/2")
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", frame="west", F=("1i",2), Vd=dbg2), "-F1i/2")
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", frame="west", F=(size=(1,2), frac=((2,3),(3,4,5))), figname="lixo.ps", Vd=dbg2), "-Ff1/2+f2,3/3,4,5")
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", frame="west", F=(width=1,height=5,fwidth=(0.5,1),fheight=(10,), fill=:red, outline=(3,:red)), Vd=dbg2), "-F1/5+f0.5,1/10+gred+p3,red")
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", frame="west", dims=(panels=((2,4),(2.5,5,1.25)),fill=:red), Vd=dbg2), "-Fs2,4/2.5,5,1.25+gred")
@test endswith(subplot(grid=(1,1), limits="0/5/0/5", F=(panels=((5,8),8),), Vd=dbg2), "-Fs5,8/8")
@test GMT.helper_sub_F((width=1,)) == "1/1"
@test GMT.helper_sub_F((width=1,height=5,fwidth=(0.5,1),fheight=(10,))) == "1/5+f0.5,1/10"
@test_throws ErrorException("SUBPLOT: when using 'fwidth' must also set 'fheight'") GMT.helper_sub_F((width=1,height=5,fwidth=(0.5,1)))
@test_throws ErrorException("'frac' option must be a tuple(tuple, tuple)") subplot(grid=(1,1),  F=(size=(1,2), frac=((2,3))), Vd=dbg2)
@test_throws ErrorException("SUBPLOT: garbage in DIMS option") GMT.helper_sub_F([1 2 3])
@test_throws ErrorException("SUBPLOT: 'grid' keyword is mandatory") subplot(F=("1i"), Vd=dbg2)
@test_throws ErrorException("Cannot call subplot(set, ...) before setting dimensions") subplot(:set, F=("1i"), Vd=dbg2)
println("    SUBPLOT2")
subplot(name="lixo", fmt=:ps, grid="1x1", limits="0/5/0/5", frame="west", F="s7/7", title="VERY VERY");subplot(:set, panel=(1,1));plot([0 0; 1 1]);subplot(:end)
gmtbegin("lixo.ps");  gmtend()
gmtbegin("lixo", fmt=:ps);  gmtend()
gmtbegin("lixo");  gmtend()
gmtbegin();  gmtend()

println("    BEGINEND")
gmtbegin("lixo.ps")
	basemap(region=(0,40,20,60), proj=:merc, figsize=16, frame=(annot=:afg, fill=:lightgreen))
	inset(D="jTR+w2.5i+o0.2i", F="+gpink+p0.5p", margins=0.6)
	basemap(region=:global360, J="A20/20/2i", frame=:afg)
	text(text_record([1 1],["INSET"]), font=18, region_justify=:TR, D="j-0.15i", noclip=true)
	inset(:end)
	text(text_record([0 0; 1 1.1],[" ";" "]), text="MAP", font=18, region_justify=:BL, D="j0.2i")
gmtend()

gmtbegin(); gmtfig("lixo.ps -Vq");	gmtend()

println("    MOVIE")
movie("main_sc.jl", pre="pre_sc.jl", C="7.2ix4.8ix100", N=:anim04, T="flight_path.txt", L="f+o0.1i", F=:mp4, A="+l+s10", Sf="", Vd=dbg2)
@test GMT.helper_fgbg("", "bla", "bla", " -Sf") == " -Sfbla"
@test_throws ErrorException("bla script has nothing useful") GMT.helper_fgbg("", rand, "bla", " -Sf")
GMT.resetGMT()		# This one is needed to reset the broken state left by calling helper_fgbg() directly
if (Sys.iswindows())
	rm("pre_script.bat");		rm("main_script.bat")
else
	rm("pre_script.sh");		rm("main_script.sh")
end
println("    EVENTS")
events("", R=:g, J="G200/5/6i", B=:af, S="E-", C=true, T="2018-05-01T", E="s+r2+d6", M=((size=5,coda=0.5),(intensity=true,)), Vd=dbg2);


text:
"""
Converts String to Integer
"""
function process_string!(message::String)
    acc::BigInt = BigInt(0)
    padding = 0
    len = length(message)

    if (len % 4 != 0)
        padding = (4 - (len % 4))
        message = repeat("\x00", padding) * message
    end

    for i in 1:4:len
        message_part = Vector{UInt8}(message[i:i + 3])
        message_part = reinterpret(UInt32, message_part)[1]
        acc = (acc << 32) + message_part
    end
    acc
end

"""
Convert Integer to String
"""
function recover_string!(number::BigInt)
    s = ""
    i = 1

    while number > 0
        s = String(to_bytes([UInt32(number & 0xffffffff)])) * s
        number = number >> 32
    end
    if (s == "")
        return (s)
    end
    while UInt8(s[i]) == 0
        i += 1
    end
    s[i:length(s)]
end