text:
module TestConstraintFunctionize

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_scalar_functionize()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.ScalarFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    ci = MOI.add_constraint(bridged_mock, x, MOI.GreaterThan(0.0))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ x
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(0.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.GreaterThan(1.0))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(1.0)
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),),
    )
    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(),
        include = [
            "test_basic_VariableIndex_GreaterThan",
            "test_basic_VariableIndex_LessThan",
        ],
    )
    for T in [Int, Float64], S in [MOI.GreaterThan{T}, MOI.LessThan{T}]
        @test MOI.Bridges.added_constraint_types(
            MOI.Bridges.Constraint.ScalarFunctionizeBridge{T,S},
        ) == [(MOI.ScalarAffineFunction{T}, S)]
    end
    return
end

function test_scalar_functionize_linear2()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.ScalarFunctionize{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0, 1],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.NONBASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC, MOI.NONBASIC],
            ],
            variable_basis_status = [MOI.BASIC, MOI.NONBASIC_AT_LOWER],
        ),
    )
    MOI.Test.test_linear_integration_2(
        bridged_mock,
        MOI.Test.Config(
            exclude = Any[MOI.VariableBasisStatus, MOI.ConstraintBasisStatus],
        ),
    )
    cis = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VariableIndex,
            MOI.GreaterThan{Float64},
        }(),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        for ci in cis
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            MOI.set(bridged_mock, attr, ci, 2.0)
            @test MOI.get(bridged_mock, attr, ci) == 2.0
        end
    end
    ci = cis[1]
    attr = MOI.Test.UnknownConstraintAttribute()
    BT = MOI.Bridges.Constraint.ScalarFunctionizeBridge{
        Float64,
        MOI.GreaterThan{Float64},
    }
    function _unsupported_attribute(attr, f)
        return MOI.UnsupportedAttribute(
            attr,
            "Bridge of type `$(BT)` does not support $(f)ting the attribute " *
            "`$attr` because `MOIB.Constraint.invariant_under_function_conversion($attr)` returns `false`.",
        )
    end
    @test_throws(
        _unsupported_attribute(attr, "set"),
        MOI.set(bridged_mock, attr, ci, 1.0),
    )
    @test_throws(
        _unsupported_attribute(attr, "get"),
        MOI.get(bridged_mock, attr, ci),
    )
    for (i, ci) in enumerate(cis)
        _test_delete_bridge(
            bridged_mock,
            ci,
            2,
            ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),),
            num_bridged = 3 - i,
        )
    end
    return
end

function test_vector_functionize()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.VectorFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    z = MOI.add_variable(bridged_mock)
    v1 = MOI.VectorOfVariables([x, y, z])
    v2 = MOI.VectorOfVariables([x, y, y])
    ci = MOI.add_constraint(bridged_mock, v1, MOI.PowerCone(0.1))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v1
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, v2)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v2
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.1)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.PowerCone(0.2))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.2)
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((MOI.VectorAffineFunction{Float64}, MOI.PowerCone{Float64}, 0),),
    )

    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(),
        include = [
            "test_basic_VectorOfVariables_Nonnegatives",
            "test_basic_VectorOfVariables_Nonpositives",
        ],
    )
    for T in [Int, Float64], S in [MOI.Nonnegatives, MOI.Nonpositives]
        @test MOI.Bridges.added_constraint_types(
            MOI.Bridges.Constraint.VectorFunctionizeBridge{T,S},
        ) == [(MOI.VectorAffineFunction{T}, S)]
    end
    return
end

function test_vector_functionize_lin1v()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.VectorFunctionize{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[0, 2, 0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-3, -1]],
        )
    MOI.Test.test_conic_linear_VectorOfVariables(
        bridged_mock,
        MOI.Test.Config(),
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
    )
    func = MOI.get(bridged_mock, MOI.ConstraintFunction(), ci)
    MOI.delete(bridged_mock, func.variables[2])
    new_func = MOI.VectorOfVariables(func.variables[[1, 3]])
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) == new_func
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.Nonnegatives(2)
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, [1.0, 2.0])
        @test MOI.get(bridged_mock, attr, ci) == [1.0, 2.0]
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),),
    )
    return
end

function test_unsupported_attribute_issue_1758()
    inner = MOI.Utilities.Model{Float64}()
    model = MOI.Bridges.Constraint.ScalarFunctionize{Float64}(inner)
    MOI.Test.test_model_copy_to_UnsupportedAttribute(model, MOI.Test.Config())
    return
end

end  # module

TestConstraintFunctionize.runtests()


text:
using Flux
println("NNPDE_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralNetDiffEq

## Example 1, 1D ode
@parameters t θ
@variables u(..)
@derivatives Dt'~t

# 1D ODE
eq = Dt(u(t,θ)) ~ t^3 + 2*t + (t^2)*((1+3*(t^2))/(1+t+(t^3))) - u(t,θ)*(t + ((1+3*(t^2))/(1+t+t^3)))

# Boundary conditions
bcs = [u(0.) ~ 1.0 , u(1.) ~ 1.202]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.1
discretization = PhysicsInformedNN(dx)

# Neural network and optimizer
opt = Flux.ADAM(0.1)
chain = FastChain(FastDense(1,12,Flux.σ),FastDense(12,1))

pde_system = PDESystem(eq,bcs,domains,[t],[u])
prob = discretize(pde_system,discretization)
alg = NNDE(chain,opt,autodiff=false)
phi,res = solve(prob,alg,verbose=true, maxiters=1000)

analytic_sol_func(t) = exp(-(t^2)/2)/(1+t+t^3) + t^2
ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(phi(t,res.minimizer)) for t in ts]

@test u_predict ≈ u_real atol = 0.1

# t_plot = collect(ts)
# plot(t_plot ,u_real)
# plot!(t_plot ,u_predict)

## Example 2, 2D Poisson equation
@parameters x y θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# 2D PDE
eq  = Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1
discretization = PhysicsInformedNN(dx)

# Neural network and optimizer
opt = Flux.ADAM(0.1)
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = discretize(pde_system,discretization)
alg = NNDE(chain,opt,autodiff=false)
phi,res  = solve(prob,alg,verbose=true, maxiters=600)

xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))

@test u_predict ≈ u_real atol = 0.3

# p1 =plot(xs, ys, u_predict, st=:surface);
# p2 = plot(xs, ys, u_real, st=:surface);
# plot(p1,p2)

## Example , 3D PDE
@parameters x y t θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y
@derivatives Dtt''~t

# 3D PDE
eq  = Dt(u(x,y,t,θ)) ~ Dxx(u(x,y,t,θ)) + Dyy(u(x,y,t,θ))
# Boundary conditions
bcs = [u(x,y,0) ~ exp(x+y)*cos(x+y) ,
       u(x,y,2) ~ exp(x+y)*cos(x+y+4*2) ,
       u(0,y,t) ~ exp(y)*cos(y+4t),
       u(2,y,t) ~ exp(2+y)*cos(2+y+4t) ,
       u(x,0,t) ~ exp(x)*cos(x+4t),
       u(x,2,t) ~ exp(x+2)*cos(x+2+4t)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]

# Discretization
dx = 0.25
discretization = PhysicsInformedNN(dx)

# Neural network and optimizer
opt = Flux.ADAM(0.08)
# opt = BFGS()
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

pde_system = PDESystem(eq,bcs,domains,[x,y,t],[u])
prob = discretize(pde_system,discretization)
alg = NNDE(chain,opt,autodiff=false)
phi,res  = solve(prob,alg,verbose=true, maxiters=2100)

xs,ys,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
analytic_sol_func(x,y,t) = exp(x+y)*cos(x+y+4t)
u_real = [reshape([analytic_sol_func(x,y,t) for x in xs  for y in ys], (length(xs),length(ys)))  for t in ts ]
u_predict = [reshape([first(phi([x,y,t],res.minimizer)) for x in xs  for y in ys], (length(xs),length(ys)))  for t in ts ]
@test u_predict ≈ u_real atol = 80.0

# p1 =plot(xs, ys, u_predict, st=:surface);
# p2 = plot(xs, ys, u_real, st=:surface);
# plot(p1,p2)

ext: jl


text:
####################
# PullRequest Type #
####################

mutable struct PullRequest <: GitHubType
    base::Union{Branch, Nothing}
    head::Union{Branch, Nothing}
    number::Union{Int, Nothing}
    id::Union{Int, Nothing}
    comments::Union{Int, Nothing}
    commits::Union{Int, Nothing}
    additions::Union{Int, Nothing}
    deletions::Union{Int, Nothing}
    changed_files::Union{Int, Nothing}
    state::Union{String, Nothing}
    title::Union{String, Nothing}
    body::Union{String, Nothing}
    merge_commit_sha::Union{String, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    closed_at::Union{Dates.DateTime, Nothing}
    merged_at::Union{Dates.DateTime, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    assignee::Union{Owner, Nothing}
    user::Union{Owner, Nothing}
    merged_by::Union{Owner, Nothing}
    milestone::Union{Dict, Nothing}
    _links::Union{Dict, Nothing}
    mergeable::Union{Bool, Nothing}
    merged::Union{Bool, Nothing}
    locked::Union{Bool, Nothing}
end

PullRequest(data::Dict) = json2github(PullRequest, data)
PullRequest(number::Real) = PullRequest(Dict("number" => number))

namefield(pr::PullRequest) = pr.number

###############
# API Methods #
###############

@api_default function pull_requests(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/pulls"; options...)
    return map(PullRequest, results), page_data
end

@api_default function pull_request(api::GitHubAPI, repo, pr; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/pulls/$(name(pr))"; options...)
    return PullRequest(result)
end

@api_default function update_pull_request(api::GitHubAPI, repo, pr; options...)
    result = gh_patch_json(api, "/repos/$(name(repo))/pulls/$(name(pr))"; options...)
    return PullRequest(result)
end

@api_default function close_pull_request(api::GitHubAPI, repo, pr; options...)
    update_pull_request(api, repo, pr; params = Dict(
        :state => "closed"
    ), options...)
end

@api_default function create_pull_request(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/pulls"; options...)
    return PullRequest(result)
end

@api_default function create_comment(api::GitHubAPI, repo, pr::PullRequest, body::AbstractString; options...)
    create_comment(api, repo, pr, :pr; params = Dict(
        :body => body
    ), options...)
end


text:
export Combiner

struct CombinerStorage <: TensorStorage
  ci::Index
end

function storage_contract(CSstore::CombinerStorage,
                          Cis::IndexSet,
                          Dstore::Dense,
                          dis::IndexSet)
  cind = Cis[1]
  if hasindex(dis, cind) # has combined index, uncombine
    cpos = findindex(dis,cind)
    dinds = inds(dis)
    Cinds = inds(Cis)
    Nis = IndexSet(vcat(dinds[1:cpos-1], Cinds[2:end], dinds[cpos+1:end]))
    return Nis, Dstore
  else # lacks combined index, combine
    # dis doesn't have cind, replace
    # Cis[1], Cis[2], ... with cind, may need to permute
    j_1_pos = findindex(dis,Cis[2])
    if j_1_pos < 1
      throw(ArgumentError("tensor missing index $(Cis[2]) in combiner-tensor product"))
    end

    # Check if Cis[2], Cis[3], ... are grouped together (contiguous)
    # and in same order as on combiner
    contig_sameord = true
    j = j_1_pos+1
    for c=3:length(Cis)
      if j > length(dis) || (dis[j] != Cis[c])
        contig_sameord = false
        break
      end
      j += 1
    end

    if contig_sameord
      lend = j_1_pos-1
      rstart = j_1_pos+length(Cis)-1
      dinds = inds(dis)
      Nis = IndexSet(vcat(dinds[1:lend], cind, dinds[rstart:end]))
    else # permutation required
      # Set P destination values to -1 to mark indices that need to be assigned destinations:
      P = fill(-1, length(dis))
      # permute combined indices to the front, in same order as in Cis:
      ni = 1
      for c in 2:length(Cis)
        j = findindex(dis, Cis[c])
        if j < 1
          throw(ArgumentError("tensor missing index $(Cis[c]) in combiner-tensor product"))
        end
        P[j] = ni
        ni += 1
      end

      Nis = IndexSet(length(dis)+2-length(Cis))
      Nis[1] = cind
      i = 2
      for j in 1:length(dis)
        if P[j] == -1
          P[j] = ni
          ni += 1
          Nis[i] = dis[j]
          i += 1
        end
      end

      ddata = vec(permutedims(reshape(data(Dstore), dims(dis)), invperm(P)))
      Dstore = Dense{eltype(Dstore)}(ddata)
    end
    return Nis, Dstore
  end
end

function storage_contract(Dstore::Dense,
                          dis::IndexSet,
                          CSstore::CombinerStorage,
                          Cis::IndexSet)
  return storage_contract(CSstore, Cis, Dstore, dis)
end


text:
"""
Colorscheme: phase plot with conformal polar grid
"""
function cs_c(; colormap=hsv_colors(), phaseres=20)
    stepfct = generate_stepfct(length(colormap))
    sawfct1 = generate_sawfct(1.0/phaseres, 0.75, 1.0)
    sawfct2 = generate_sawfct(2*π/phaseres, 0.75, 1.0)

    let cm = RGB.(colormap)
      return function(z, fz)
          anglefz = angle_zero_one(fz)
          absfz = abs(fz)
          factor = sawfct1(anglefz) * (
            absfz == 0.0 ? 0.75 : absfz == Inf ? 1.0 : sawfct2(log(absfz)))
          color = cm[stepfct(anglefz)]
          return brightenRGB(
            RGB(factor*color.r, factor*color.g, factor*color.b), 0.1)
      end
    end
end

"""
Colorscheme: colored phase plot with module jumps
"""
function cs_m(; colormap=hsv_colors(), logabsres=20)
    sawfct = generate_sawfct(2*π/logabsres, 0.7, 1.0)
    stepfct = generate_stepfct(length(colormap))
    let cm = RGB.(colormap)
      return function(z, fz)
        absfz = abs(fz)
        factor = absfz == 0.0 ? 0.7 :
                 absfz == Inf ? 1.0 : sawfct(log(absfz))
        color = cm[stepfct(angle_zero_one(fz))]
        return RGB(factor*color.r, factor*color.g, factor*color.b)
      end
    end
end

"""
Colorscheme: enhanced domain coloring
"""
cs_e(; colormap=hsv_colors(), logabsres=20) = cs_mod_abs_brightness(
    cs_m(colormap=colormap, logabsres=logabsres))

"""
Colorscheme: colored phase plot with specific phase jumps
"""
function cs_j(; colormap=hsv_colors(),
                jumps=collect(LinRange(-π, π, 21))[1:20], a=0.8, b=1.0)
    sawfct = generate_sawfctx(
      angle_zero_one.(exp.(1im .* jumps)); a=a, b=b)
    stepfct = generate_stepfct(length(colormap))
    let cm = RGB.(colormap)
      return function(z, fz)
          angle_fz = angle_zero_one(fz)
          factor = sawfct(angle_fz)
          color = cm[stepfct(angle_fz)]
          return RGB(factor*color.r, factor*color.g, factor*color.b)
      end
    end
end

"""
Colorscheme: proper phase plot
"""
function cs_p(; colormap=hsv_colors())
    stepfct = generate_stepfct(length(colormap))
    let cm = copy(colormap)
      return (z, fz) -> cm[stepfct(angle_zero_one(fz))]
    end
end

"""
Colorscheme: phase plot colored in steps
"""
cs_q(; colormap=hsv_colors(20)) = cs_p(; colormap=colormap)

"""
Colorscheme: standard domain coloring
"""
cs_d(; colormap=hsv_colors()) = cs_mod_abs_brightness(
  cs_p(colormap=colormap))

"""
Colorscheme: for exams; cs_j with NIST colors and more pronounced shadow jumps.
"""
cs_exam(; colormap=nist_colors(), jumps=collect(LinRange(-π, π, 21))[1:20],
          a=0.6, b=1.0) = cs_j(colormap=colormap, jumps=jumps, a=a, b=b)

"""
Colorscheme: Stripes in (real(z), imag(z)) plane

Directions is a vector of complex numbers.
For each such number d in directions the stripes are perpendicular
to (real(d), imag(d)). And abs(d) is the number of stripes per unit
length.
This devides the complex plane in two sets Sd and ℂ\\Sd.

Each z ∊ ℂ has n=length(directions) flags (z ∊ S1, z ∊ S1, ..., z ∊ Sn).
This make 2^n possible flag-combinations. For every combination there
need to be a color in colors.
"""
function cs_stripes(directions, colors)
  let no_dirs = length(directions)
    if 2^no_dirs != length(colors)
      error("for $(no_dirs) directions $(2^no_dirs) colors are needed.")
    end
    rectfct = generate_rectfct(1.0)
    return function(z, fz)
      color_index = 0
      color_pos = 1
      for dir in directions
        color_index += color_pos*Int(rectfct(real(fz*dir)))
        color_pos *= 2
      end
      return colors[1 + color_index]
    end
  end
end

"""
Colorscheme: Stripes in (angle(z), abs(z)) plane

Use `cs_stripes` for (angle(z), abs(z)).
"""
function cs_angle_abs_stripes(directions, colors)
  col_fct = cs_stripes(directions, colors)
  return (z, fz) -> col_fct(z,angle(fz) + 1im*abs(fz))
end

"""
Colorscheme: uniform low-brightness grid lines in real- and imag-part and
             low saturation grid lines in lowsat(abs(fz)).
"""
function cs_gridReIm_lowsat(; reim_a=-9.0, reim_b=0.5, abs_a=-9.0, abs_b=0.8,
                              lowsat=absfz -> log(exp(0.5)*absfz),
                              vcorr=(s,v) -> max(1.0 - s, v))
  sp_reim = generate_spike(a=reim_a, b=reim_b)
  sp_abs = generate_spike(a=abs_a, b=abs_b)
  return function(z, fz)
    absfz = abs(fz)
    h = nist_trafo(4*angle_zero_one(fz))
    s = absfz == 0.0 ? 1.0 : absfz == Inf ? 0.0 : sp_abs(lowsat(absfz))
    v = absfz == 0.0 ? 0.0 :
        absfz == Inf ? 1.0 : min(sp_reim(real(fz)-0.5), sp_reim(imag(fz)-0.5))
    return HSV(h, s, vcorr(s,v))
  end
end

"""
Colorscheme: uniform low-brightness grid lines in real- and imag-part and
             low saturation grid lines in log(abs(fz)).
"""
cs_gridReIm_logabs(; reim_a=-9.0, reim_b=0.5, abs_a=-9.0, abs_b=0.8,
                     vcorr=(s,v) -> max(1.0 - s, v)) = cs_gridReIm_lowsat(
  reim_a=reim_a, reim_b=reim_b, abs_a=abs_a, abs_b=abs_b, vcorr=vcorr,
  lowsat=absfz -> log(exp(0.5)*absfz))

"""
Colorscheme: uniform low-brightness grid lines in real- and imag-part and
             low saturation grid lines in abs(fz).
"""
cs_gridReIm_abs(; reim_a=-9.0, reim_b=0.5, abs_a=-9.0, abs_b=0.8,
                  vcorr=(s,v) -> max(1.0 - s, v)) = cs_gridReIm_lowsat(
  reim_a=reim_a, reim_b=reim_b, abs_a=abs_a, abs_b=abs_b, vcorr=vcorr,
  lowsat=absfz -> absfz+0.5)

"""
Colorscheme: uniform low-brightness grid lines in real- and imag-part.
"""
cs_gridReIm(; reim_a=-9.0, reim_b=0.5) = cs_gridReIm_abs(
  reim_a=reim_a, reim_b=reim_b, abs_a=-Inf, vcorr=(s,v) -> v)

"""
Colorscheme: use the pixels of an image as color source.
"""
function cs_useImage(image;
                      img_upperleft=0.0 + 1.0im, img_lowerright=1.0 + 0.0im,
                      is_hv_periodic=(true, true))
  let (img_height, img_width) = size(image),
      x_left = real(img_upperleft),
      x_width = real(img_lowerright) - x_left,
      y_top = imag(img_upperleft),
      y_height = y_top - imag(img_lowerright),
      epsilon = eps()

      width_fct = generate_stepfct(img_width)
      height_fct = generate_stepfct(img_height)

      return function (z,fz)
        x = (real(fz) - x_left)/x_width
        y = (y_top - imag(fz))/y_height
        if !is_hv_periodic[1]
          x = min(max(0.0, x), 1.0-epsilon)
        end
        if !is_hv_periodic[2]
          y = min(max(0.0, y), 1.0-epsilon)
        end
        return image[height_fct(y), width_fct(x)]
      end
  end
end

# vim:syn=julia:cc=79:fdm=indent:sw=2:ts=2:


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Empty <: ProtoType
    Empty(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Empty

export Empty


text:
p = 2
n = 10000

function ftrefethen(p)
  x = p[1]; y = p[2]
  return exp(sin(50 * x)) + sin(60 * exp(y)) + sin(70 * sin(x)) +
             sin(sin(80 * y)) - sin(10 * (x + y)) + (x^2 + y^2)/4
end
f = ftrefethen

alpha = 0.602
gamma = 0.101

theta_lo = -100*ones(p,1)
theta_hi = 100*ones(p,1)
theta = 0.40 * (theta_hi - theta_lo) * rand(p, 1)
theta=min(theta,theta_hi)
theta=max(theta,theta_lo)

# Input to gain settings calculations
initial_step_length = 0.40 * maximum(theta_hi - theta_lo)
gavg = 2 # How many gradient estimates per iteration
n_for_gain_calc = n/20
estimated_std_dev_measurement_noise = 0.0

A = 0.10 * n / (2*gavg)
c = max(estimated_std_dev_measurement_noise/(gavg^0.5), 0.0001)

gbar=zeros(p, 1)
for i=1:int(n_for_gain_calc/(2*gavg))
   ghat=zeros(p,1)
   for j=1:gavg
      delta=2*round(rand(p,1))-1;
      thetaplus=theta+c*delta;
      thetaminus=theta-c*delta;
      yplus=f(thetaplus);
      yminus=f(thetaminus);
      ghat=(yplus-yminus)./(2*c*delta)+ghat;
   end
   gbar=gbar+abs(ghat/gavg);
end
meangbar=mean(gbar)/(n_for_gain_calc/(2*gavg))
a = initial_step_length * ((A+1)^alpha)/meangbar;

println("a = $(a), c = $(c), A = $(A)")

#a = 0.0017
#c = 0.00001 # No noise so we can select c small
#A = 5

min_f = Inf
best = 1

for k=0:n-1

    ak = a/(k+1+A)^alpha
    ck = c/(k+1)^gamma

    delta = 2*round(rand(p,1))-1

    thetaplus = theta + ck*delta
    thetaminus = theta - ck*delta

    yplus=f(thetaplus)
    yminus=f(thetaminus)

    ghat = (yplus - yminus)./(2*ck*delta)
    theta = theta - ak * ghat

    theta=min(theta,theta_hi)
    theta=max(theta,theta_lo)

    tf = f(theta)

    if tf < min_f
      min_f = tf
      best = theta
    end

    println("theta = ", theta, ", fitness = ", )

end
println("a = $(a), c = $(c), A = $(A)")
println("best = ", best, ", fitness = ", min_f)


text:
export  teleportationBasedSympGenerator

# Calculate the effective Symplectic Matrix output of the adaptive Gaussian control scheme
function teleportationBasedSympGenerator(S::Symp, in_modes::Vector, out_modes::Vector)::Symp
    all_modes = [1:nModes(S);]
    in_ = vcat( [ [2*i - 1, 2 * i] for i in in_modes ]... )
    out = vcat( [ [2*i - 1, 2 * i] for i in out_modes ]... )
    ancilliary_modes = [i for i in all_modes if !(i in in_modes)]
    idle_modes = [i for i in all_modes if !(i in out_modes)]
    usq = 2*ancilliary_modes
    hm = 2*idle_modes
    S_out_in = S.S[ out, in_  ]
    S_hm_in = S.S[ hm, in_ ]
    S_out_usq = S.S[ out, usq ]
    S_hm_usq = S.S[ hm, usq ]
    return S_out_in - S_out_usq * S_hm_usq^-1 * S_hm_in
end


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/bkamins/ABCDGraphGenerator.jl"))


text:
export Integer, Node, InputNode, InteriorNode, OutputNode, print_node_tuple, gate_tuple

abstract type Node end
mutable struct Func
    func::Function
    arity::Integer
    name::AbstractString
end
Func(f::Function, a::Integer) = Func(f, a, string(f))

mutable struct InputNode <: Node
    index::Integer
    active::Bool
    cache::MyInt
end

function InputNode(index::Integer)
    return InputNode(index, false, 0 )
end

mutable struct InteriorNode <: Node
    func::Func
    inputs::Vector{Integer}
    active::Bool
    cache::MyInt
end

function InteriorNode(func::Func, inputs::Vector{Int64})
    return InteriorNode(func, inputs, false, MyInt(0) )
end

function InteriorNode(func::Func, inputs::Vector{Integer})
    return InteriorNode(func, inputs, false, MyInt(0) )
end

mutable struct OutputNode <: Node
    input::Integer
end

#
function print_node_tuple( f::IO, node_vect::Vector{InputNode} )
  print(f,"(")
  len = length(node_vect)
  for i in 1:(len-1)
    print(f,node_vect[i].index,",")
  end
  if len > 1
    print(f,node_vect[len].index,")")
  else
    print(f,node_vect[len].index,",)") # add comma so 1-element tuple will be read as a tuple
  end
end

function print_node_tuple( node_vect::Vector{InputNode} )
  print_node_tuple( Base.stdout, node_vect )
end

function print_node_tuple( f::IO, node_vect::Vector{InteriorNode} )
  print(f," (")
  len = length(node_vect)
  for i in 1:(len-1)
    print(f,"(",node_vect[i].func.name,",",node_vect[i].inputs,"),")
  end
  if len > 1
    print(f,"(",node_vect[len].func.name,",",node_vect[len].inputs,"))")
  else
    # add comma so 1-element tuple will be read as a tuple
    print(f,"(",node_vect[len].func.name,",",node_vect[len].inputs,"),)")
  end
end

function print_node_tuple( node_vect::Vector{InteriorNode} )
  print_node_tuple( Base.stdout, node_vect )
end

function print_node_tuple( f::IO, node_vect::Vector{OutputNode} )
  print(f," (")
  len = length(node_vect)
  for i in 1:(len-1)
    print(f,node_vect[i].input,",")
  end
  if len > 1
    print(f,node_vect[len].input,")")
  else
    print(f,node_vect[len].input,",)")
  end
end

function print_node_tuple( node_vect::Vector{OutputNode} )
  print_node_tuple( Base.stdout, node_vect )
end

# Outputs an vector of InteriorNodes to IO stream f in the format needed by circuit() and print_circuit() in Chromosome.jl
function gate_tuple( f::IO, node_vect::Vector{InteriorNode}, numinputs::Int64 )
  print(f," (")
  len = length(node_vect)
  for i in 1:(len-1)
    print(f,"(",numinputs+i,",",node_vect[i].func.name,",",node_vect[i].inputs[1],",",node_vect[i].inputs[2],"), ")
  end
  if len > 1
    print(f,"(",numinputs+len,",",node_vect[len].func.name,",",node_vect[len].inputs[1],",",node_vect[len].inputs[2],"))")
  else
    # add comma so 1-element tuple will be read as a tuple
    print(f,"(",numinputs+len,",",node_vect[len].func.name,",",node_vect[len].inputs[1],",",node_vect[len].inputs[2],"),)")
  end
end

#= moved to Chromosome.jl
function print_build_chromosome( f::IO, c::Chromosome )
  println(f, "build_chromosome(")
  print_node_tuple(f, c.inputs )
  print(f,",")
  print_node_tuple(f, c.interiors )
  print(f,",")
  print_node_tuple(f, c.outputs )
  println(f, ")")
end
=#


text:
CONFIG = Dict()
CONFIG["wordvec_file"] = ".data/word2vec.NYT.100d.h5"
CONFIG["train_file"] = ".data/CoNLL2009-ST-English-train.txt"
CONFIG["test_file"] = ".data/CoNLL2009-ST-evaluation-English.txt"
CONFIG["nepochs"] = 50
CONFIG["learning_rate"] = 0.01
CONFIG["batchsize"] = 25
CONFIG["training"] = true


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module RiskAdjustedLinearizations

import Base.show
using ForwardDiff, LinearAlgebra, OrderedCollections, Printf, SparseArrays# , SparseDiffTools, SparsityDetection, UnPack
using UnPack
using DiffEqBase: DiffCache, get_tmp, dualcache

# Utilities
include("util.jl")

# RiskAdjustedLinearization
include("risk_adjusted_linearization.jl")

# Solution Algorithms
# include("solve.jl")
# include("qzdecomp.jl")
include("solution_algorithms/blanchard_kahn.jl")

export RiskAdjustedLinearization, update!, nonlinear_system, linearized_system
end # module


text:
using InvariantMeasures
using ValidatedNumerics

using Plots
using LaTeXStrings

D = Mod1Dynamic(x -> 4*x + 0.01*InvariantMeasures.sinpi(8*x))
#D = Mod1Dynamic(x -> 16*x + 0.01*InvariantMeasures.sinpi(32*x))
#D = Mod1Dynamic(x->2*x+0.5*x*(1-x))
# D = PwMap(
# [x -> 17x/5,
#  x -> 34(x-5//17)^2/25 + 3(x-5//17),
#  x -> 34(x-10//17)^2/25 + 3(x-10//17),
#  x -> 17(x-15//17)/5
# ],
# [0, @interval(5/17), @interval(10/17), @interval(15/17), 1]
# )

num_norms = 30

B = Hat(1024)
Q = DiscretizedOperator(B, D)

trivial_norms = norms_of_powers_trivial(weak_norm(B), Q, num_norms)
computed_norms = norms_of_powers(weak_norm(B), num_norms, Q, integral_covector(B))

(dfly_strongs, dfly_norms) = norms_of_powers_dfly(B, D, num_norms)

norms = min.(trivial_norms, computed_norms, dfly_norms)
better_norms = refine_norms_of_powers(norms, num_norms)

pgfplotsx()
p = plot(trivial_norms,
    label = L"$\|Q\|^k$",
    yscale = :log10,
    legend = :bottomleft,
    title = "Various bounds for norms of powers",
    xlabel = L"$k$",
    ylabel = L"bound to $\|Q^k|_U\|$"
    )
plot!(p, computed_norms,
    label = "computational bounds")

plot!(p, dfly_norms,
    label = L"DFLY $2\times 2$ matrix bounds")

plot!(p, better_norms,
    label = "min(previous) + refinement")

# savefig(p, "norm_bounds_kinds.tikz")


text:
@reexport module ArgumentParsers

using TSML
using TSML.TSMLTypes
using TSML.Utils

using TSML.ValDateFilters: CSVDateValReader, CSVDateValWriter
using TSML.Statifiers
using TSML.Monotonicers
using TSML.Outliernicers
using TSML.Plotters

using Dates
using DataFrames
using CSV
using ArgParse

export tsmlmain

const COMMONARG = Dict(
        :dateformat=>"dd/mm/yyyy HH:MM",
        :dateinterval=>Dates.Hour(1),
        :processmissing=>true
      )

const DATEINTERVAL = Dict(
        "minutely" => Dates.Minute(1),
        "hourly" => Dates.Hour(1),
        "daily"=> Dates.Day(1),
        "weekly" => Dates.Week(1),
        "monthly"=> Dates.Month(1),
        "yearly"=> Dates.Year(1)
       )

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--aggregate"
            help = "aggregate interval such as: minutely, hourly, weekly, monthly, Dates.Minute(30),Dates.Hour(2)"
            arg_type = String
            default = "hourly"
        "--dateformat"
            help = "date and time format"
            arg_type = String
            default = "dd/mm/yyyy HH:MM"
        "--impute"
            help = "impute with NN"
            action = :store_true
        "--mono"
            help = "normalize monotonic"
            action = :store_true
        "--clean"
            help = "remove outliers"
            action = :store_true
        "--stat"
            help = "get statistics"
            action = :store_true
        "--save"
            help = "save results"
            action = :store_true
        "csvfilename"
            help = "csv file to process"
            required = true
    end

    return parse_args(s,as_symbols=true)
end

function tsmlmain()
    parsed_args = parse_commandline()
    createrunpipeline(parsed_args)
end

function createrunpipeline(args::Dict)
    dateformat = args[:dateformat]
    fname = args[:csvfilename]
    csvreader = CSVDateValReader(Dict(:filename=>fname,:dateformat=>dateformat))

    csvwriter = CSVDateValWriter()
    fnameoutput = ""
    if args[:save] == true
        name = split(fname,".csv")[1]
        fnameoutput = name*"_output"*".csv"
        csvwriter.args[:filename] = fnameoutput
    end

    # extract date interval
    dateintervalstring = args[:aggregate]
    dateinterval = Dates.Hour(1)

    if dateintervalstring in keys(DATEINTERVAL)
       dateinterval  = DATEINTERVAL[dateintervalstring]
    else
        # parse the string into expression and evaluate
        # example: "Dates.Hour(1)" or "Dates.Minute(30)"
        dateinterval = eval(Meta.parse(dateintervalstring))
    end
    commonarg = Dict(:dateformat=>dateformat,:dateinterval=>dateinterval)
    valgator = DateValgator(commonarg)
    imputer = DateValNNer(commonarg)
    mono = Monotonicer()
    outliernicer = Outliernicer(commonarg)
    statifier = Statifier(Dict(:processmissing=>true))

    transformers = [csvreader,valgator]
    if args[:impute] == true
        transformers = [transformers;imputer]
    end
    if args[:mono] == true
        transformers = [transformers;mono]
    end
    if args[:clean] == true
        transformers = [transformers;outliernicer]
    end
    if args[:stat] == true
        transformers = [transformers;statifier]
    end
    if args[:save] == true
        transformers = [transformers;csvwriter]
    end

    pipeline = Pipeline()
    pipeline.args[:transformers]=transformers
    fit!(pipeline)
    transform!(pipeline)
end

function tsmlrun(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM",otype::AbstractString="table")
    if otype == "table"
        imputedoutput(inputname,outputname,datefmt)
    elseif otype == "stat"
        imputedstat(inputname,outputname,datefmt)
    else
        error("wrong output type")
    end
end

function rawstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = CSVDateValReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function aggregatedstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = CSVDateValReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function aggregatedoutput(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = CSVDateValReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function imputedstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = CSVDateValReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    valnner = DateValNNer(COMMONARG)
    mononicer = Monotonicer()
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,valnner,mononicer,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function imputedoutput(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = CSVDateValReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    valnner = DateValNNer(COMMONARG)
    mononicer = Monotonicer()
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,valnner,mononicer]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

end


text:
#= Worker functions=#
low(i::Int, mid::Int, bpc::Int) = Int(mid + (i-1)*(1-mid)/bpc)
high(i::Int, mid::Int, bpc::Int, slen::Int, L::Int) = Int(slen - (slen-mid)*(i-L+bpc)/(bpc))

# used in general dirichlet BC. To simulate a constant value beyond the boundary
limit(i, N) = N>=i>=1 ? i : (i<1 ? 1 : N)

# used in Neumann 0 BC
function reflect(idx, L)
    abs1 = abs(L-idx)
    if L - abs1 > 0
        return L-abs1
    else
        return abs(L-abs1)+2
    end
end

# gives the index for periodic BC
function rem1(idx,L)
    r = idx%L
    if r > 0
        return r
    else
        return r+L
    end
end

#= LEFT BOUNDARY CONDITIONS =#
function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,:Dirichlet0,RBC}) where {T<:Real,S<:SVector,RBC}
    mid = div(A.stencil_length,2) + 1
    bpc = A.stencil_length - mid
    for i in 1 : A.boundary_point_count[1]
        dirichlet_0!(x_temp, x, A.stencil_coefs, mid, bpc, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,:Dirichlet0,RBC}) where {T<:Real,S<:SVector,RBC}
    mid = div(A.stencil_length,2) + 1
    bpc = A.stencil_length - mid
    for i in 1 : A.boundary_point_count[1]
        dirichlet_0!(x_temp, x, A.stencil_coefs[1], mid, bpc, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,:Dirichlet0,RBC}) where {T<:Real,S<:SVector,RBC}
    stencil_length = A.stencil_length
    stencil_rem = 1-stencil_length%2
    for i in 1 : A.boundary_point_count[1]
        A.directions[][i] ? start_idx = stencil_length - stencil_rem : start_idx = 1 + stencil_rem
        # we have to modify the number of boundary points to be considered as with upwind operators
        # the number of bpc is only 0 or 1 depending on the order
        A.directions[][i] ? bpc = A.boundary_point_count[1] : bpc = stencil_rem
        # println("***D0 i = $i, start_idx/mid = $start_idx, bpc = $bpc, stencil_length = $stencil_length ***")
        dirichlet_1!(x_temp, x, A.directions[][i] ? A.down_stencil_coefs : A.up_stencil_coefs, start_idx, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,:Dirichlet,RBC}) where {T<:Real,S<:SVector,RBC}
    x[1] = A.boundary_condition[][1][3](A.t)
    mid = div(A.stencil_length,2)+1
    for i in 1 : A.boundary_point_count[1]
        dirichlet_1!(x_temp, x, A.stencil_coefs, mid, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,:Dirichlet,RBC}) where {T<:Real,S<:SVector,RBC}
    x[1] = A.boundary_condition[][1][3](A.t)
    mid = div(A.stencil_length,2)+1
    for i in 1 : A.boundary_point_count[1]
        dirichlet_1!(x_temp, x, A.stencil_coefs[1], mid, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,:periodic,RBC}) where {T<:Real,S<:SVector,RBC}
    for i in 1 : A.boundary_point_count[1]
        periodic!(x_temp, x, A.stencil_coefs, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,:periodic,RBC}) where {T<:Real,S<:SVector,RBC}
    for i in 1 : A.boundary_point_count[1]
        periodic!(x_temp, x, A.stencil_coefs[1], i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,:Neumann0,RBC}) where {T<:Real,S<:SVector,RBC}
    for i in 1 : A.boundary_point_count[1]
        neumann0!(x_temp, x, A.stencil_coefs, i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,:Neumann0,RBC}) where {T<:Real,S<:SVector,RBC}
    for i in 1 : A.boundary_point_count[1]
        neumann0!(x_temp, x, A.stencil_coefs[1], i)
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,:Neumann,RBC},FiniteDifference{T,S,:Neumann,RBC}}) where {T<:Real,S<:SVector,RBC}
    @inbounds for i in 1 : A.boundary_point_count[1]
        bc = A.low_boundary_coefs[][i]
        tmp = zero(T)
        @inbounds for j in 1 : length(bc)
            tmp += bc[j] * x[j]
        end
        x_temp[i] = tmp
    end
    x_temp[1] += A.boundary_condition[][1][3]
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,:Robin,RBC},FiniteDifference{T,S,:Robin,RBC}}) where {T<:Real,S<:SVector,RBC}
    @inbounds for i in 1 : A.boundary_point_count[1]
        bc = A.low_boundary_coefs[][i]
        tmp = zero(T)
        @inbounds for j in 1 : length(bc)
            tmp += bc[j] * x[j]
        end
        x_temp[i] = tmp
    end
    x_temp[1] += A.boundary_condition[][1][3]
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,:None,RBC},UpwindOperator{T,S,:None,RBC},FiniteDifference{T,S,:None,RBC}}) where {T<:Real,S<:SVector,RBC}
    halfstencil = div(A.stencil_length, 2)
    for i in 1 : A.boundary_point_count[1]
        @inbounds bc = A.low_boundary_coefs[][i]
        tmp = zero(T)
        startid = max(0,i-1-halfstencil)
        @inbounds for j in 1 : length(bc)
            tmp += bc[j] * x[startid+j]
        end
        @inbounds x_temp[i] = tmp
    end
end

function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,:nothing,RBC}) where {T<:Real,S<:SVector,RBC}
    stencil_length = length(A.down_stencil_coefs)
    stencil_rem = 1 - stencil_length%2
    start_idx = 1

    # this case when our stencil uses 1 point at the left, so we cant use the upwind stencil
    # as it will spill over. So we use a special boundary stencil.

    if stencil_rem == 1
        x_temp[1] = sum(A.low_boundary_coefs[][1].*x[1:length(A.low_boundary_coefs[][1])])
        start_idx = 2
    end

    for i in start_idx : A.boundary_point_count[1]
        xtempi = zero(T)
        # startid = max(0,i-1-halfstencil)
        @inbounds for j in 1:length(A.up_stencil_coefs)
            xtempi += A.up_stencil_coefs[j] * x[i+j-1-stencil_rem]
        end
        @inbounds x_temp[i] = xtempi
    end
end

#= INTERIOR CONVOLUTION =#
function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC}
    N = length(x)
    coeffs = A.stencil_coefs
    mid = div(A.stencil_length, 2) + 1

    Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2]
        # dirichlet_0!(x_temp,x,A.stencil_coefs, i)
        xtempi = zero(T)
        @inbounds for idx in 1:A.stencil_length
            xtempi += coeffs[idx] * x[i - (mid-idx)]
        end
        x_temp[i] = xtempi
    end
end

function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC}
    N = length(x)
    coeffs = A.stencil_coefs
    mid = div(A.stencil_length, 2) + 1

    Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2]
        # dirichlet_0!(x_temp,x,A.stencil_coefs, i)
        xtempi = zero(T)
        @inbounds for idx in 1:A.stencil_length
            xtempi += coeffs[i-mid+1][idx] * x[i - (mid-idx)]
        end
        x_temp[i] = xtempi
    end
end

#= INTERIOR CONVOLUTION =#
function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC}
    N = length(x)
    stencil_length = length(A.up_stencil_coefs)
    stencil_rem = 1 - stencil_length%2
    Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2]
        xtempi = zero(T)
        if A.directions[][i] == false
            @inbounds for j in 1 : length(A.up_stencil_coefs)
                xtempi += A.up_stencil_coefs[j] * x[i+j-1-stencil_rem]
            end
        else
            @inbounds for j in -length(A.down_stencil_coefs)+1 : 0
                xtempi += A.down_stencil_coefs[j+stencil_length] * x[i+j+stencil_rem]
                # println("i = $i, j = $j, s_idx = $(j+stencil_length), x_idx = $(i+j+stencil_rem), $(A.down_stencil_coefs[j+stencil_length]) * $(x[i+j+stencil_rem]), xtempi = $xtempi")
            end
        end

        x_temp[i] = xtempi
    end
end

#= RIGHT BOUNDARY CONDITIONS =#
function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,LBC,:Dirichlet0}) where {T<:Real,S<:SVector,LBC}
    # Dirichlet 0 means that the value at the boundary is 0
    N = length(x)
    mid = div(A.stencil_length,2) + 1
    bpc = A.stencil_length - mid
    for i in 1 : A.boundary_point_count[2]
        dirichlet_0!(x_temp, x, A.stencil_coefs, mid, bpc, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,LBC,:Dirichlet0}) where {T<:Real,S<:SVector,LBC}
    # Dirichlet 0 means that the value at the boundary is 0
    N = length(x)
    mid = div(A.stencil_length,2) + 1
    bpc = A.stencil_length - mid
    for i in 1 : A.boundary_point_count[2]
        dirichlet_0!(x_temp, x, A.stencil_coefs[end], mid, bpc, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,:Dirichlet0}) where {T<:Real,S<:SVector,LBC}
    # Dirichlet 0 means that the value at the boundary is 0
    N = length(x)
    bpc = A.boundary_point_count[2]
    stencil_length = A.stencil_length
    stencil_rem = 1 - stencil_length%2

    for i in 1 : A.boundary_point_count[2]
        pt_idx = N - A.boundary_point_count[2] + i
        A.directions[][pt_idx] ? start_idx = stencil_length - stencil_rem : start_idx = 1 + stencil_rem
        dirichlet_1!(x_temp, x, A.directions[][pt_idx] ? A.down_stencil_coefs : A.up_stencil_coefs, start_idx, pt_idx)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,LBC,:Dirichlet}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    mid = div(A.stencil_length,2) + 1
    x[end] = A.boundary_condition[][2][3](A.t)

    for i in 1 : A.boundary_point_count[2]
        dirichlet_1!(x_temp, x, A.stencil_coefs, mid, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,LBC,:Dirichlet}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    mid = div(A.stencil_length,2) + 1
    x[end] = A.boundary_condition[][2][3](A.t)

    for i in 1 : A.boundary_point_count[2]
        dirichlet_1!(x_temp, x, A.stencil_coefs[end], mid, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,LBC,:periodic}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    for i in 1 : A.boundary_point_count[2]
        periodic!(x_temp, x, A.stencil_coefs, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,LBC,:periodic}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    for i in 1 : A.boundary_point_count[2]
        periodic!(x_temp, x, A.stencil_coefs[end], N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S,LBC,:Neumann0}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    for i in 1 : A.boundary_point_count[2]
        neumann0!(x_temp, x, A.stencil_coefs, N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::FiniteDifference{T,S,LBC,:Neumann0}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    for i in 1 : A.boundary_point_count[2]
        neumann0!(x_temp, x, A.stencil_coefs[end], N - A.boundary_point_count[2] + i)
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,LBC,:Neumann}, FiniteDifference{T,S,LBC,:Neumann}}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    @inbounds for i in 1 : A.boundary_point_count[2]
        bc = A.high_boundary_coefs[][A.boundary_point_count[2] - i + 1]
        tmp = zero(T)
        @inbounds for j in 1 : length(bc)
            # our coefficients and points are aligned so as we have not reversed anything in the stencil
            tmp += bc[j] * x[N-j+1]
        end
        x_temp[N-i+1] = tmp
    end
    x_temp[end] += A.boundary_condition[][2][3]
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,LBC,:Robin}, FiniteDifference{T,S,LBC,:Robin}}) where {T<:Real,S<:SVector,LBC}
    N = length(x)
    @inbounds for i in 1 : A.boundary_point_count[2]
        bc = A.high_boundary_coefs[][A.boundary_point_count[2] - i + 1]
        tmp = zero(T)
        @inbounds for j in 1 : length(bc)
            tmp += bc[j] * x[N-j+1]
        end
        x_temp[N-i+1] = tmp
    end
    x_temp[end] += A.boundary_condition[][2][3]
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::Union{DerivativeOperator{T,S,LBC,:None},UpwindOperator{T,S,LBC,:None}, FiniteDifference{T,S,LBC,:None}}) where {T<:Real,S<:SVector,LBC}
    # halfstencil = div(A.stencil_length, 2)
    for i in 1 : A.boundary_point_count[2] # the first stencil is for the last point ie. in reverse order
        @inbounds bc = A.high_boundary_coefs[][i]
        tmp = zero(T)
        # startid = max(0,i-1-halfstencil)
        @inbounds for j in 1 : length(bc)
            tmp += bc[j] * x[end-length(bc)+j]
        end
        @inbounds x_temp[end-i+1] = tmp
    end
end

function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,:nothing}) where {T<:Real,S<:SVector,LBC}
    stencil_length = length(A.down_stencil_coefs)
    stencil_rem = 1 - stencil_length%2
    start_idx = 1

    # this case when our stencil uses 1 point at the right, so we cant use the downwind stencil
    # as it will spill over. So we use a special boundary stencil.
    if stencil_rem == 1
        x_temp[end] = sum(A.high_boundary_coefs[][1].*x[end-length(A.high_boundary_coefs[][1])+1:end])
        start_idx = 2
    # else
    #     x_temp[end] = sum(A.down_stencil_coefs.*x[end-length(A.down_stencil_coefs)+1:end])
    end

    # println(start_idx)
    for i in start_idx : A.boundary_point_count[2] # the first stencil is for the last point ie. in reverse order
        xtempi = zero(T)
        # startid = max(0,i-1-halfstencil)
        @inbounds for j in -length(A.down_stencil_coefs)+1 : 0
            # println("$(A.down_stencil_coefs[j+stencil_length]) * $(x[end+1-i+j+stencil_rem])")
            xtempi += A.down_stencil_coefs[j+stencil_length] * x[end+1-i+j+stencil_rem]
        end
        @inbounds x_temp[end-i+1] = xtempi
    end
end

#= DIFFERENT BOUNDARIES =#
function dirichlet_0!(x_temp::AbstractVector{T}, x::AbstractVector{T}, coeffs::SVector, mid::Int, bpc::Int, i::Int) where T<:Real
    #=
        The high and low functions determine the starting and ending indices of the weight vector.
        As we move along the input vector to calculate the derivative at the pointhe weights which
        are to be considered to calculate the derivative are to be chosen carefully. eg. at the boundaries,
        only half of the stencil is going to be used to calculate the derivative at that point.
        So, observing that the left index grows as:-
                  i ^
                    |       mid = ceil(stencil_length/2)
               mid--|       bpc = boundary_point_count
                    |\
                    | \
               0  <_|__\________>
                    |  bpc      i

        And the right index grows as:-
                  i ^       mid = ceil(stencil_length/2)
                    |       bpc = boundary_point_count
               mid--|_______
                    |       \
                    |        \
               0  <_|_________\___>
                    |        bpc  i
        The high and low functions are basically equations of these graphs which are used to calculate
        the left and right index of the stencil as a function of the index i (where we need to find the derivative).
    =#
    stencil_length = length(coeffs)
    # mid = div(stencil_length, 2) + 1 # generalizing to any mid for upwind operators
    N = length(x)
    wndw_low = i>bpc ? 1 : max(1, low(i, mid, bpc))
    wndw_high = i>N-bpc ? min(stencil_length, high(i, mid, bpc, stencil_length, N)) : stencil_length

    # println(wndw_low," ",wndw_high, " mid = ", mid)
    # println("#####")

    #=
        Here we are taking the weighted sum of a window of the input vector to calculate the derivative
        at the middle point. This requires choosing the end points carefully which are being passed from above.
    =#
    xtempi = zero(T)
    @inbounds for idx in wndw_low:wndw_high
        xtempi += coeffs[idx] * x[(i - (mid-idx))]
        # println("i = $i, idx = $((i - (mid-idx))), $(coeffs[idx]) * $(x[(i - (mid-idx))]), xtempi = $xtempi")
    end
    x_temp[i] = xtempi
end

function dirichlet_1!(x_temp::AbstractVector{T}, x::AbstractVector{T}, coeffs::SVector, mid::Int, i::Int) where T<:Real
    stencil_length = length(coeffs)
    N = length(x)
    #=
        Here we are taking the weighted sum of a window of the input vector to calculate the derivative
        at the middle point. Once the stencil goes out of the edge, we assume that it has a constant
        value outside for all points.
    =#
    xtempi = zero(T)
    @inbounds for idx in 1:stencil_length
        xtempi += coeffs[idx] * x[limit(i - (mid-idx), N)]
    end
    x_temp[i] = xtempi
end

function periodic!(x_temp::AbstractVector{T}, x::AbstractVector{T}, coeffs::SVector, i::Int) where T<:Real
    stencil_length = length(coeffs)
    mid = div(stencil_length, 2) + 1
    bpc = stencil_length - mid
    wndw_low = 1
    wndw_high = length(coeffs)
    L = length(x)
    #=
        Here we are taking the weighted sum of a window of the input vector to calculate the derivative
        at the middle point. Instead of breaking the stencil we loop it over from the other side if it
        doesn't fit the column of the transformation matrix to simulate the periodic boundary condition.
    =#
    xtempi = zero(T)
    @inbounds for idx in wndw_low:wndw_high
        xtempi += coeffs[idx] * x[rem1(i - (mid-idx), L)]
    end
    x_temp[i] = xtempi
end

function neumann0!(x_temp::AbstractVector{T}, x::AbstractVector{T}, coeffs::SVector, i::Int) where T<:Real
    stencil_length = length(coeffs)
    mid = div(stencil_length, 2) + 1
    bpc = stencil_length - mid
    wndw_low = 1
    wndw_high = length(coeffs)
    L = length(x)
    #=
        Here we are taking the weighted sum of a window of the input vector to calculate the derivative
        at the middle point. Instead of breaking the stencil we loop it over from the other side if it
        doesn't fit the column of the transformation matrix to simulate the periodic boundary condition.
    =#
    xtempi = zero(T)
    @inbounds for idx in wndw_low:wndw_high
        xtempi += coeffs[idx] * x[reflect(i - (mid-idx), L)]
    end
    x_temp[i] = xtempi
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

import Dates
import TimeZones

struct WrappedDate
    x::Dates.Date
end
Arrow.ArrowTypes.registertype!(WrappedDate, WrappedDate)

struct WrappedTime
    x::Dates.Time
end
Arrow.ArrowTypes.registertype!(WrappedTime, WrappedTime)

struct WrappedDateTime
    x::Dates.DateTime
end
Arrow.ArrowTypes.registertype!(WrappedDateTime, WrappedDateTime)

struct WrappedZonedDateTime
    x::TimeZones.ZonedDateTime
end
Arrow.ArrowTypes.registertype!(WrappedZonedDateTime, WrappedZonedDateTime)

@testset "Date and time wrappers with missing" begin
    for T in (WrappedDate, WrappedTime, WrappedDateTime, WrappedZonedDateTime)
        if T == WrappedZonedDateTime
            time = T(Dates.now(TimeZones.tz"UTC"))
        else
            time = T(Dates.now())
        end
        table = (; x = [missing, missing, time, missing, time])
        io = IOBuffer()
        Arrow.write(io, table)
        seekstart(io)
        tbl = Arrow.Table(io)
        @test isequal(collect(tbl.x), table.x)
    end
end

@testset "`default(T) isa T`" begin
    for T in (Dates.Date, Dates.Time, Dates.DateTime, TimeZones.ZonedDateTime, Dates.Nanosecond, Dates.Millisecond, Dates.Second, Dates.Day, Dates.Month, Dates.Year)
        @test Arrow.ArrowTypes.default(T) isa T
    end
end

max_stars_repo_name: timkpaine/arrow

max_issues_repo_name: timkpaine/arrow

max_forks_repo_name: XpressAI/arrow


text:
using BasicStats
using Test

@testset "BasicStats.jl" begin
    @test xlogx(20)==59.914645471079815
    @test avg(1:10)==5.5
    include("fisherindex.jl")
end


text:
using BenchmarkTools

day16data = """E0525D9802FA00B80021B13E2D4260004321DC648D729DD67B2412009966D76C0159ED274F6921402E9FD4AC1B0F652CD339D7B82240083C9A54E819802B369DC0082CF90CF9280081727DAF41E6A5C1B9B8E41A4F31A4EF67E2009834015986F9ABE41E7D6080213931CB004270DE5DD4C010E00D50401B8A708E3F80021F0BE0A43D9E460007E62ACEE7F9FB4491BC2260090A573A876B1BC4D679BA7A642401434937C911CD984910490CCFC27CC7EE686009CFC57EC0149CEFE4D135A0C200C0F401298BCF265377F79C279F540279ACCE5A820CB044B62299291C0198025401AA00021D1822BC5C100763A4698FB350E6184C00A9820200FAF00244998F67D59998F67D5A93ECB0D6E0164D709A47F5AEB6612D1B1AC788846008780252555097F51F263A1CA00C4D0946B92669EE47315060081206C96208B0B2610E7B389737F3E2006D66C1A1D4ABEC3E1003A3B0805D337C2F4FA5CD83CE7DA67A304E9BEEF32DCEF08A400020B1967FC2660084BC77BAC3F847B004E6CA26CA140095003900BAA3002140087003D40080022E8C00870039400E1002D400F10038C00D100218038F400B6100229500226699FEB9F9B098021A00800021507627C321006E24C5784B160C014A0054A64E64BB5459DE821803324093AEB3254600B4BF75C50D0046562F72B1793004667B6E78EFC0139FD534733409232D7742E402850803F1FA3143D00042226C4A8B800084C528FD1527E98D5EB45C6003FE7F7FCBA000A1E600FC5A8311F08010983F0BA0890021F1B61CC4620140EC010100762DC4C8720008641E89F0866259AF460C015D00564F71ED2935993A539C0F9AA6B0786008D80233514594F43CDD31F585005A25C3430047401194EA649E87E0CA801D320D2971C95CAA380393AF131F94F9E0499A775460"""

function bint(bits, offset, len)
    return parse(Int, bits[(begin + offset):(begin + offset + len - 1)]; base=2)
end
const part = [0, 0]

function parsepacket(bits::String, idx=0)
    results, result = Int[], 0
    version = bint(bits, idx, 3)
    part[1] += version
    idx += 3
    typeid = bint(bits, idx, 3)
    idx += 3
    if typeid == 4 # literal value
        numbits, firstwas1 = "", true
        while firstwas1
            numbits *= bits[(begin + idx + 1):(begin + idx + 4)]
            firstwas1 = (bits[begin + idx] == '1')
            idx += 5
        end
        result = bint(numbits, 0, length(numbits))
    else  # or else typeid is one of the operators
        lengthtypeid = bits[begin + idx]
        idx += 1
        if lengthtypeid == '1'  # vector of nresults child packets follows
            nresults = bint(bits, idx, 11)
            idx += 11
            for _ in 1:nresults
                child, idx = parsepacket(bits, idx)
                push!(results, child)
            end
        else  # length of subresults format when length type id is 0
            subpacketlength = bint(bits, idx, 15)
            idx += 15
            subpacketend = idx + subpacketlength
            while idx < subpacketend
                child, idx = parsepacket(bits, idx)
                push!(results, child)
            end
        end
    end
    # part 2 computations on packets which depend on operator (the typeid)
    if typeid == 0
        result = sum(results)
    elseif typeid == 1
        result = prod(results)
    elseif typeid == 2
        result = minimum(results)
    elseif typeid == 3
        result = maximum(results)
    elseif typeid == 5
        result = first(results) > last(results)
    elseif typeid == 6
        result = first(results) < last(results)
    elseif typeid == 7
        result = first(results) == last(results)
    end
    return result, idx
end

function day16(s)
    part[1] = 0
    part[2], _ = parsepacket(string(parse(BigInt, s; base=16); base=2))
    return part
end

par = day16(day16data)
println("Part 1: ", par[1])
println("Part 2: ", par[2])

@btime day16(day16data)


text:
using POMDPGym, POMDPs, Test
mdp = AtariPOMDP(:Pong, version = :v0)
@test length(actions(mdp)) == 6

s = rand(initialstate(mdp))
@test size(s) == (210, 160, 3)

s, o, r = gen(mdp, s, 1)

@test preproc_atari_frame(s) isa Array{UInt8}

@test size(o) == (80,80,4,1)

@test o[1] isa UInt8
@test s[1] isa UInt8
render(mdp)

max_forks_repo_name: ancorso/POMDPGym


text:
export birth_death_pushing, birth_death_pushing!, DataFrame

########---FORMAT---########

"""
    Cell(index, position, mutations, parent, b, t_birth, p_birth)

Mutable structure used in the simulation algorithm. function. The output of the simulation
method `birth_death_pushing` is a `Vector{Cell}` and parameters.
"""
mutable struct Cell
    index :: Int64
    position :: Vector{Float64}
    mutations :: Vector{Int64}
    parent :: Int64
    b :: Float64
    t_birth :: Float64
	p_birth :: SVector{dim,Float64} where dim
end
"""
    Cell(;index, position, mutations, parent, b, t_birth, p_birth)

Constructor of type `Cell` with kwrgs for fields.
"""
Cell(;index, position, mutations, parent, b, t_birth, p_birth) = Cell(index, position, mutations, parent, b, t_birth, p_birth)

struct Mutation
	origin :: Int64
	ancestor :: Int64
	t_birth :: Float64
	N_birth :: Int64
	p_birth :: SVector{dim,Float64} where dim
end
Mutation(; origin, ancestor, t_birth, N_birth, p_birth) :: Mutation = Mutation(origin, ancestor, t_birth, N_birth, p_birth)

"""
    DataFrame(tumor::Vector{T})

Generates a DataFrame out of a Vector of type T (Cell). Columnnames of the output
are fieldnames of the type T (Cell).
"""
function DataFrame(tumor::Vector{T}) where T
    df = DataFrame()
    for field in fieldnames(T)
        ftype = fieldtype(T, field)
        df[!,field] = if ftype <: SVector
            collect.(reinterpret.(eltype(ftype),getfield.(tumor, field)))
        else
            getfield.(tumor, field)
        end
    end
    return df
end

############################
######---BIRTH RATE---######

w(r::Float64; σ=1.)  = exp(-r^2/(2*σ^2))/√(2π*σ^2)

"""
    b_linear()
Call before running the simulation to set the birth rate profile.
Linear is set by default. A custom profile can also be set by overwriting TumorGrowth.b_curve(ρ; bup, ρc).
"""
function b_linear()
	@eval b_curve(ρ::Float64; bup::Float64, ρc::Float64) = max(0., bup*(1. - ρ/ρc))
end

"""
    b_hill(n=2)
Call before running the simulation to set the birth rate profile to a hill curve o f order `n`.
A custom profile can also be set by overwriting TumorGrowth.b_curve(ρ; bup, ρc).
"""
function b_hill(n=2)
	@eval b_curve(ρ::Float64; bup::Float64, ρc::Float64) = bup/(1+exp(ρ-ρc)^($n))
end
b_linear()

"""
    update_birthrate!(cell::Cell, tumor; bup::Float64, ρ = Inf)

Determines a cells birth rate based on its surrounding from the set birth rate profile `b_curve`. For each cell in the variable `tumor` a gaussian weight (hardcoded width `σ=4`) is calculated on its distance. The sum of weights is the density that gives the cells birth rate at set denstiy threshold `ρ` and maximal birth rate `bup`.
"""
function update_birthrate!(cell::Cell, tumor; bup::Float64, ρ = Inf)
	isempty(tumor) && return cell.b = bup
    ws = getfield.(tumor, :position) .|> p -> w(norm(cell.position - p); σ=4.)
    return cell.b = b_curve(sum(ws); bup = bup, ρc = ρ)
end

#############################
#######---SIMULATION---######

"""
    birth!(tumor::Vector{Cell}, parent, cur_id, cur_mutation, μ, cellbox, t, dimv::Val{dim}) where dim
Birth event handler in simulation algorithm.
At division the birth! function is called and parent or daughter may or may not gain a new mutation.
"""
function birth!(tumor::Vector{Cell}, mutations::Vector{Mutation}, parent, cur_id, cur_mutation, μ, cellbox, t, dimv::Val{dim}) where dim

    δ = rand(dim) .- 0.5
    pos = parent.position + 2.1 .* δ/norm(δ)
    push!(cellbox, pos2box(pos, dimv))

	new = Cell(cur_id, copy(pos), copy(parent.mutations), parent.index, parent.b, t, SVector{dim,Float64}(pos))
	push!(tumor, new)

	m = 0
	if rand() < μ/2
		m += 1
        push!(mutations, Mutation(parent.index, isempty(parent.mutations) ? 0 : last(parent.mutations), t, length(tumor), SVector{dim,Float64}(parent.position)))
		push!(parent.mutations, cur_mutation+m)
	end
	if rand() < μ/2
		m += 1
		push!(mutations, Mutation(new.index, isempty(new.mutations) ? 0 : last(new.mutations), t, length(tumor), SVector{dim,Float64}(pos)))
		push!(new.mutations, cur_mutation+m)
	end
	return m
end

"""
    birth_death_pushing!( tumor::Vector{Cell}, until;
    b, d, μ, ρ=Inf, dim=length(tumor[1].position), seed=abs(rand(Int)),
    cur_id = 1, cur_mutation = 0, t = 0.0, showprogress=true)

Evolves an existing collection of cells `Vector{Cell}`.
mutates `tumor` variable and returns current `index`, `mutation` and time`t` as a `Dict`.

Rejection algorithm for density dependent tumor growth:
At each iteration a parent is uniformly chosen and its birth rate updated. It then either devides, dies (removed) or nothing happens. Time is incremented appropriatly.
For division the birth! function is called and parent or daughter may or may not gain a new mutation. The pushing! algorithm then removes overlaps.

Call `birth_death_pushing` to run a simulation from a single cell
returns:    `Dict{:tumor, :index, :mutation, :time}`
arg:        `until`       loopcondition, if Float final time, if Int final size
kwargs:     `b, d, μ, ρ=Inf, cur_mutation = 0, dim, seed=nothing`

To further evolve a tumor call `birth_death_pushing!` on existing `Cell` array and set kwargs `cur_id` (next cell index), `cur_mutation` (next mutation), `t` (time) appropriatly.
"""
function birth_death_pushing!( tumor::Vector{Cell}, mutations::Vector{Mutation}, until;
	b, d, μ, ρ=Inf, dim=length(tumor[1].position), seed=nothing,
	cur_id = 1, cur_mutation = 0, t = 0.0, showprogress=true)
    dimv = Val(dim)

    isnothing(seed) || Random.seed!(seed)

    N = length(tumor)

    cellbox = [pos2box(p, dimv) for p in getfield.(tumor,:position)]

    prog = showprogress ? ProgressUnknown("Progress: Tumor size ") : nothing
    while loop_condition(N,t, until)

        N==0 && error("Tumor died")
        row = rand(1:N)
        parent = tumor[row]

        p = rand()*(b+d) - update_birthrate!(parent, view(tumor, find_neighbors(cellbox, row; s=4)); bup=b, ρ = ρ)

        t += randexp()/((d+b)*N)

        if p < 0.
            cur_id += 1
            N += 1
            cur_mutation += birth!(tumor, mutations, parent, cur_id, cur_mutation, μ, cellbox, t, dimv)

            pushing!(tumor, N, cellbox, dimv)
        elseif p < d
            N -= 1
            deleteat!.( (tumor, cellbox), row)
        end
        showprogress && ProgressMeter.update!(prog, N )
    end
	for (i, cell) in enumerate(tumor)
	    update_birthrate!(cell, view(tumor, find_neighbors(cellbox, i; s=4)); bup=b, ρ = ρ)
	end
    return Dict{Symbol, Any}(:index => cur_id, :mutation => cur_mutation, :time => t)
end
birth_death_pushing!( until; tumor::Vector{Cell}, mutations::Vector{Mutation}, simparams...) =  birth_death_pushing!( tumor, mutations, until; simparams...)

loop_condition(N::Int,t::Float64, Nfinal::Int) = N < Nfinal
loop_condition(N::Int,t::Float64, tfinal::Float64) = t < tfinal

"""
    function birth_death_pushing( until; b, d, μ, ρ=Inf, cur_mutation = 0, dim , seed=abs(rand(Int)), showprogress=true)

Rejection algorithm for density dependent tumor growth:
At each iteration a parent is uniformly chosen and its birth rate updated. It then either devides, dies (removed) or nothing happens. Time is incremented appropriatly.
For division the birth! function is called and parent or daughter may or may not gain a new mutation. The pushing! algorithm then removes overlaps.

Call `birth_death_pushing` to run a simulation from a single cell
returns:    `Dict{:tumor, :index, :mutation, :time}`
arg:        `until`       loopcondition, if Float final time, if Int final size
kwargs:     `b, d, μ, ρ=Inf, cur_mutation = 0, dim, seed=nothing`

The output `:tumor` is a `Vector{Cell}` which can be passed to `DataFrame` for analysis or plotting.

To further evolve a tumor call `birth_death_pushing!` on existing `Cell` vector and set kwargs `cur_id` (next cell index), `cur_mutation` (next mutation), `t` (time) appropriatly.
"""
function birth_death_pushing( until; b, cur_mutation = 0, dim, simparams...)
    p₀ = zeros(Float64,dim)
    tumor = [Cell(; index=1, position=p₀, parent=0, mutations=collect(1:cur_mutation), b=b, t_birth=0.0, p_birth=SVector{dim,Float64}(p₀))]

	mutations = [Mutation(; origin=1, ancestor=0, t_birth=0., N_birth=1, p_birth=SVector{dim,Float64}(p₀)) for m=1:cur_mutation]

    output = birth_death_pushing!(tumor, mutations, until; b=b, cur_mutation=cur_mutation, dim=dim, simparams...)
	output[:tumor] = tumor
	output[:mutations] = mutations

	return output
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
using Test
using CartesianGeneticProgramming
import YAML

# Configuration file used for all tests
test_filename = string(@__DIR__, "/test.yaml")

function test_ind(ind::CGPInd, cfg::NamedTuple)
    @test length(ind.nodes) == cfg.rows * cfg.columns + cfg.n_in
    @test length(ind.chromosome) == cfg.rows * cfg.columns * (3 + cfg.n_parameters) + cfg.n_out
    @test size(ind.genes) == (cfg.rows, cfg.columns, 3 + cfg.n_parameters)
    for node in ind.nodes
        if node.active
            @test node.x >= 1
            @test node.x <= length(ind.nodes)
            @test node.y >= 1
            @test node.y <= length(ind.nodes)
        end
        if node.x == 0 && node.y == 0 # input node, no parameters
            @test length(node.p) == 0
        else # non-input node, test parameters
            @test length(node.p) == cfg.n_parameters
            @test all(0.0 .<= node.p .<= 1.0)
        end
        # test that stringifying works
        @test typeof(string(node)) == String
    end
    @test typeof(ind.genes) == Array{Float64,3}
    # assert that genes encoding function, x and y are integers
    @test ind.genes[:, :, 1:3] == Int16.(ind.genes[:, :, 1:3])
    # assert that genes encoding parameters are floating numbers in [0, 1]
    if cfg.n_parameters > 0
        @test all(0.0 .<= ind.genes[:, :, 4:end] .<= 1.0)
    end
end

@testset "CGPInd construction" begin
    cfg = get_config(test_filename)
    ind = CGPInd(cfg)
    test_ind(ind, cfg)
end

# TODO: remove or convert to new fgen
"""
A minimal function module example.
Note that one can provide any function names, these are just to keep consistency
with the `test.yaml` configuration file.
The `arity` dictionary is necessary though.
"""
module MinimalFunctionModuleExample
    global arity = Dict()
    function fgen(name::Symbol, ar::Int, s1::Union{Symbol, Expr})
        @eval function $name(x::Float64, y::Float64, p::Array{Float64})::Float64
            $s1
        end
        arity[String(name)] = ar
    end
    fgen(:f_add, 2, :(p[1] * 0.5 * (x + y)))
    fgen(:f_subtract, 2, :(p[1] * (x - y)))
    fgen(:f_mult, 2, :(p[1] * x * y))
    fgen(:f_div, 1, :(p[1] * x / 2))
    fgen(:f_abs, 2, :(p[1] * abs(x)))
end

# Similar to CGPInd construction but uses custom set of CGP functions
@testset "CGPInd construction with custom functions" begin
    cfg = get_config(test_filename, function_module=MinimalFunctionModuleExample)
    ind = CGPInd(cfg)
    test_ind(ind, cfg)
end

# Similar to CGPInd construction but uses a custom buffer
@testset "CGPInd construction with custom buffer" begin
    cfg = get_config(test_filename)
    my_buffer = zeros(Int64, cfg.rows * cfg.columns + cfg.n_in)
    ind = CGPInd(cfg; buffer=my_buffer)
    test_ind(ind, cfg)
    @test typeof(ind.buffer) == Array{Int64,1}
end

"""
using random values, sort individuals that are different
"""
function select_random(pop::Array{CGPInd}, elite::Int; n_in=113, n_sample=100)
    actions = zeros(Int, length(pop))
    dists = zeros(n_sample, length(pop))
    inputs = rand(n_in, n_sample)

    for i in 1:n_sample
        for j in eachindex(pop)
            actions[j] = argmax(process(pop[j], inputs[:, i]))
        end
        for j in eachindex(pop)
            dists[i, j] = sum(actions[j] .!= actions)
        end
    end
    d = sum(dists, dims=1)[:]
    ds = sortperm(d)[1:elite]
    pop[ds]
end

@testset "Node genes" begin
    cfg = get_config(test_filename)
    ind = CGPInd(cfg)
    ind2 = CGPInd(cfg)
    @test any(ind.chromosome .!= ind2.chromosome)

    for i in 1:ind.n_in
        genes = get_genes(ind, i)
        @test all(genes .== 0)
    end

    for i in 1:(length(ind.nodes)-ind.n_in)
        genes = get_genes(ind, ind.n_in+i)
        @test all(0.0 .<= genes .<= 1.0)
    end

    for i in 1:(length(ind.nodes)-ind.n_in)
        ind2 = CGPInd(cfg)  # re-init
        genes = get_genes(ind, ind.n_in+i)
        set_genes!(ind2, ind.n_in+i, genes)
        @test any(ind.chromosome .== ind2.chromosome)
    end

    for i in 1:length(ind.nodes)
        genes = get_genes(ind, i)
        set_genes!(ind2, i, genes)
    end
    o = length(ind.chromosome) - cfg.n_out
    @test all(ind.chromosome[1:o] .== ind2.chromosome[1:o])

    all_genes = get_genes(ind, collect((ind.n_in+1):length(ind.nodes)))
    # @test all(ind.chromosome[1:o] .== all_genes)
    for g in all_genes
        @test g in ind.chromosome
    end
end

@testset "Processing" begin
    cfg = get_config(test_filename; functions=["f_abs", "f_add", "f_mult"])
    ind = CGPInd(cfg)

    # test that f(0, 0, 0, 0) = 0
    inputs = zeros(cfg.n_in)
    set_inputs(ind, inputs)
    for i in 1:cfg.n_in
        @test ind.buffer[i] == 0.0
    end
    output = process(ind)
    @test output[1] == 0.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] == 0.0
        end
    end

    # test that f(1, 1, 1, 1) = 1
    for i in eachindex(ind.nodes)
        ind.buffer[i] = 1.0 # requires that buffer is 1
    end
    output = process(ind, ones(cfg.n_in))
    @test output[1] == 1.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] == 1.0
        end
    end

    cfg = get_config(test_filename; functions=["f_abs", "f_add", "f_mult"], recur=1.0)
    ind = CGPInd(cfg)
    output = process(ind, rand(cfg.n_in))
    @test output[1] <= 1.0 && output[1] >= -1.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] <= 1.0 && ind.buffer[i] >= -1.0
        end
    end

    # test parametric functions
    cfg = get_config(test_filename; functions=["f_add"], n_parameters=1)
    ind = CGPInd(cfg)
    output = process(ind, rand(cfg.n_in))
    @test typeof(output[1]) == Float64
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test typeof(ind.buffer[i]) == Float64
        end
    end

    pop = [CGPInd(cfg) for i in 1:10]
    sp = select_random(pop, 2; n_in=cfg.n_in, n_sample=5)
    @test length(sp) == 2
    @test sp[1].buffer[1] != 0.0

    reset!(ind)
    @test all(ind.buffer .== 0.0)

    f_conns = forward_connections(ind)
    @test any(map(x->length(x)>1, f_conns))

    ot = get_output_trace(ind, 1)
    @test length(ot) > 0

    all_traces = get_output_trace(ind)
    @test issubset(ot, all_traces)
end


text:
using ControlBenchmarks
using ControlSystems

benchmarkProb = controlbenchmark( BallOnPlate() )

@test ControlSystems.nstates( benchmarkProb.sys ) == 2
@test ControlSystems.ninputs( benchmarkProb.sys ) == 1


text:
#= src/util.jl
=#

crandn(::Type{T}, dims::Integer...) where {T<:Real} = convert(Array{T}, randn(dims...))
for (elty, relty) in ((:ComplexF32, :Float32), (:ComplexF64, :Float64))
  @eval begin
    crandn(::Type{$elty}, dims::Integer...) =
      reshape(reinterpret($elty, vec(crandn($relty, 2*dims[1], dims[2:end]...))),
              (dims...,))
  end
end

function findmaxabs(x::AbstractVecOrMat{T}) where T
  m  = zero(real(T))
  mi = 0
  @inbounds for i = 1:length(x)
    t = abs(x[i])
    t < m && continue
    m  = t
    mi = i
  end
  m, mi
end

function getcols(trans::Symbol, A::AbstractMatrix, cols)
  if     trans == :n  return A[:,cols]
  elseif trans == :c  return collect(A[cols,:]')
  end
end
function getcols(trans::Symbol, A::AbstractLinearOperator, cols)
  if     trans == :n  return A[:,cols]
  elseif trans == :c  return collect(A'[:, cols])
  end
end

function hermitianize!(A::AbstractMatrix, uplo::Symbol=:U)
  checksquare(A)
  n = size(A, 2)
  if uplo == :U
    for j = 1:n, i = 1:j
      @inbounds A[i,j] = 0.5*(A[i,j] + conj(A[j,i]))
    end
  elseif uplo == :L
    for i = 1:n, j = 1:i
      @inbounds A[i,j] = 0.5*(A[i,j] + conj(A[j,i]))
    end
  end
  Hermitian(A, uplo)
end

function iscale!(A::AbstractMatrix, b::AbstractVector)
  m, n = size(A)
  length(b) == n || throw(DimensionMismatch)
  @inbounds for j = 1:n
    bj = b[j]
    @simd for i = 1:m
      A[i,j] /= bj
    end
  end
  A
end
function iscale!(b::AbstractVector, A::AbstractMatrix)
  m, n = size(A)
  length(b) == m || throw(DimensionMismatch)
  @inbounds for j = 1:n
    @simd for i = 1:m
      A[i,j] /= b[i]
    end
  end
  A
end

function orthcols!(
    A::AbstractMatrix{T}, tau::Vector{T}, work::Vector{T}; thin::Bool=true) where T<:BlasFloat
  m, n = size(A)
  k = min(m, n)
  A, tau, work = _LAPACK.geqrf!(A, tau, work)
  A, tau, work = _LAPACK.orgqr!(A, tau, k, work)
  if thin  A = k < n ? A[:,1:k] : A
  else     A[:,k+1:end] .= 0
  end
  A, tau, work
end
orthcols!(A::AbstractMatrix{T}; thin::Bool=true) where {T} =
  orthcols!(A, Array{T}(undef, 1), Array{T}(undef, 1), thin=thin)[1]

function orthrows!(
    A::AbstractMatrix{T}, tau::Vector{T}, work::Vector{T}; thin::Bool=true) where T<:BlasFloat
  m, n = size(A)
  k = min(m, n)
  A, tau, work = _LAPACK.gelqf!(A, tau, work)
  A, tau, work = _LAPACK.orglq!(A, tau, k, work)
  if thin  A = k < m ? A[1:k,:] : A
  else     A[k+1:end,:] .= 0
  end
  A, tau, work
end
orthrows!(A::AbstractMatrix{T}; thin::Bool=true) where {T} =
  orthrows!(A, Array{T}(undef, 1), Array{T}(undef, 1), thin=thin)[1]

function scalevec!(s::AbstractVector, x::AbstractVector)
  n = length(x)
  length(s) == n || throw(DimensionMismatch)
  @inbounds @simd for i = 1:n
    x[i] *= s[i]
  end
  x
end
function iscalevec!(s::AbstractVector, x::AbstractVector)
  n = length(x)
  length(s) == n || throw(DimensionMismatch)
  @inbounds @simd for i = 1:n
    x[i] /= s[i]
  end
  x
end

symrelerr(x, y) = 2*abs((x - y)/(x + y))