text:
module kara_gui
using Test
using JuliaKara

path = joinpath(@__DIR__,"..","test","example.world")

cleanup = !JuliaKara.Blink.AtomShell.isinstalled()
cleanup && JuliaKara.Blink.AtomShell.install()

@World wtest (10,10)

@testset "JuliaKara Example 01" begin

    kara = place_kara(wtest,4,4)

    function nextTest(world,kara)
        turnRight(world,kara)
        move(world,kara)
        turnLeft(world,kara)
    end

    @wtest place_tree(4,5)
    @wtest place_tree(5,6)
    @wtest place_mushroom(5,5)
    @wtest place_mushroom(6,5)
    @wtest place_mushroom(6,6)
    @wtest place_mushroom(7,5)
    @wtest place_leaf(8,5)

    @test_throws JuliaKara.JuliaKara_noGUI.ActorInvalidMovementError move(wtest,kara)
    nextTest(wtest,kara)
    @test_throws JuliaKara.JuliaKara_noGUI.ActorInvalidMovementError move(wtest,kara)
    nextTest(wtest,kara)
    @test_throws JuliaKara.JuliaKara_noGUI.ActorInvalidMultipleMovementError move(wtest,kara)
    nextTest(wtest,kara)
    move(wtest,kara)
    @test (mushroomFront(wtest,kara)) == true
    turnLeft(wtest,kara)
    turnLeft(wtest,kara)
    move(wtest,kara)
    turnLeft(wtest,kara)
    turnLeft(wtest,kara)
    @test (mushroomFront(wtest,kara)) == false
    nextTest(wtest,kara)
    move(wtest,kara)
    @test (onLeaf(wtest,kara)) == true
    removeLeaf(wtest,kara)
    @test (onLeaf(wtest,kara)) == false
    @test_throws JuliaKara.JuliaKara_noGUI.ActorGrabNotFoundError removeLeaf(wtest,kara)
    putLeaf(wtest,kara)
    @test (onLeaf(wtest,kara)) == true
    @test (treeLeft(wtest,kara)) == false
    @test (treeRight(wtest,kara)) == false
    @test (treeFront(wtest,kara)) == false
    place_tree(wtest,7,5)
    place_tree(wtest,9,5)
    place_tree(wtest,8,6)
    @test (treeLeft(wtest,kara)) == true
    @test (treeRight(wtest,kara)) == true
    @test (treeFront(wtest,kara)) == true
    # Test get kara
    @test get_kara(wtest) == kara
    place_kara(wtest,10,10)
    @test length(get_kara(wtest)) == 2

end

@World path
@World wcompare path
lara = @wcompare get_kara()

@testset "Loading Types" begin
    # Test if loaded worlds are equal
    for i in 1:length(world.world.actors)
        @test world.world.actors[i].actor_definition == wcompare.world.actors[i].actor_definition
        @test world.world.actors[i].location == wcompare.world.actors[i].location
    end
    # Move kara in global world using wrapper
    @test move(kara) == nothing
    @test JuliaKara.JuliaKara_noGUI.get_actors_at_location(
        world.world,
        JuliaKara.JuliaKara_noGUI.Location(4,8 )
    )[1].actor_definition == JuliaKara.JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
    # Test all other functions from global scope
    @test turnLeft(kara) == nothing
    @test turnRight(kara) == nothing
    @test putLeaf(kara) == nothing
    @test removeLeaf(kara) == nothing
    @test onLeaf(kara) == false
    @test treeFront(kara) == false
    @test treeRight(kara) == false
    @test treeLeft(kara) == false
    @test mushroomFront(kara) == false
    # Test reset
    reset!(world)
    @test JuliaKara.JuliaKara_noGUI.get_actors_at_location(
        world.world,
        JuliaKara.JuliaKara_noGUI.Location(3,8)
    )[1].actor_definition == JuliaKara.JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
    move(kara)
    store!(world)
    reset!(world)
    @test JuliaKara.JuliaKara_noGUI.get_actors_at_location(
        world.world,
        JuliaKara.JuliaKara_noGUI.Location(4,8)
    )[1].actor_definition == JuliaKara.JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
    # Test other loaded world

    @wcompare move(lara)
    @test JuliaKara.JuliaKara_noGUI.get_actors_at_location(
        wcompare.world,
        JuliaKara.JuliaKara_noGUI.Location(4,8)
    )[1].actor_definition == JuliaKara.JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
end

JuliaKara.close(wtest.window)
JuliaKara.close(wcompare.window)
JuliaKara.close(world.window)
end


text:
module Fixtures

    import Log4jl
    import Log4jl: level, marker, append!

    export IncompleteMessage, InvalidConfiguration,
           TestEvent, level, marker,
           InvalidAppender, TestAppender, append!

    const TESTNAME   = "TEST"
    const TESTLEVEL  = Log4jl.Level.INFO
    const TESTLEVELN  = Log4jl.LEVEL(TESTLEVEL)
    const TESTMARKER = :TESTMARKER

    type IncompleteMessage <: Log4jl.Message
    end
    TESTMSG = Log4jl.Messages.SimpleMessage(Fixtures.TESTNAME)

    type InvalidConfiguration <: Log4jl.Configuration
    end

    type TestEvent <: Log4jl.Event
    end
    level(evnt::TestEvent) = TESTLEVEL
    marker(evnt::TestEvent) = Log4jl.MARKER(TESTMARKER)

    type TestFilter <: Log4jl.Filter
    end
    TESTFILTER = TestFilter()

    type InvalidAppender <: Log4jl.Appender
    end

    type TestAppender <: Log4jl.Appender
        name::AbstractString
        layout::Log4jl.LAYOUT
        filter::Log4jl.FILTER
        state::Log4jl.LifeCycle.State
    end
    TestAppender(cfg::Dict{Symbol, Any}) = TestAppender(TESTNAME, Log4jl.LAYOUT(), Log4jl.FILTER(), Log4jl.LifeCycle.INITIALIZED)
    TestAppender() = TestAppender(Dict{Symbol, Any}())
    append!(apnd::TestAppender, evnt::Log4jl.Event) = evnt

end

ext: jl


text:
export linear, quadratic, quaddiff, cubic,
       neville, vandermonde, vandermonde!, cubichermite,
       LinearInterpolator, CubicInterpolator,
       BilinearInterpolator, BicubicInterpolator

"""
    linear(x, xₚ, yₚ)

Perform simple linear interpolation of the points defined by coordinates `xₚ` and values `yₚ`, at the coordinate `x`. `xₚ` and `yₚ` must both contain two points.
"""
function linear(x, xₚ, yₚ)
    @inbounds x₁, x₂ = xₚ[1], xₚ[2]
    @inbounds y₁, y₂ = yₚ[1], yₚ[2]
    return y₁ + (x - x₁)*(y₂ - y₁)/(x₂ - x₁)
end

"""
    quadratic(x, xₚ, yₚ)

Perform quadratic polynomial interpolation of the points defined by coordinates `xₚ` and values `yₚ`, at the coordinate `x`, using Neville's algorithm. `xₚ` and `yₚ` must both contain three points.
"""
function quadratic(x, xₚ, yₚ)
    @assert length(xₚ) == length(yₚ) == 3 "need 3 points for quadratic interpolation"
    #short names
    @inbounds x₁, x₂, x₃ = xₚ[1], xₚ[2], xₚ[3]
    @inbounds y₁, y₂, y₃ = yₚ[1], yₚ[2], yₚ[3]
    #first stage
    p₁₂ = ((x - x₂)*y₁ + (x₁ - x)*y₂)/(x₁ - x₂)
    p₂₃ = ((x - x₃)*y₂ + (x₂ - x)*y₃)/(x₂ - x₃)
    #final stage
    ((x - x₃)*p₁₂ + (x₁ - x)*p₂₃)/(x₁ - x₃)
end

"""
    cubic(x, xₚ, yₚ)

Perform cubic polynomial interpolation of the points defined by coordinates `xₚ` and values `yₚ`, at the coordinate `x`, using Neville's algorithm. `xₚ` and `yₚ` must both contain four points.
"""
function cubic(x, xₚ, yₚ)
    @assert length(xₚ) == length(yₚ) == 4 "need 4 points for cubic interpolation"
    #short names
    @inbounds x₁, x₂, x₃, x₄ = xₚ[1], xₚ[2], xₚ[3], xₚ[4]
    @inbounds y₁, y₂, y₃, y₄ = yₚ[1], yₚ[2], yₚ[3], yₚ[4]
    #first stage
    p₁₂ = ((x - x₂)*y₁ + (x₁ - x)*y₂)/(x₁ - x₂)
    p₂₃ = ((x - x₃)*y₂ + (x₂ - x)*y₃)/(x₂ - x₃)
    p₃₄ = ((x - x₄)*y₃ + (x₃ - x)*y₄)/(x₃ - x₄)
    #second stage
    p₁₂₃ = ((x - x₃)*p₁₂ + (x₁ - x)*p₂₃)/(x₁ - x₃)
    p₂₃₄ = ((x - x₄)*p₂₃ + (x₂ - x)*p₃₄)/(x₂ - x₄)
    #final stage
    ((x - x₄)*p₁₂₃ + (x₁ - x)*p₂₃₄)/(x₁ - x₄)
end

"""
    neville(x, xₚ, yₚ)

Perform polynomial interpolation of the points defined by coordinates `xₚ` and values `yₚ`, at the coordinate `x`, using Neville's algorithm with as many points as are provided. `xₚ` and `yₚ` must have the same length. With only 3 or 4 points the [`quadratic`](@ref) and [`cubic`](@ref) functions will be considerably faster.
"""
function neville(x, xₚ, yₚ)
    @assert length(xₚ) == length(yₚ) "can't Neville with vectors of different lengths"
    n = length(xₚ)
    P = zeros(n, n)
    P[:,1] = yₚ
    for i ∈ 2:n, j ∈ i:n
        @inbounds P[j,i] = ((x - xₚ[j-i+1])*P[j,i-1] - (x - xₚ[j])*P[j-1,i-1])/(xₚ[j] - xₚ[j-i+1])
    end
    return P[n,n]
end

"""
    vandermonde(x, y)

Generate the coefficients of an arbitrary order polynomial passing through the ponts defined by coordinates `x` and value `y`. For n points, n coefficients ``[c_0, c_1, ..., c_{n-1}]`` are returned forming the polynomial ``c_0 + c_1x + ... + c_{n-1}x^{n-1}``

!!! warning

    Solving for the the coefficients of a high-order polynomial is a notoriously ill-conditioned problem. It is not recommended for orders greater than 5 or 6, although it depends on the application. If you must interpolate with a high-order polynomial, it's better to use the [`neville`](@ref) function instead of computing coefficients.
"""
function vandermonde(x, y)
    @assert length(x) == length(y) "length of x must equal length of y"
    c = zeros(length(x))
    vandermonde!(c, x, y)
    return c
end

function vandermonde!(c, x, y)
    @assert length(c) == length(x) == length(y) "vandermonde! requires length(c) == length(x) == length(y)"
    n = length(x)
    s = similar(c)
    c[:] .= zero(eltype(c))
    s[n] = -x[1]
    @inbounds for i = 2:n
        for j = n-i+1:n-1
            s[j] -= x[i]*s[j+1]
        end
        s[n] -= x[i]
    end
    @inbounds for j = 1:n
        ϕ = n
        for k = n:-1:2
            ϕ = (k-1)*s[k] + x[j]*ϕ
        end
        f = y[j]/ϕ
        b = 1.0
        for k = n:-1:1
            c[k] += b*f
            b = s[k] + x[j]*b
        end
    end
end

"""
    cubichermite(x, x₁, x₂, y₁, y₂, y₁′, y₂′)

Interpolate a cubic polynomial between two points, given its values and first derivatives at the points.
"""
function cubichermite(x, x₁, x₂, y₁, y₂, y₁′, y₂′)
    Δx = (x₂ - x₁)
    ξ = (x - x₁)/Δx
    u = Δx*y₁′
    v = Δx*y₂′
    ξ^3*(2*y₁ + u - 2*y₂ + v) + ξ^2*(-3*y₁ - 2*u + 3*y₂ - v) + ξ*u + y₁
end

#-------------------------------------------------------------------------------
# simple piecewise linear interpolator

struct LinearInterpolator{T,B} <: OneDimensionalInterpolator
    r::InterpolatorRange{T}
    boundaries::B
    i::RefValue{Int64} #previous cell index
end

"""
    LinearInterpolator(x, y, boundaries=StrictBoundaries())

Construct a `LinearInterpolator` for the points defined by coordinates `x` and values `y`
"""
function LinearInterpolator(x, y, boundaries::AbstractBoundaries=StrictBoundaries())
    LinearInterpolator(InterpolatorRange(x, y), boundaries, Ref(1))
end

"""
    LinearInterpolator(f, xa, xb, n, boundaries=StrictBoundaries())

Construct a `LinearInterpolator` for the function `f` using `n` evenly spaced function evaluations in the range [`xa`,`xb`]
"""
function LinearInterpolator(f, xa, xb, n::Int, boundaries::AbstractBoundaries=StrictBoundaries())
    linstruct(LinearInterpolator, f, xa, xb, n, boundaries)
end

function (ϕ::LinearInterpolator)(x)
    #enforce boundaries
    ϕ.boundaries(x, ϕ.r.xa, ϕ.r.xb)
    #find the interpolation cell
    i = findcell(x, ϕ)
    #short names
    @inbounds x₁, x₂ = ϕ.r.x[i], ϕ.r.x[i+1]
    @inbounds y₁, y₂ = ϕ.r.y[i], ϕ.r.y[i+1]
    #interpolate
    (x - x₁)*(y₂ - y₁)/(x₂ - x₁) + y₁
end

Base.getindex(ϕ::LinearInterpolator, i) = ϕ.r.y[i]
Base.firstindex(::LinearInterpolator) = 1
Base.lastindex(ϕ::LinearInterpolator) = ϕ.r.n
function Base.setindex!(ϕ::LinearInterpolator, v, i)
    ϕ.r.y[i] = v
end
function Base.copy(ϕ::LinearInterpolator)
    LinearInterpolator(ϕ.r, ϕ.boundaries, Ref(1))
end

#-------------------------------------------------------------------------------
# piecewise cubic interpolator without continuous derivatives (not splines)

struct CubicInterpolator{T,B} <: OneDimensionalInterpolator
    r::InterpolatorRange{T}
    boundaries::B
    i::RefValue{Int64} #previous cell index
end

"""
    CubicInterpolator(x, y, boundaries=StrictBoundaries())

Construct a `CubicInterpolator` for the points defined by coordinates `x` and values `y`
"""
function CubicInterpolator(x, y, boundaries::AbstractBoundaries=StrictBoundaries())
    @assert length(x) > 3 "can't do cubic interpolation with < 4 points"
    CubicInterpolator(InterpolatorRange(x, y), boundaries, Ref(1))
end

"""
    CubicInterpolator(f, xa, xb, n, boundaries=StrictBoundaries())

Construct a `CubicInterpolator` for the function `f` using `n` evenly spaced function evaluations in the range [`xa`,`xb`]
"""
function CubicInterpolator(f, xa, xb, n::Int, boundaries::AbstractBoundaries=StrictBoundaries())
    linstruct(CubicInterpolator, f, xa, xb, n, boundaries)
end

function (ϕ::CubicInterpolator)(x)
    #enforce boundaries if desired
    ϕ.boundaries(x, ϕ.r.xa, ϕ.r.xb)
    #find the interpolation point
    i = findcell(x, ϕ)
    #determine which points to neville
    if i == 1
        I = 1:4
    elseif i == ϕ.r.n - 1
        I = ϕ.r.n-3:ϕ.r.n
    else
        I = i-1:i+2
    end
    #interpolate
    @inbounds cubic(x, view(ϕ.r.x,I), view(ϕ.r.y,I))
end

Base.getindex(ϕ::CubicInterpolator, i) = ϕ.r.y[i]
Base.firstindex(::CubicInterpolator) = 1
Base.lastindex(ϕ::CubicInterpolator) = ϕ.r.n
function Base.setindex!(ϕ::CubicInterpolator, v, i)
    ϕ.r.y[i] = v
end
function Base.copy(ϕ::CubicInterpolator)
    CubicInterpolator(ϕ.r, ϕ.boundaries, Ref(1))
end

#-------------------------------------------------------------------------------
# bilinear interpolator

struct BilinearInterpolator{T,B} <: TwoDimensionalInterpolator
    G::InterpolatorGrid{T}
    boundaries::B
    i::RefValue{Int64} #previous cell index
    j::RefValue{Int64} #previous cell index
end

"""
    BilinearInterpolator(x, y, Z, boundaries=StrictBoundaries())

Construct a `BilinearInterpolator` for the grid of points points defined by coordinates `x`,`y` and values `Z`.
"""
function BilinearInterpolator(x, y, Z, boundaries::AbstractBoundaries=StrictBoundaries())
    BilinearInterpolator(InterpolatorGrid(x, y, Z), boundaries, Ref(1), Ref(1))
end

"""
    BilinearInterpolator(f, xa, xb, nx, ya, yb, ny, boundaries=StrictBoundaries())

Construct a `BilinearInterpolator` for the function `f` using a grid of `nx` points evenly spaced on the first axis in [`xa`,`xb`] and `ny` points evenly spaced on the second axis in [`ya`,`yb`].
"""
function BilinearInterpolator(f,
                              xa, xb, nx::Int,
                              ya, yb, ny::Int,
                              boundaries::AbstractBoundaries=StrictBoundaries())
    linstruct(BilinearInterpolator, f, xa, xb, nx, ya, yb, ny, boundaries)
end

function (Φ::BilinearInterpolator)(x, y)
    #enforce boundaries if desired
    Φ.boundaries(x, Φ.G.xa, Φ.G.xb, y, Φ.G.ya, Φ.G.yb)
    #find the proper grid box to interpolate inside
    i, j = findcell(x, y, Φ)
    #clear names
    @inbounds x₁, x₂ = Φ.G.x[i], Φ.G.x[i+1]
    @inbounds y₁, y₂ = Φ.G.y[j], Φ.G.y[j+1]
    @inbounds Z₁₁, Z₂₁ = Φ.G.Z[i,j],     Φ.G.Z[i+1,j]
    @inbounds Z₂₂, Z₁₂ = Φ.G.Z[i+1,j+1], Φ.G.Z[i,j+1]
    #reused differences
    xx₁ = x - x₁
    x₂x₁ = x₂ - x₁
    #interpolate along axis 1 first
    t = xx₁*(Z₂₁ - Z₁₁)/x₂x₁ + Z₁₁
    u = xx₁*(Z₂₂ - Z₁₂)/x₂x₁ + Z₁₂
    #then finish by interpolating between the interpolated values along axis 2
    (y - y₁)*(u - t)/(y₂ - y₁) + t
end

Base.getindex(Φ::BilinearInterpolator, i, j) = Φ.G.Z[i,j]
function Base.setindex!(Φ::BilinearInterpolator, v, i, j)
    Φ.G.Z[i,j] = v
end
function Base.copy(Φ::BilinearInterpolator)
    BilinearInterpolator(Φ.G, Φ.boundaries, Ref(1), Ref(1))
end

#-------------------------------------------------------------------------------
# bicubic interpolator

struct BicubicInterpolator{T,B} <: TwoDimensionalInterpolator
    G::InterpolatorGrid{T}
    boundaries::B
    i::RefValue{Int64} #previous cell index
    j::RefValue{Int64} #previous cell index
end

"""
    BicubicInterpolator(x, y, Z, boundaries=StrictBoundaries())

Construct a `BicubicInterpolator` for the grid of points points defined by coordinates `x`,`y` and values `Z`.
"""
function BicubicInterpolator(x, y, Z, boundaries::AbstractBoundaries=StrictBoundaries())
    #insist on at least 4 points in each dimension
    @assert (length(x) > 3) & (length(y) > 3) "bicubic interpolation requires at least 4 points in each dimension"
    BicubicInterpolator(InterpolatorGrid(x, y, Z), boundaries, Ref(1), Ref(1))
end

"""
    BicubicInterpolator(f, xa, xb, nx, ya, yb, ny, boundaries=StrictBoundaries())

Construct a `BicubicInterpolator` for the function `f` using a grid of `nx` points evenly spaced on the first axis in [`xa`,`xb`] and `ny` points evenly spaced on the second axis in [`ya`,`yb`].
"""
function BicubicInterpolator(f,
                             xa, xb, nx::Int,
                             ya, yb, ny::Int,
                             boundaries::AbstractBoundaries=StrictBoundaries())
    linstruct(BicubicInterpolator, f, xa, xb, nx, ya, yb, ny, boundaries)
end

function (Φ::BicubicInterpolator)(x, y)
    #enforce boundaries if desired
    Φ.boundaries(x, Φ.G.xa, Φ.G.xb, y, Φ.G.ya, Φ.G.yb)
    #find the proper grid box to interpolate inside
    i, j = findcell(x, y, Φ)
    #get indices of points along axis 1 to use for interpolation
    if i == 1
        I = 1:4
    elseif i == Φ.G.nx - 1
        I = Φ.G.nx-3:Φ.G.nx
    else
        I = i-1:i+2
    end
    #get indices along axis 2 where initial 4 interpolations occur
    if j == 1
        J = 1:4
    elseif j == Φ.G.ny - 1
        J = Φ.G.ny-3:Φ.G.ny
    else
        J = j-1:j+2
    end
    #view the proper chunks of the arrays
    xᵣ = view(Φ.G.x, I)
    yᵣ = view(Φ.G.y, J)
    #perform initial 4 interpolations
    zₓ = (
        cubic(x, xᵣ, @view Φ.G.Z[I,J[1]]),
        cubic(x, xᵣ, @view Φ.G.Z[I,J[2]]),
        cubic(x, xᵣ, @view Φ.G.Z[I,J[3]]),
        cubic(x, xᵣ, @view Φ.G.Z[I,J[4]])
    )
    #final interpolation
    cubic(y, yᵣ, zₓ)
end

Base.getindex(Φ::BicubicInterpolator, i, j) = Φ.G.Z[i,j]
function Base.setindex!(Φ::BicubicInterpolator, v, i, j)
    Φ.G.Z[i,j] = v
end
function Base.copy(Φ::BicubicInterpolator)
    BicubicInterpolator(Φ.G, Φ.boundaries, Ref(1), Ref(1))
end


text:
using Constants

function predictor!(hydro, nq, jj)

    nq::Int = nq
	lm_ = dtm_ / dr

	#take care of ions
	for i = 1:neqi
		hydro.U1D_p[2:nq-1,i] = hydro.U1D[2:nq-1,i] -
                lm_ * (hydro.F1D[3:nq,i] - hydro.F1D[2:nq-1,i]) +
				dtm_ * (phi * hydro.G1D[3:nq,i] + (1-phi) * hydro.G1D[2:nq-1,i])
		if (geom=="spherical") #add correction for spherical geometry
			hydro.U1D_p[2:nq-1,i] = hydro.U1D_p[2:nq-1,i] - 2. * dtm_ / r[2:nq-1] .*
						(phi * (hydro.F1D[3:nq,i] - hydro.C1D[3:nq,i]) +
							(1-phi) * (hydro.F1D[2:nq-1,i] - hydro.C1D[2:nq-1,i]))
		end
	end

	#take care of electrons
	if mod(jj,2) == 1
		hydro.U1D_p[2:nq-1,neqi+1] = hydro.U1D[2:nq-1,neqi+1] -
                lm_ * (hydro.F1D[3:nq,neqi+1] - hydro.F1D[2:nq-1,neqi+1]) +
				dtm_ * (phi * hydro.G1D[3:nq,neqi+1] + (1-phi) * hydro.G1D[2:nq-1,neqi+1])
		if geom=="spherical" #add correction for spherical geometry
			hydro.U1D_p[2:nq-1,neqi+1] = hydro.U1D_p[2:nq-1,neqi+1] - 2. * dtm_ / r[2:nq-1] .*
						(phi * hydro.F1D[3:nq,neqi+1] + (1-phi) * hydro.F1D[2:nq-1,neqi+1])
		end
	else
		hydro.U1D_p[2:nq-1,neqi+1] = hydro.U1D[2:nq-1,neqi+1] -
                lm_ * (hydro.F1D[2:nq-1,neqi+1] - hydro.F1D[1:nq-2,neqi+1]) +
				dtm_ * ((1-phi) * hydro.G1D[2:nq-1,neqi+1] + phi * hydro.G1D[1:nq-2,neqi+1])
		if (geom=="spherical") #add correction for spherical geometry
			hydro.U1D_p[2:nq-1,neqi+1] = hydro.U1D_p[2:nq-1,neqi+1] - 2. * dtm_ / r[2:nq-1] .*
						((1-phi) * hydro.F1D[2:nq-1,neqi+1] + phi * hydro.F1D[1:nq-2,neqi+1])
		end
	end
    return hydro
end

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function corrector!(hydro, nq, jj)

    nq::Int = nq
	lm_ = dtm_ / dr

	for i = 1:neqi
		hydro.U1D_c[2:nq-1,i] = hydro.U1D_p[2:nq-1,i] -
                lm_ * (hydro.F1D[2:nq-1,i] - hydro.F1D[1:nq-2,i]) +
				dtm_ * ((1-phi) * hydro.G1D[2:nq-1,i] + phi * hydro.G1D[1:nq-2,i])
		if geom=="spherical"
			hydro.U1D_c[2:nq-1,i] = hydro.U1D_c[2:nq-1,i] - 2. * dtm_ ./ r[2:nq-1] .*
						((1-phi) * (hydro.F1D[2:nq-1,i] - hydro.C1D[2:nq-1,i]) +
						phi * (hydro.F1D[1:nq-2,i] - hydro.C1D[1:nq-2,i]))
		end
	end

	#take care of electrons
	if mod(jj,2) == 1
		hydro.U1D_c[2:nq-1,neqi+1] = hydro.U1D_p[2:nq-1,neqi+1] -
                lm_ * (hydro.F1D[2:nq-1,neqi+1] - hydro.F1D[1:nq-2,neqi+1])  +
				dtm_ * ((1-phi) * hydro.G1D[2:nq-1,neqi+1] + phi * hydro.G1D[1:nq-2,neqi+1])
		if geom=="spherical" #add correction for spherical geometry
			hydro.U1D_c[2:nq-1,neqi+1] = hydro.U1D_c[2:nq-1,neqi+1] - 2. * dtm_ ./ r[2:nq-1]  .*
						((1-phi) * hydro.F1D[2:nq-1,neqi+1] + phi * hydro.F1D[1:nq-2,neqi+1])
		end
	else
		hydro.U1D_c[2:nq-1,neqi+1] = hydro.U1D_p[2:nq-1,neqi+1] -
                lm_ *  (hydro.F1D[3:nq,neqi+1] - hydro.F1D[2:nq-1,neqi+1]) +
				dtm_ * (phi * hydro.G1D[3:nq,neqi+1] + (1-phi) * hydro.G1D[2:nq-1,neqi+1])
		if geom=="spherical" #add correction for spherical geometry
			hydro.U1D_c[2:nq-1,neqi+1] = hydro.U1D_c[2:nq-1,neqi+1] - 2. * dtm_ ./ r[2:nq-1] .*
						(phi * hydro.F1D[3:nq,neqi+1] + (1-phi) * hydro.F1D[2:nq-1,neqi+1])
		end
	end

    return hydro
end

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function update_variables!(UU, hydro)
	# need to initialize electron density and speed to zero
	hydro.rho[:,nspec+1] = 0.
	hydro.ne = zeros(Float64, nz)
	hydro.u[:,nspec+1] = 0.

	for i = 1:nspec
		hydro.rho[:,i] = UU[:,3*(i-1)+1]
		hydro.u[:,i] = UU[:,3*(i-1)+2] ./ hydro.rho[:,i]
		hydro.T[:,i] = (g-1) * mi[i] * (UU[:,3*(i-1)+3] ./ hydro.rho[:,i] - 0.5*hydro.u[:,i].^2)
		hydro.p[:,i] = hydro.rho[:,i] .* hydro.T[:,i] / mi[i]

		hydro.F1D[:,3*(i-1)+1] = hydro.rho[:,i] .* hydro.u[:,i]
		hydro.F1D[:,3*(i-1)+2] = hydro.rho[:,i] .* hydro.u[:,i].^2 + hydro.p[:,i]
		hydro.F1D[:,3*(i-1)+3] = hydro.rho[:,i] .* hydro.u[:,i] .*
            UU[:,3*(i-1)+3] ./ hydro.rho[:,i] + hydro.p[:,i] .* hydro.u[:,i]

		if geom=="spherical"  # define correcting matrix
			hydro.C1D[:,3*(i-1)+1] = 0.
			hydro.C1D[:,3*(i-1)+2] = copy(hydro.p[:,i])
			hydro.C1D[:,3*(i-1)+3] = 0.
		end

	    # electrons: quasi-neutrality + zero-current equation
	    hydro.rho[:,nspec+1] = hydro.rho[:,nspec+1] + me * Zi[i] * hydro.rho[:,i] / mi[i]  # quasi-neutrality
	    hydro.u[:,nspec+1] = hydro.u[:,nspec+1] +
                Zi[i] * hydro.rho[:,i] .* hydro.u[:,i] / mi[i]  # zero-current condition

		hydro.ni[:,i] = hydro.rho[:,i] / mi[i]
		hydro.ne = hydro.ne + Zi[i] * hydro.ni[:,i]
	end
	hydro.u[:,nspec+1] = hydro.u[:,nspec+1] ./ hydro.ne # correct with ne

	# electron temperature and pressure
	hydro.T[:,nspec+1] = (g-1) * me * (UU[:,neqi+1] ./ hydro.rho[:,nspec+1] - 0.5*hydro.u[:,nspec+1].^2 )
	hydro.p[:,nspec+1] = hydro.rho[:,nspec+1] .* hydro.T[:,nspec+1] / me
	hydro.F1D[:,neqi+1] = hydro.u[:,nspec+1] .* UU[:,neqi+1]  + hydro.p[:,nspec+1] .* hydro.u[:,nspec+1]

# 	limit max temperature by Tmax and avoid negative temperatures
	for j = 1:nspec+1
# 		hydro.T[:,j] = min.(hydro.T[:,j], Tmax)
	    hydro.T[:,j] = max.(hydro.T[:,j], 1.e-5*qe_C) # avoid negative temperature
	end

    return hydro
end
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function source_terms!(hydro, k_DT, ke, Qextra, nq)

    Q_DT = zeros(Float64, nz, nspec + 1)
    hydro.G1D[:, :] = 0.0
    q_diff = zeros(Float64, nz)
    q_FS = zeros(Float64, nz)

	if heattransfer_switch & electron_switch
		for i = 1:nspec+1
			for j = 1:nspec+1
				if j != i
					Q_DT[:,i] = Q_DT[:,i] + k_DT[:,i,j] .* (hydro.T[:,j] - hydro.T[:,i])
				end
			end
		end
		Q_DT[:,nspec+1] = Q_DT[:,nspec+1] + Qextra  #electrons
	elseif heattransfer_switch & ~electron_switch   #only ions
		for i = 1:nspec
			for j = 1:nspec
				if j != i
					Q_DT[:,i] = Q_DT[:,i] + k_DT[:,i,j] .* (hydro.T[:,j] - hydro.T[:,i])
				end
			end
		end
	end

	if efield_switch & electron_switch
		hydro.Efield[2:nq-1] = - 0.5 * (hydro.p[3:nq,nspec+1] - hydro.p[1:nq-2,nspec+1]) /
		                  dr ./ (qe_C * hydro.ne[2:nq-1] )
	else
		hydro.Efield[:] = 0.
	end

	for i = 1:nspec
		hydro.G1D[:,3*(i-1)+2] = Zi[i] * qe_C * hydro.ni[:,i] .* hydro.Efield  #source term
		hydro.G1D[:,3*(i-1)+2] = hydro.G1D[:,3*(i-1)+2] + hydro.R1D[:,3*(i-1)+2]
		hydro.G1D[:,3*(i-1)+3] = Q_DT[:,i] + Zi[i] * qe_C * hydro.ni[:,i] .* hydro.Efield .* hydro.u[:,i]
		hydro.G1D[:,3*(i-1)+3] = hydro.G1D[:,3*(i-1)+3] + hydro.R1D[:,3*(i-1)+3]
	end

	q_diff[2:nq-1] = 1. / dr * ke[2:nq-1] .* (hydro.T[3:nq,nspec+1] - hydro.T[2:nq-1,nspec+1])

	vth_e = sqrt.(hydro.T[:,nspec+1] / me)
	q_FS[2:nq-1] = hydro.rho[3:nq,nspec+1] / me .* hydro.T[3:nq,nspec+1] .* vth_e[3:nq]  #ensures that it is forward-differencing later
	q_diff[:] = sign.(q_diff) .* min.(flimit * q_FS, abs.(q_diff))

	hydro.G1D[:,neqi+1] = Q_DT[:,nspec+1] - qe_C * hydro.ne .* hydro.Efield .* hydro.u[:,nspec+1]

	if electron_heat_conduction
		hydro.G1D[2:nq-1,neqi+1] = hydro.G1D[2:nq-1,neqi+1] + 1. / dr * (q_diff[2:nq-1] - q_diff[1:nq-2])
		if geom=="spherical"  #add correction due to spherical geometry
			hydro.G1D[:,neqi+1] = hydro.G1D[:,neqi+1] + 2. * q_diff ./ r
		end
	end

	#Lindl, ch. 3 eq. (26) - rho in g/cm3, T in keV
	if bremsstrahlung
		hydro.G1D[:,neqi+1] = hydro.G1D[:,neqi+1] - 3.0e16 * (1.e3 * me * hydro.ne).^2 .*
			             (1.e-3 * hydro.T[:,nspec+1] / qe_C)
	end

    return hydro
end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
#!/usr/bin/env julia

using StochasticCubicDescent
import Distributions
using Base.Test

srand(42)

function gen_predictors(N, n)
    const auto_precision = -0.3
    const sd_range = 2.

    precision_mat = eye(n)
    for i in 2:(n-1)
        precision_mat[i-1, i] = precision_mat[i+1,i] = auto_precision
    end
    if n > 1
        precision_mat[2, 1] = precision_mat[n - 1, n] = auto_precision
    end
    cov_mat = precision_mat^-1
    cov_mat = (cov_mat + cov_mat') / 2
    mean_vec = [sqrt(cov_mat[j,j]) * ((j - (n/2)) * sd_range / (n / 2)) for j in 1:n]
    design_dist_unrounded = Distributions.MvNormal(mean_vec, cov_mat)

    design_unrounded = Distributions.rand(design_dist_unrounded, N)
    (sign(design_unrounded) + 1) / 2
end

function gen_responses(x, w_true)
    n, N = size(x)
    ret = Array(Float64, N)

    for i in 1:N
        c_i = 1 / (1  + exp(-x[:, i]' * w_true))
        ret[i] = rand() > c_i ? 0.0 : 1.0
    end

    ret
end

function get_logistic_objective(x, z, b, bH, reg)
    n, N = size(x)

    ξ!(w, ξ) = ξ[:] = rand(1:N, b)

    val(w, ξ) = begin
        y = 0.0
        for i in ξ
            c_i = 1 / (1  + exp(-x[:, i]' * w))
            y -= z[i] > 0.5 ? log(c_i) : log(1 - c_i)
        end

        # make these unbiased estimates of the full objective
        y /= length(ξ)

        # a per-example regularizer
        y + 0.5 * reg * (w' * w)
    end

    grad!(w, ξ, gr) = begin
        fill!(gr, 0.0)

        for i in ξ
            c_i = 1 / (1  + exp(-x[:, i]' * w))
            ci_offset = c_i - z[i]
            gr .+= ci_offset * x[:, i]
        end

        # make these unbiased estimates of the full objective
        gr ./= length(ξ)
        gr .+= reg * w
    end

    hv!(w, ξ, v, hv) = begin
        fill!(hv, 0.0)

        for i in ξ
            c_i = 1 / (1  + exp(-x[:, i]' * w))
            ci_ent = c_i * (1 - c_i)
            d_i = ci_ent * (x[:, i]' * v)
            hv .+= d_i * x[:, i]
        end

        hv ./= length(ξ)
        hv .+= reg * v
    end

    ObjectiveFunction(val, grad!, hv!, ξ!, 45.)
end

function run()
    const N = 7000
    const n = 40

    const reg = 1e2

    const b = 50
    const bH = 600

    const L = 1
    const M = 1

    x = gen_predictors(N, n)
    w_true = randn(n)
    z = gen_responses(x, w_true)
    f = get_logistic_objective(x, z, b, bH, reg)

    # verifying gradient
    grad1 = zeros(n)
    f1 = f.val(w_true, 44:55)
    f.grad!(w_true, 44:55, grad1)

    grad2 = zeros(n)
    f2 = f.val(w_true + 1e-5, 44:55)
    f.grad!(w_true + 1e-5, 44:55, grad2)

    finite_diff_deriv = (f2 - f1) / 1e-5
    @test sum(grad1) ≈ finite_diff_deriv atol=1e-1

    # verifying hessian
    finite_diff_hv = (grad2 - grad1) / 1e-5
    hv1 = zeros(n)
    f.hv!(w_true, 44:55, ones(n), hv1)
    @test hv1 ≈ finite_diff_hv atol=1e-1
end

run()


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    esomTrain(data::AbstractMatrix{Float64})

Train an ESOM for the given data set.

Uses on-line learning per default. (Batch is not yet implemented)
"""
function esomTrain(data::AbstractMatrix{Float64}, settings::Settings = defaultSettings)
    return esomTrainOnline(data, settings)
end

"""
    esomTrainOnline(data::AbstractMatrix{Float64})

Train an ESOM on-line for the given data set.
"""
function esomTrainOnline(data::AbstractMatrix{Float64}, settings::Settings = defaultSettings)
    weights = esomInit(data, settings)
    return esomTrainOnline!(data, weights, settings)
end

"""
    esomTrainOnline!(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64})

Train the ESOM weights on-line for the given data set.

This function mutates `weights`
"""
function esomTrainOnline!(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                         settings::Settings = defaultSettings)
    @assert size(data, 2) == size(weights, 1)
    s = settings
    coolDownRadius = coolDown(s.radiusCooling, s.radius, s.epochs)
    coolDownLearningrate = coolDown(s.learningRateCooling, s.learningRate, s.epochs)
    for i in 1:s.epochs
        radius = coolDownRadius(i)
        learningRate = coolDownLearningrate(i)
        esomTrainEpoch!(data, weights, radius, learningRate, settings)
    end
    if settings.shiftToHighestDensity && settings.toroid
        weights = shiftToHighestDensity(data, weights)
    end
    return weights
end

"""
    esomInit(data::AbstractMatrix{Float64}, settings::Settings = defaultSettings)

Initilize an ESOM with values based on `settings.initMethod`.
"""
function esomInit(data::AbstractMatrix{Float64}, settings::Settings = defaultSettings)
    f = initMethod(settings)
    result::Matrix{Float64} = hcat(f.(Slices(data, 1))...)
    return reshape(permutedims(result), (size(data,2), settings.latticeSize...))
end

"""
    esomTrainEpoch!(data, weights, radius, learningRate)

Train the ESOM for a single epoch.

This function mutates `weights`
"""
function esomTrainEpoch!(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                         radius::Float64, learningRate::Float64,
                         settings::Settings = defaultSettings)
    data_view = view(data, randperm(size(data, 1)), :)
    for dataPoint in eachrow(data_view)
        esomTrainStep!(dataPoint, weights, radius, learningRate, settings)
    end
    return weights
end

"""
    esomTrainStep!(dataPoint, weights, radius, learningRate)

Train the ESOM with a single data point.

This function mutates `weights`
"""
function esomTrainStep!(dataPoint::AbstractVector{Float64}, weights::EsomWeights{Float64},
                        radius::Float64, learningRate::Float64,
                        settings::Settings = defaultSettings)
    @assert size(dataPoint, 1) == size(weights, 1)
    bestMatch_index = bestMatch(dataPoint, weights, settings)
    neighbours = neighbourhood(bestMatch_index, radius, settings)
    kernel = neighbourhoodKernel(settings.neighbourhoodKernel, radius)
    distances = map(i -> latticeDistance(i, bestMatch_index, settings), neighbours)
    @. weights[:,neighbours] += (learningRate * kernel(distances))' *
                             (dataPoint - weights[:, neighbours])
    return weights
end

"""
    projection(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64})

Generate a projection from each data point to the best matching ESOM neuron.
"""
function projection(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                    settings::Settings = defaultSettings;
                    key::AbstractVector{Int} = 1:size(data,1))
    @assert length(key) == size(data, 1)
    @assert allunique(key)
    @assert all(i -> 1 <= i <= size(data, 1), key)
    return (i -> i => bestMatch(data[i,:], weights)).(key) |> Dict
end

# compatability with DataIO.LRNData
function projection(data::LRNData, args...; kwargs...)
    projection(data.data, args...; key=data.key, kwargs...)
end

"""
    bestMatch(dataPoint::AbstractVector{Float64}, weights::EsomWeights{Float64})

Search for the best matching ESOM neuron for the given data point.
"""
function bestMatch(dataPoint::AbstractVector{Float64}, weights::EsomWeights{Float64},
                   settings::Settings = defaultSettings)
    @assert size(dataPoint, 1) == size(weights, 1)
    dist = settings.distance
    slice = Slices(weights, 1)
    index::CartesianIndex{2} = _findmin(weight -> dist(weight, dataPoint), slice)[2]
    return index
end

"""
    shiftWeights(weights::EsomWeights{Float64}, pos::CartesianIndex{2})

Shift the ESOM so the given `pos` is in the middle.

This makes only sense for a toroidal map.
"""
function shiftWeights(weights::EsomWeights{Float64}, pos::CartesianIndex{2},
                      settings::Settings = defaultSettings)
    # since the plot shos the matrix four times, the midpoint of the plot is
    # equal to the latticeSize
    offset = CartesianIndex(settings.latticeSize...) - pos
    indices = CartesianIndices(Slices(weights, 1))
    new_indices = (indices .- offset)[:]
    new_indices = wrapCoordsOnToroid(new_indices)
    result = similar(weights)
    for (old, new) in zip(indices, new_indices)
        result[:,new] = weights[:,old]
    end
    return result
end

"""
    shiftToHighestDensity(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64})

Shift the ESOM so the point of highest density is centered.
"""
function shiftToHighestDensity(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                               settings = defaultSettings)
    if !settings.toroid return weights end
    radius = filter(!iszero, pairwise(Euclidean(), data, dims=1)) |> mean
    p = pmatrix(data, weights, settings, radius = radius)
    pos = findfirst(isequal(maximum(p)), p)
    return shiftWeights(weights, pos, settings)
end

# compatability with DataIO.LRNData
for f in (:esomTrain, :esomTrainOnline, :esomTrainOnline!, :esomInit,
          :shiftToHighestDensity)
    @eval @inline ($f)(data::LRNData, args...; kwargs...) = ($f)(data.data, args...; kwargs...)
end


text:
using Jecco
using Jecco.AdS5_3_1
using FFTW
using Test
include("../scripts/test_utils_2.jl")

#We are going to manipulate an analytic configuration and test the output of create_new_initial. Perturbed conformal BB.
#Miguel and Thanasis: if you change the arguments of the function create_new_initial
#you'll have to change line 75 of this script.
dirname = "/home/mikel/Documentos/Jecco.jl/data/test/"
outdir  = "/home/mikel/Documentos/Jecco.jl/data/test_new/"

grid = SpecCartGrid3D(
    x_min            = -15.0,
    x_max            =  15.0,
    x_nodes          =  32,
    y_min            = -15.0,
    y_max            =  15.0,
    y_nodes          =  32,
    u_outer_min      =  0.1,
    u_outer_max      =  1.005,
    u_outer_domains  =  1,
    u_outer_nodes    =  64,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

new_grid = SpecCartGrid3D(
    x_min            = -30.0,
    x_max            =  30.0,
    x_nodes          =  32,
    y_min            = -30.0,
    y_max            =  30.0,
    y_nodes          =  32,
    u_outer_min      =  0.1,
    u_outer_max      =  1.005,
    u_outer_domains  =  1,
    u_outer_nodes    =  64,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

parameters_collision =AdS5_3_1.new_parameters_coll(
    dirname1  = dirname,
    dirname2  = dirname,
    x1_center = -15.0,
    y1_center = 0.0,
    x2_center = 15.0,
    y2_center = 0.0,
    fx21      = -0.1,
    fy21      = 0.0,
    fx22      = 0.1,
    fy22      = 0.0,
    u_AH      = 1.0,
)

analytic = analytic_parameters(
    energy   = 1.0,
    phi0     = 1.0,
    u_AH     = 1.0,
    a4_amp   = 0.01,
    phi2_amp = 0.01,
    b14_amp  = 0.01,
    b24_amp  = 0.01,
)

tolerance = 5.e-15 #tolerance in the error of our routine.

potential = Phi8Potential(oophiM2=-1.0, oophiQ=0.1,)#It doesn't matter for the purpose of changing a config. it will matter for the gauge change
io = InOut(out_dir = dirname, checkpoint_dir = dirname, out_boundary_every=1,
                  out_gauge_every=1,out_bulk_every=1,remove_existing = true,)

#We generate the data, and the file according to analytic expressions, and modify the configuration analytically
boundary_new, gauge_new, bulkevols_new = test_design_collision(grid,new_grid,analytic,parameters_collision,io,potential)

# We call the function here
io = InOut(recover_dir = dirname, out_dir = outdir, checkpoint_dir = outdir,
               out_boundary_every=1,out_gauge_every=1,out_bulk_every=1,remove_existing = true,)

AdS5_3_1.design_collision(new_grid, io, parameters_collision)

#We compare the analytical result with the routine
Nsys = grid.u_outer_domains+1
for i in 1:Nsys
    B1_aux  = BulkTimeSeries(outdir,:B1,i)
    B2_aux  = BulkTimeSeries(outdir,:B2,i)
    G_aux   = BulkTimeSeries(outdir,:G,i)
    phi_aux = BulkTimeSeries(outdir,:phi,i)

    difB1  = maximum(abs.(B1_aux[end,:,:,:]-bulkevols_new[i].B1))
    difB2  = maximum(abs.(B2_aux[end,:,:,:]-bulkevols_new[i].B2))
    difG   = maximum(abs.(G_aux[end,:,:,:]-bulkevols_new[i].G))
    difphi = maximum(abs.(phi_aux[end,:,:,:]-bulkevols_new[i].phi))

    @testset "Bulk grid $i" begin
        @test difB1  < tolerance
        @test difB2  < tolerance
        @test difG   < tolerance
        @test difphi < tolerance
    end
end

a4_aux  = BoundaryTimeSeries(outdir,:a4)
fx2_aux = BoundaryTimeSeries(outdir,:fx2)
fy2_aux = BoundaryTimeSeries(outdir,:fy2)

difa   = maximum(abs.(a4_aux[end,:,:]-boundary_new.a4[1,:,:]))
diffx2 = maximum(abs.(fx2_aux[end,:,:]-boundary_new.fx2[1,:,:]))
diffy2 = maximum(abs.(fy2_aux[end,:,:]-boundary_new.fy2[1,:,:]))

@testset "Boundary" begin
    @test difa   < tolerance
    @test diffx2 < tolerance
    @test diffy2 < tolerance
end

xi    = XiTimeSeries(outdir)
difxi = maximum(abs.(xi[end,:,:]-gauge_new.xi[1,:,:]))

@testset "Gauge" begin
    @test difxi < tolerance
end
nothing


text:
import Base: isempty,
             rand,
             convert

export HPolyhedron,
       dim, σ, ∈,
       addconstraint!,
       constraints_list,
       tohrep, tovrep,
       convex_hull,
       cartesian_product,
       vertices_list,
       singleton_list,
       isempty,
       remove_redundant_constraints,
       remove_redundant_constraints!,
       constrained_dimensions

"""
    HPolyhedron{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a convex polyhedron in H-representation.

### Fields

- `constraints` -- vector of linear constraints
"""
struct HPolyhedron{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
    constraints::Vector{LinearConstraint{N, VN}}

    function HPolyhedron(constraints::Vector{LinearConstraint{N, VN}}) where {N<:Real,
                                                                              VN<:AbstractVector{N}}
        return new{N, VN}(constraints)
    end
end

isoperationtype(::Type{<:HPolyhedron}) = false
isconvextype(::Type{<:HPolyhedron}) = true

# constructor for an HPolyhedron with no constraints
function HPolyhedron{N, VN}() where {N<:Real, VN<:AbstractVector{N}}
    HPolyhedron(Vector{LinearConstraint{N, VN}}())
end

# constructor for an HPolyhedron with no constraints and given numeric type
function HPolyhedron{N}() where {N<:Real}
    HPolyhedron(Vector{LinearConstraint{N, Vector{N}}}())
end

# constructor for an HPolyhedron without explicit numeric type, defaults to Float64
function HPolyhedron()
    HPolyhedron{Float64}()
end

# constructor from a simple H-representation
HPolyhedron(A::AbstractMatrix{N}, b::AbstractVector{N}) where {N<:Real} =
    HPolyhedron(constraints_list(A, b))

# convenience union type
const HPoly{N} = Union{HPolytope{N}, HPolyhedron{N}}

# --- LazySet interface functions ---

"""
    dim(P::HPoly{N}) where {N<:Real}

Return the dimension of a polyhedron in H-representation.

### Input

- `P`  -- polyhedron in H-representation

### Output

The ambient dimension of the polyhedron in H-representation.
If it has no constraints, the result is ``-1``.
"""
function dim(P::HPoly{N}) where {N<:Real}
    return length(P.constraints) == 0 ? -1 : length(P.constraints[1].a)
end

"""
    ρ(d::AbstractVector{N}, P::HPoly{N};
      solver=default_lp_solver(N)) where {N<:Real}

Evaluate the support function of a polyhedron (in H-representation) in a given
direction.

### Input

- `d`      -- direction
- `P`      -- polyhedron in H-representation
- `solver` -- (optional, default: `default_lp_solver(N)`) the backend used to
              solve the linear program

### Output

The support function of the polyhedron.
If a polytope is unbounded in the given direction, we throw an error.
If a polyhedron is unbounded in the given direction, the result is `Inf`.
"""
function ρ(d::AbstractVector{N}, P::HPoly{N};
           solver=default_lp_solver(N)) where {N<:Real}
    lp, unbounded = σ_helper(d, P, solver)
    if unbounded
        if P isa HPolytope
            error("the support function in direction $(d) is undefined " *
                  "because the polytope is unbounded")
        end
        return N(Inf)
    end
    return dot(d, lp.sol)
end

"""
    σ(d::AbstractVector{N}, P::HPoly{N}; solver=default_lp_solver(N)
     ) where {N<:Real}

Return the support vector of a polyhedron (in H-representation) in a given
direction.

### Input

- `d`      -- direction
- `P`      -- polyhedron in H-representation
- `solver` -- (optional, default: `default_lp_solver(N)`) the backend used to
              solve the linear program

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector{N}, P::HPoly{N}; solver=default_lp_solver(N)
          ) where {N<:Real}
    lp, unbounded = σ_helper(d, P, solver)
    if unbounded
        if P isa HPolytope
            error("the support vector in direction $(d) is undefined because " *
                  "the polytope is unbounded")
        end
        # construct the solution from the solver's ray result
        if lp == nothing
            ray = d
        elseif haskey(lp.attrs, :unboundedray)
            ray = lp.attrs[:unboundedray]
        else
            error("LP solver did not return an infeasibility ray")
        end
        res = Vector{N}(undef, length(ray))
        @inbounds for i in 1:length(ray)
            if ray[i] == zero(N)
                res[i] = zero(N)
            elseif ray[i] > zero(N)
                res[i] = N(Inf)
            else
                res[i] = N(-Inf)
            end
        end
        return res
    else
        return lp.sol
    end
end

function σ_helper(d::AbstractVector{N}, P::HPoly{N}, solver) where {N<:Real}
    # let c = -d as a Vector since GLPK does not accept sparse vectors
    # (see #1011)
    c = to_negative_vector(d)

    (A, b) = tosimplehrep(P)
    if length(b) == 0
        unbounded = true
        lp = nothing
    else
        sense = '<'
        l = -Inf
        u = Inf
        lp = linprog(c, A, sense, b, l, u, solver)
        if lp.status == :Unbounded
            unbounded = true
        elseif lp.status == :Infeasible
            error("the support vector is undefined because the polyhedron is " *
                  "empty")
        else
            unbounded = false
        end
    end
    return (lp, unbounded)
end

"""
    isbounded(P::HPolyhedron)

Determine whether a polyhedron in constraint representation is bounded.

### Input

- `P` -- polyhedron in constraint representation

### Output

`true` iff the polyhedron is bounded.

### Algorithm

We first check if the polyhedron has more than `max(dim(P), 1)` constraints,
which is a necessary condition for boundedness.
If so, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(P::HPolyhedron)
    if length(P.constraints) <= max(dim(P), 1)
        return false
    end
    return isbounded_unit_dimensions(P)
end

"""
    rand(::Type{HPolyhedron}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a polyhedron.

### Input

- `HPolyhedron` -- type for dispatch
- `N`           -- (optional, default: `Float64`) numeric type
- `dim`         -- (optional, default: 2) dimension (is ignored)
- `rng`         -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`        -- (optional, default: `nothing`) seed for reseeding

### Output

A polyhedron.

### Algorithm

We first create a random polytope and then randomly remove some of the
constraints.
"""
function rand(::Type{HPolyhedron};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    P = rand(HPolytope; N=N, dim=dim, rng=rng)
    constraints_P = constraints_list(P)
    constraints_Q = Vector{eltype(constraints_P)}()
    for i in 1:length(constraints_P)
        if rand(Bool)
            push!(constraints_Q, constraints_P[i])
        end
    end
    return HPolyhedron(constraints_Q)
end

# ===========================================
# HPolyhedron and HPolytope's shared methods
# ===========================================

"""
    addconstraint!(P::HPoly{N}, constraint::LinearConstraint{N}) where {N<:Real}

Add a linear constraint to a polyhedron in H-representation.

### Input

- `P`          -- polyhedron in H-representation
- `constraint` -- linear constraint to add

### Notes

It is left to the user to guarantee that the dimension of all linear constraints
is the same.
"""
function addconstraint!(P::HPoly{N},
                        constraint::LinearConstraint{N}) where {N<:Real}
    push!(P.constraints, constraint)
    return nothing
end

"""
    constraints_list(P::HPoly{N}) where {N<:Real}

Return the list of constraints defining a polyhedron in H-representation.

### Input

- `P` -- polyhedron in H-representation

### Output

The list of constraints of the polyhedron.
"""
function constraints_list(P::HPoly{N}) where {N<:Real}
    return P.constraints
end

"""
    tohrep(P::HPoly{N}) where {N<:Real}

Return a constraint representation of the given polyhedron in constraint
representation (no-op).

### Input

- `P` -- polyhedron in constraint representation

### Output

The same polyhedron instance.
"""
function tohrep(P::HPoly{N}) where {N<:Real}
    return P
end

"""
    normalize(P::HPoly{N}, p=N(2)) where {N<:Real}

Normalize a polyhedron in constraint representation.

### Input

- `P` -- polyhedron in constraint representation
- `p` -- (optional, default: `2`) norm

### Output

A new polyhedron in constraint representation whose normal directions ``a_i``
are normalized, i.e., such that ``‖a_i‖_p = 1`` holds.
"""
function normalize(P::HPoly{N}, p=N(2)) where {N<:Real}
    constraints = [normalize(hs, p) for hs in constraints_list(P)]
    T = basetype(P)
    return T(constraints)
end

"""
    remove_redundant_constraints(P::HPoly{N};
                                 backend=default_lp_solver(N)
                                ) where {N<:Real}

Remove the redundant constraints in a polyhedron in H-representation.

### Input

- `P`       -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the backend used to
               solve the linear program

### Output

A polyhedron equivalent to `P` but with no redundant constraints, or an empty
set if `P` is detected to be empty, which may happen if the constraints are
infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints(P::HPoly{N};
                                      backend=default_lp_solver(N)
                                     ) where {N<:Real}
    Pred = copy(P)
    if remove_redundant_constraints!(Pred, backend=backend)
        return Pred
    else # the polyhedron P is empty
        return EmptySet{N}(dim(P))
    end
end

"""
    remove_redundant_constraints!(P::HPoly{N};
                                  backend=default_lp_solver(N)) where {N<:Real}

Remove the redundant constraints in a polyhedron in H-representation; the
polyhedron is updated in-place.

### Input

- `P`       -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the backend used to
               solve the linear program

### Output

`true` if the method was successful and the polyhedron `P` is modified by
removing its redundant constraints, and `false` if `P` is detected to be empty,
which may happen if the constraints are infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints!(P::HPoly{N};
                                       backend=default_lp_solver(N)
                                      ) where {N<:Real}
    remove_redundant_constraints!(P.constraints, backend=backend)
end

"""
    translate(P::HPoly{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:Real}

Translate (i.e., shift) a polyhedron in constraint representation by a given
vector.

### Input

- `P`     -- polyhedron in constraint representation
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated polyhedron in constraint representation.

### Notes

The normal vectors of the constraints (vector `a` in `a⋅x ≤ b`) are shared with
the original constraints if `share == true`.

### Algorithm

We translate every constraint.
"""
function translate(P::HPoly{N}, v::AbstractVector{N}; share::Bool=false
                  ) where {N<:Real}
    @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    constraints = [translate(c, v; share=share) for c in constraints_list(P)]
    T = basetype(P)
    return T(constraints)
end

# ========================================================
# External methods that require Polyhedra.jl to be loaded
# ========================================================

"""
    convex_hull(P1::HPoly{N}, P2::HPoly{N};
               [backend]=default_polyhedra_backend(P1, N)) where {N<:Real}

Compute the convex hull of the set union of two polyhedra in H-representation.

### Input

- `P1`         -- polyhedron
- `P2`         -- another polyhedron
- `backend`    -- (optional, default: `default_polyhedra_backend(P1, N)`)
                  the polyhedral computations backend

### Output

The `HPolyhedron` (resp. `HPolytope`) obtained by the concrete convex hull of
`P1` and `P2`.

### Notes

For performance reasons, it is suggested to use the `CDDLib.Library()` backend
for the `convex_hull`.

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function convex_hull(P1::HPoly{N},
                     P2::HPoly{N};
                     backend=default_polyhedra_backend(P1, N)) where {N<:Real}
    require(:Polyhedra; fun_name="convex_hull")
    Pch = Polyhedra.convexhull(polyhedron(P1; backend=backend),
                               polyhedron(P2; backend=backend))
    removehredundancy!(Pch)
    return convert(basetype(P1), Pch)
end

"""
    cartesian_product(P1::HPoly{N}, P2::HPoly{N};
                      [backend]=default_polyhedra_backend(P1, N)
                     ) where {N<:Real}

Compute the Cartesian product of two polyhedra in H-representaion.

### Input

- `P1`         -- polyhedron
- `P2`         -- another polyhedron
- `backend`    -- (optional, default: `default_polyhedra_backend(P1, N)`)
                  the polyhedral computations backend

### Output

The polyhedron obtained by the concrete cartesian product of `P1` and `P2`.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function cartesian_product(P1::HPoly{N},
                           P2::HPoly{N};
                           backend=default_polyhedra_backend(P1, N)
                          ) where {N<:Real}
    require(:Polyhedra; fun_name="`cartesian_product")
    Pcp = Polyhedra.hcartesianproduct(polyhedron(P1; backend=backend),
                                      polyhedron(P2; backend=backend))
    return convert(basetype(P1), Pcp)
end

"""
    tovrep(P::HPoly{N};
          [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Transform a polyhedron in H-representation to a polytope in V-representation.

### Input

- `P`       -- polyhedron in constraint representation
- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`) the
               backend for polyhedral computations

### Output

The `VPolytope` which is the vertex representation of the given polyhedron
in constraint representation.

### Notes

The conversion may not preserve the numeric type (e.g., with `N == Float32`)
depending on the backend.
For further information on the supported backends see [Polyhedra's
documentation](https://juliapolyhedra.github.io/).
"""
function tovrep(P::HPoly{N};
                backend=default_polyhedra_backend(P, N)) where {N<:Real}
    require(:Polyhedra; fun_name="tovrep")
    P = polyhedron(P; backend=backend)
    return VPolytope(P)
end

"""
    vertices_list(P::HPolyhedron{N}) where {N<:Real}

Return the list of vertices of a polyhedron in constraint representation.

### Input

- `P` -- polyhedron in constraint representation

### Output

This function returns an error because the polyhedron is possibly unbounded.
If `P` is known to be bounded, try converting to `HPolytope` first:

```jldoctest
julia> P = HPolyhedron([HalfSpace([1.0, 0.0], 1.0),
                        HalfSpace([0.0, 1.0], 1.0),
                        HalfSpace([-1.0, 0.0], 1.0),
                        HalfSpace([0.0, -1.0], 1.0)]);

julia> P_as_polytope = convert(HPolytope, P);
```
"""
function vertices_list(P::HPolyhedron{N}) where {N<:Real}
    throw(ArgumentError("the list of vertices of a (possibly unbounded) " *
        "polyhedron is not defined; if the polyhedron is bounded, try " *
        "converting to `HPolytope` first"))
end

"""
    singleton_list(P::HPolyhedron{N}) where {N<:Real}

Return the vertices of a polyhedron in H-representation as a list of singletons.

### Input

- `P` -- polytope in constraint representation

### Output

This function returns an error because the polyhedron is possibly unbounded.
If `P` is known to be bounded, try converting to `HPolytope` first.
"""
function singleton_list(P::HPolyhedron{N}) where {N<:Real}
    throw(ArgumentError("the list of singletons of a (possibly unbounded) " *
        "polyhedron is not defined; if the polyhedron is bounded, try " *
        "converting to `HPolytope` first"))
end

"""
   isempty(P::HPoly{N}, witness::Bool=false;
           [use_polyhedra_interface]::Bool=false, [solver]=default_lp_solver(N),
           [backend]=nothing) where {N<:Real}

Determine whether a polyhedron is empty.

### Input

- `P`       -- polyhedron
- `witness` -- (optional, default: `false`) compute a witness if activated
- `use_polyhedra_interface` -- (optional, default: `false`) if `true`, we use
               the `Polyhedra` interface for the emptiness test
- `solver`  -- (optional, default: `default_lp_solver(N)`) LP-solver backend
- `backend` -- (optional, default: `nothing`) backend for polyhedral
               computations in `Polyhedra`; its value is set internally (see the
               Notes below for details)

### Output

* If `witness` option is deactivated: `true` iff ``P = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``P = ∅``
  * `(false, v)` iff ``P ≠ ∅`` and ``v ∈ P``

### Notes

The default value of the `backend` is set internally and depends on whether the
`use_polyhedra_interface` option is set or not.
If the option is set, we use `default_polyhedra_backend(P, N)`.

Witness production is not supported if `use_polyhedra_interface` is `true`.

### Algorithm

The algorithm sets up a feasibility LP for the constraints of `P`.
If `use_polyhedra_interface` is `true`, we call `Polyhedra.isempty`.
Otherwise, we set up the LP internally.
"""
function isempty(P::HPoly{N},
                 witness::Bool=false;
                 use_polyhedra_interface::Bool=false,
                 solver=default_lp_solver(N),
                 backend=nothing) where {N<:Real}
    if length(constraints_list(P)) < 2
        # catch corner case because of problems in LP solver for Rationals
        return witness ? (false, an_element(P)) : false
    end
    if use_polyhedra_interface
        require(:Polyhedra; fun_name="isempty", explanation="with the active " *
            "option `use_polyhedra_interface`")
        if backend == nothing
            backend = default_polyhedra_backend(P, N)
        end
        result = Polyhedra.isempty(polyhedron(P; backend=backend), solver)
        if result
            return witness ? (true, N[]) : true
        elseif witness
            error("witness production is not supported yet")
        else
            return false
        end
    else
        A, b = tosimplehrep(P)
        lbounds, ubounds = -Inf, Inf
        sense = '<'
        obj = zeros(N, size(A, 2))
        lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
        if lp.status == :Optimal
            return witness ? (false, lp.sol) : false
        elseif lp.status == :Infeasible
            return witness ? (true, N[]) : true
        end
        error("LP returned status $(lp.status) unexpectedly")
    end
end

convert(::Type{HPolytope}, P::HPolyhedron{N}) where {N<:Real} =
    HPolytope(copy(constraints_list(P)))
convert(::Type{HPolyhedron}, P::HPolytope{N}) where {N<:Real} =
    HPolyhedron(copy(constraints_list(P)))

# ==========================================
# Lower level methods that use Polyhedra.jl
# ==========================================

function load_polyhedra_hpolyhedron() # function to be loaded by Requires
return quote
# see the interface file AbstractPolytope.jl for the imports

function convert(::Type{HPolyhedron}, P::HRep{N}) where {N}
    VN = Polyhedra.hvectortype(P)
    constraints = Vector{LinearConstraint{N, VN}}()
    for hi in Polyhedra.allhalfspaces(P)
        a, b = hi.a, hi.β
        if isapproxzero(norm(a))
            continue
        end
        push!(constraints, HalfSpace(a, b))
    end
    return HPolyhedron(constraints)
end

"""
     HPolyhedron(P::HRep{N}) where {N}

Return a polyhedron in H-representation given a `HRep` polyhedron
from `Polyhedra.jl`.

### Input

- `P` -- `HRep` polyhedron

### Output

An `HPolyhedron`.
"""
function HPolyhedron(P::HRep{N}) where {N}
    convert(HPolyhedron, P)
end

"""
    polyhedron(P::HPoly{N};
               [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Return an `HRep` polyhedron from `Polyhedra.jl` given a polytope in
H-representation.

### Input

- `P`       -- polytope
- `backend` -- (optional, default: call `default_polyhedra_backend(P, N)`)
                the polyhedral computations backend

### Output

An `HRep` polyhedron.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function polyhedron(P::HPoly{N};
                    backend=default_polyhedra_backend(P, N)) where {N<:Real}
    A, b = tosimplehrep(P)
    return Polyhedra.polyhedron(Polyhedra.hrep(A, b), backend)
end

end # quote
end # function load_polyhedra_hpolyhedron()


text:
module jPCA

export fit, JPCA

using MultivariateStats,Optim

include("reshape_skew.jl")
include("getRealVs.jl")

import MultivariateStats.fit

type JPCA
    pca
    M
    Mskew
    proj
end

"""
    fit(jpca::Type{JPCA}, X, k=10 or num features)
NB: X should be a time X features array.
"""
function fit(jpca::Type{JPCA}, X, k=size(X,2) < 10 ? size(X,2) : 10)
    size(X,2) > size(X,1) && error("Must have more rows than columns")
    k > size(X,2) && error("k must be <= the number of columns in X")
	isodd(k) && error("k must be even")

    pca   = fit(PCA, X, pratio=1., maxoutdim=k)
    X2    = pca.proj[:,1:k]

    T     = (1:size(X2,1))[2:end]
    Tpre  = (1:size(X2,1))[1:end-1]
    dX2 = X2[T,:].-X2[Tpre,:]
    M = dX2/X2[Tpre,:]

    Mskew = skewsymregress(dX2, X2[Tpre,:])
    D,V   = eig(Mskew)
    V     .= V[:,sortperm(abs.(D))]

    jPCs = zeros(size(V)...)
    for pair in 1:div(k,2)
		vi1 = 1+2(pair-1)
		vi2 = 2pair

		VconjPair  = V[:, [vi1,vi2]]
		evConjPair = D[[vi1, vi2]]
		jPCs[:,vi1] .= getRealVs(VconjPair, evConjPair)[1]
		jPCs[:,vi2] .= getRealVs(VconjPair, evConjPair)[2]
    end
	proj = pca.proj * jPCs

    JPCA(pca, M, Mskew, proj)
end

function skewsymregress(dX,X)
    M0 = X\dX
    M0 .= 0.5*(M0 - M0')
    m0 = reshape_skew(M0)

    function skewsymderiv(m)
        # Evaluate objective function
        f = vecnorm( dX - X*reshape_skew(m))^2;

        # Evaluate derivative
        # D = (dX - X*reshape_skew(m))'*X;
        # df = 2*reshape_skew(D - D')
        f
    end

    M = optimize(skewsymderiv, m0)
    return reshape_skew(M.minimizer)
end

end

max_stars_repo_name: alexmorley/jPCA

max_issues_repo_name: alexmorley/jPCA

max_forks_repo_name: alexmorley/jPCA


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

using Documenter, ZhuyinPinyin

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

deploydocs(
    repo = "github.com/tmthyln/ZhuyinPinyin.jl.git",
    devbranch = "main",
    devurl="latest"
    )


text:
function model_tyc3559_norv(time,param)
# Input parameters (x=param) are:
# x[1] = P  (units of day)
# x[2] = inc = inclination angle
# x[3] = p = R_p/R_* = radius of planet in units of radius of star
# x[4] = t0 = mid-point of transit
# x[5] = u1 = linear limb-darkening coefficient
# x[6] = u2 = quadratic limb-darkening coefficient
# x[7] = f0 = uneclipsed flux
# x[8] = a/R_* = semi-major axis divided by R_*
# x[9] = e = eccentricity
# x[10] = omega = longitude of pericentre
# x[11]= a1c
# x[12]= a1s
# x[13]= a2c
# x[14]= a2s
# x[15]= ratio of flux of 3rd star to total flux
# x[16]= fraction of EB flux due to 2nd star
nt = length(time)
model=zeros(nt)
dt = median(time[2:nt]-time[1:nt-1])
# Photometric data points have positive times, while
# RV data points have negative times:
period=param[1]
t0=param[4]
a1c=param[11]
a1s=param[12]
a2c=param[13]
a2s=param[14]
param2=zeros(10)
for i=1:10
  param2[i]=param[i]
end
param2[3]=1.0/param2[3]
param2[8]=param2[8]*param2[3]
param2[7]=1.0
param2[4]=param2[4]+param2[1]/2.0
param2[10]=param2[10]+180.0 # (for now we're assuming circular, but if we add eccentricity, we'll need this)
param1=zeros(10)
for i=1:10
  param1[i] = param[i]
end
param1[7] = 1.0
tdur = 2*(1+param[3])*period/(2pi*param[8])
#println("T_dur: ",tdur)
#print,'Impact parameter: ',cos(param[2]*pi/180.)*param[8]
nsub = 10
time_phot1=0.0
model_phot2 =0.0
for i=1:nt
  time_phot=time[i]
  if abs(mod(time_phot-t0+period/2,period)-period/2) < (0.51*tdur)
    model_phot1 = transit_orb(time_phot,param1,dt,nsub)
  else
    model_phot1 = 1.0
  end
  model_phot1 += a1c*cos((time_phot-t0)*2.0*pi/period)+a1s*sin((time_phot-t0)*2.0*pi/period)
  model_phot1 += a2c*cos((time_phot-t0)*4.0*pi/period)
  if abs(mod(time_phot-t0,period)-period/2) < (0.51*tdur)
    model_phot2 = transit_orb(time_phot,param2,dt,nsub)
  else
    model_phot2 = 1.0
  end
  model_phot2 += a1s*sin((time_phot-t0)*2.0*pi/period)
  model[i]=((model_phot1*(1.0-param[16])+model_phot2*param[16])* (1.0-param[15])+param[15])*param[7]
end
#if(total(finite(model) eq 0) gt 0) then begin
#   print,'Model is non-finite: ',where(finite(model) eq 0)
#   char=get_kbrd(1)
#endif
return model
end


text:
"""
This function reduces the input tasks to a single task by applying `func` to each input task pair.
"""
function reducetasks(func, tasks)
    @assert !isempty(tasks)
    todo = collect(tasks)
    while length(todo) > 1
        a, b = popfirst!(todo), popfirst!(todo)
        newtask = Threads.@spawn func(fetch(a), fetch(b))
        push!(todo, newtask)
    end
    return todo[1]
end

function verify(::Type{SpawnThreads}, sig::Model.Signature, pkhashpairs)
    # kick off workers in parallel
    task_rhside = Threads.@spawn Curve.curve_miller(Config.G2, sig.sig)
    tasks_lhside = [Threads.@spawn(Curve.curve_miller(pair.first.pk, pair.second.hash)) for pair in pkhashpairs]
    task_lhside = reducetasks(*, tasks_lhside)
    res = Curve.field_final_exp(fetch(task_lhside) // fetch(task_rhside))
    return res == one(typeof(res))
end


text:
"""
    Levels(:a => ["yes", "no"])

Return a copy of the table with specified levels and orders for categorical columns
allowing only changing the order of the column.

# Examples

```julia
Levels(:a => ["yes, "no"], :c => [1, 2, 4], :d => ["a", "b", "c"])
Levels("a" => ["yes", "no"], "c" => [1, 2, 4], ordered = ["a", "c"])
Levels(:a => ["yes", "no"], :c => [1, 2, 4], :d => [1, 23, 5, 7], ordered = [:a, :b, :c])
```
"""
struct Levels{K} <: Stateless
  levelspec::K
  ordered::Vector{Symbol}
end

Levels(pairs::Pair{Symbol}...; ordered=Symbol[]) =
  Levels(NamedTuple(pairs), ordered)

Levels(pairs::Pair{K}...; ordered=K[]) where {K<:AbstractString} =
  Levels(NamedTuple(Symbol(k) => v for (k,v) in pairs), Symbol.(ordered))

isrevertible(transform::Levels) = true

# when the col is already a categorical array and wanna change levels and order
_categorify(l::AbstractVector, x::CategoricalVector, o) =
  categorical(x, levels=l, ordered=o), levels(x)

# when the col is normal array and want to change to categorical array
_categorify(l::AbstractVector, x::AbstractVector, o) =
  categorical(x, levels=l, ordered=o), unwrap

# when the col is not need for change or convert back to normal array
_categorify(f::Function, x::AbstractVector, o) =
  o ? (categorical(x, ordered=true), levels(x)) : (f.(x), f)

function apply(transform::Levels, table)
  cols = Tables.columns(table)
  names = Tables.columnnames(cols)

  result = map(names) do nm
    x = Tables.getcolumn(cols, nm)
    l = get(transform.levelspec, nm, identity)
    o = nm ∈ transform.ordered
    _categorify(l, x, o)
  end

  categ = first.(result)
  cache = last.(result)

  𝒯 = (; zip(names, categ)...)
  newtable = 𝒯 |> Tables.materializer(table)

  newtable, cache
end

function revert(transform::Levels, newtable, cache)
  cols = Tables.columns(newtable)
  names = Tables.columnnames(cols)

  ocols = map(zip(cache, names)) do (f, nm)
    x = Tables.getcolumn(cols, nm)
    c, _ = _categorify(f, x, false)
    c
  end

  𝒯 = (; zip(names, ocols)...)
  𝒯 |> Tables.materializer(newtable)
end


text:
using Test
using Random
using UnitTestDesign

### total_combinations
tc_trials = [
    [[2, 2, 2], 2, 12],
    [[2, 3, 2], 2, 16],
    [[2, 3, 2, 2], 2, 6+4+4+6+6+4],
    [[2, 3, 2, 2], 3, 12+12+8+12]
]
for trial in tc_trials
    res = UnitTestDesign.total_combinations(trial[1], trial[2])
    @test res == trial[3]
end

### combination_number

for cn_a in 1:6
    for cn_b in 1:cn_a
        @test UnitTestDesign.combination_number(cn_a, cn_b) == factorial(cn_a) ÷ (factorial(cn_a - cn_b) * factorial(cn_b))
    end
end

### next_multiplicative
arity0 = [3, 7, 9, 4]
vv0 = copy(arity0)
UnitTestDesign.next_multiplicative!(vv0, arity0)
@test vv0 == [1, 1, 1, 1]
vv = [1,1,1]
nmarity = [2,3,2]
for i in 1:prod(nmarity)
    UnitTestDesign.next_multiplicative!(vv, nmarity)
end
@test vv == [1, 1, 1]

### all_combinations
rng = MersenneTwister(9237425)
for ac_trial_idx in 1:5
    n_way = [2, 3, 4][rand(rng, 1:3)]
    param_cnt = rand(rng, (n_way + 1):(n_way + 3))
    arity = rand(rng, 2:4, param_cnt)
    coverage = UnitTestDesign.all_combinations(arity, n_way)
    # It has the right column dimnsion.
    @test size(coverage, 1) == length(arity)
    # Every combination is nonzero.
    @test sum(sum(coverage, dims = 1) == 0) == 0
    # The total number of combinations agrees with expectations.
    @test UnitTestDesign.total_combinations(arity, n_way) == size(coverage, 2)
    # Generate some random combinations and check that they are in there.
    for comb_idx in 1:100
        comb = [rand(rng, 1:arity[cj]) for cj in 1:param_cnt]
        comb[randperm(rng, param_cnt)[1:(param_cnt - n_way)]] .= 0
        @test sum(comb .!= 0) == n_way
        found = false
        for sidx in 1:size(coverage, 2)
            if coverage[:, sidx] == comb
                found = true
            end
        end
        @test found
    end
end

### one_parameter_combinations(arity, n_way)

minimal = UnitTestDesign.one_parameter_combinations([2, 3], 1)
@test minimal == [0 1; 0 2; 0 3]'

paired = UnitTestDesign.one_parameter_combinations([2, 3], 2)
@test paired == [1 1; 2 1; 1 2; 2 2; 1 3; 2 3]'

singled = UnitTestDesign.one_parameter_combinations([2, 2, 3], 2)
@test singled == [
    1 0 1; 2 0 1; 0 1 1; 0 2 1;
    1 0 2; 2 0 2; 0 1 2; 0 2 2;
    1 0 3; 2 0 3; 0 1 3; 0 2 3
    ]'

again = UnitTestDesign.one_parameter_combinations([2, 2, 3], 3)
@test again == [
    1 1 1; 1 2 1; 2 1 1; 2 2 1;
    1 1 2; 1 2 2; 2 1 2; 2 2 2;
    1 1 3; 1 2 3; 2 1 3; 2 2 3
    ]'


text:
#Function to calculate intrinsic mode functions
#N is the maximum number of modes to find
"""
IMF(y; toldev=0.01, tolzero = 0.01, maxorder=4, N=5)
Calculate the intrinsic mode functions of the sequence y along timespan t.
"""
function IMF(y::Array{T,1};
             toldev::Float64   = 0.0001,
             tolzero::Float64  = 0.0001,
             maxorder::Int64   = 4,
             N::Int64          = 5,
             use_fortran::Bool = true,
             verbose::Bool     = false) where T<:Real

    n = length(y)
    t = collect(1:n)
    f = zeros(n,N)

    residue  = zeros(Float64, size(y)...)
    copy!(residue, y)
    eps = 0.00001
    avg = zeros(n)     # pre-allocation for mean envelope
    d   = zeros(n-1,1) # pre-allocation for diff(residue)

    n_modes = 0;
    for i = 1:N # for each IMF
        verbose && println("Computing $i of $N...")
        fill!(avg, 1.0)
        sd = 2*toldev

        while (mean(abs.(avg)) > tolzero && sd > toldev) # do sift?
            sift!(view(f,:,i), avg, residue, y, t, eps, use_fortran)
        end

        copy!(residue, f[:,i])

        n_modes += 1
        # Check to see if it's worth continuing or if the remainder is monotone
        is_monotone!(d, residue, tolzero) && break
    end

    C = zeros(n,N)
    C[:,1] = y - f[:,1]
    for i = 2:n_modes
        C[:,i] = f[:,i-1]-f[:,i]
    end

    return C
end

minorder(x) = length(x) > 1 ? length(x) - 1 : 1

"""
    function sift!(f, avg, residue, eps, use_fortran = false)
Interpolate a spline through the maxima and minima then get the average of the
two. Everything is done in place and nothing is returned.

If use_fortrain is true then we use Dierckx (a fortran library) for interpolation. Otherwise we use (a fork of_ interpolations.jl - which is very beta at the moment).
"""
function sift!(fi, avg, residue, data, t, eps, use_fortran=true)
    # Interpolate a spline through the maxima and minima
    tmax, tmin, max_ar, min_ar  = findExtrema(residue)

    debug = get_mean_envelope!(avg,t,tmax,max_ar,tmin,min_ar,use_fortran)

    residue .= residue .- avg

    sd = mean((avg.^2) ./ ((data-fi).^2 + eps) )

    fi .= fi .+ avg
    return nothing
end

function get_mean_envelope!(avg,t,tmax,max_ar,tmin,min_ar,use_fortran)
    maxorder = 4
	if use_fortran
		p_max = (length(max_ar) > maxorder) ? maxorder : minorder(max_ar)
		p_min = (length(min_ar) > maxorder) ? maxorder : minorder(min_ar)

		S1 = Spline1D(tmax, max_ar, k = p_max)
		S2 = Spline1D(tmin, min_ar, k = p_min)
		avg   .= (S1(t) .+ S2(t)) / 2
        return avg,S1(t),S2(t)
	else
        S1 = interpolate((tmax,), max_ar, Gridded(Cubic(Natural())))
        S2 = interpolate((tmin,), min_ar, Gridded(Cubic(Natural())))
        avg   .= (S1[t] .+ S2[t]) / 2
        return avg,S1[t],S2[t]
    end
end

function is_monotone!(d, residue, tolzero)
    c  = mean(abs.(residue))
    d .= diff(residue)
    is_m = all(d .+ c .* tolzero .> 0) || all(d .- c .* tolzero .< 0)
    return is_m
end

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition


text:
include("lib_for_robot.jl")
function task6(r)
    x,y=tostart(r)
    side=Ost
    bside=Nord
    side = moveandcheck(r,side)
    markerBorder(r,side,bside)
    bside=reverse(side)
    side=Nord
    markerBorder(r,side,bside)
    side=bside
    bside=Sud
    markerBorder(r,side,bside)
    bside=reverse(side)
    side=Sud
    markerBorder(r,side,bside)
    tostart(r)
    while x>0
        if isborder(r,Ost)
            x-=obojti(r,Ost)
        else
            move!(r,Ost)
            x-=1
        end
    end
    while y>0
        if isborder(r,Nord)
            y-=obojti(r,Nord)
        else
            move!(r,Nord)
            y-=1
        end
    end
end