text:
using JuMP
using CPLEX

include("parser.jl")

start = time()

function Mymodel(MyFileName::String)
  n, P, alpha, Survie, indice_danger, cout = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, x[1:n,1:n], Bin)
  @variable(m, y[1:n,1:n], Bin)
  @variable(m, z[1:P,1:n,1:n] <= 0)

  ## Constraints
  @constraint(m, [k in 1:P, i in 1:n, j in 1:n ; indice_danger[k] == 0], z[k,i,j] >= x[i,j]*log(1-Survie[k,i,j]))
  @constraint(m, [k in 1:P, i in 1:n, j in 1:n ; indice_danger[k] == 1], z[k,i,j] >= y[i,j]*log(1-Survie[k,i,j]))
  @constraint(m, [k in 1:P], sum(z[k,i,j] for i in 1:n, j in 1:n) <= log(1-alpha[k]))
  @constraint(m, [i in 1:n, j in 1:n ; i == 1 || j == 1 || i == n || j == n], y[i] <= 0)    # bordure => y = 0
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i,j])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i-1,j])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i+1,j])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i,j-1])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i,j+1])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i-1,j-1])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i-1,j+1])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i+1,j-1])
  @constraint(m, [i in 2:n-1, j in 2:n-1], y[i,j] <= x[i+1,j+1])

  ## Objective
  @objective(m, Min, sum(x[i,j]*cout[i,j] for i in 1:n, j in 1:n))

  #resolution
  optimize!(m)

  # solution_summary(m, verbose=true)
  z_vals = JuMP.getvalue.( m[:z] )
  # println(z_vals)

  for k in 1:P
    println(-(exp(sum(z_vals[k,i,j] for i in 1:n, j in 1:n))-1))
  end

  # println(-(exp(sum(z[k,i,j] for i in 1:n, j in 1:n))-1))

end

Mymodel("outputn10P6.txt")

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
function gpu_correlate(
    correlator::EarlyPromptLateCorrelator,
    downconverted_signal_re::Array,
    downconverted_signal_im::Array,
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where NC
    late = zero(Complex{Int32})
    prompt = zero(Complex{Int32})
    early = zero(Complex{Int32})
    @inbounds for i = start_sample:num_samples_left + start_sample - 1
        late = late + (downconverted_signal_re[i] + 1im * downconverted_signal_im[i]) * code[i]
    end
    @inbounds for i = start_sample:num_samples_left + start_sample - 1
        prompt = prompt + (downconverted_signal_re[i] + 1im * downconverted_signal_im[i]) * code[i + early_late_sample_shift]
    end
    @inbounds for i = start_sample:num_samples_left + start_sample - 1
        early = early + (downconverted_signal_re[i] + 1im * downconverted_signal_im[i]) * code[i + 2 * early_late_sample_shift]
    end
    EarlyPromptLateCorrelator(
        Tracking.get_early(correlator) + early,
        Tracking.get_prompt(correlator) + prompt,
        Tracking.get_late(correlator) + late
    )
end

function correlate(
    correlator::EarlyPromptLateCorrelator{<: SVector{N}},
    downconverted_signal::StructArray{Complex{Float32},2,NamedTuple{(:re, :im),Tuple{Array{Float32,2},Array{Float32,2}}},Int64},
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where {N, NC}
    late = zero(MVector{N, Complex{Int32}})
    prompt = zero(MVector{N, Complex{Int32}})
    early = zero(MVector{N, Complex{Int32}})
    @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1
        late[j] = late[j] + downconverted_signal[i,j] * code[i]
    end
    @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1
        prompt[j] = prompt[j] + downconverted_signal[i,j] * code[i + early_late_sample_shift]
    end
    @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1
        early[j] = early[j] + downconverted_signal[i,j] * code[i + 2 * early_late_sample_shift]
    end
    EarlyPromptLateCorrelator(
        get_early(correlator) + early .* agc_attenuation / 1 << (agc_bits + NC),
        get_prompt(correlator) + prompt .* agc_attenuation / 1 << (agc_bits + NC),
        get_late(correlator) + late .* agc_attenuation / 1 << (agc_bits + NC)
    )
end

function sgpu_correlate(
    correlator::EarlyPromptLateCorrelator,
    downconverted_signal_re::CuArray,
    downconverted_signal_im::CuArray,
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where NC
    late = zero(ComplexF32)
    prompt = zero(ComplexF32)
    early = zero(ComplexF32)
    @views late = dot(
                (downconverted_signal_re[start_sample:num_samples_left + start_sample - 1]
                + 1im*downconverted_signal_im[start_sample:num_samples_left + start_sample - 1]),
                code[start_sample:num_samples_left + start_sample - 1])
    @views prompt = dot(
                (downconverted_signal_re[start_sample:num_samples_left + start_sample - 1]
                + 1im*downconverted_signal_im[start_sample:num_samples_left + start_sample - 1]),
                code[(start_sample:num_samples_left + start_sample - 1) .+ early_late_sample_shift])
    @views early = dot(
                (downconverted_signal_re[start_sample:num_samples_left + start_sample - 1]
                + 1im*downconverted_signal_im[start_sample:num_samples_left + start_sample - 1]),
                code[(start_sample:num_samples_left + start_sample - 1) .+ 2 * early_late_sample_shift])
    EarlyPromptLateCorrelator(
        Tracking.get_early(correlator) + early,
        Tracking.get_prompt(correlator) + prompt,
        Tracking.get_late(correlator) + late
    )
end

# StructArray GPU
function sgpu_correlate(
    correlator::EarlyPromptLateCorrelator{<: SVector{N}},
    downconverted_signal_re::CuMatrix,
    downconverted_signal_im::CuMatrix,
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where {N, NC}
    late = zero(MVector{N, ComplexF32})
    prompt = zero(MVector{N, ComplexF32})
    early = zero(MVector{N, ComplexF32})
    for i in 1:N
        late[i] = @views dot(
            code[start_sample:num_samples_left + start_sample - 1],(
            downconverted_signal_re[(start_sample:num_samples_left + start_sample - 1),i] +
            1im*downconverted_signal_im[(start_sample:num_samples_left + start_sample - 1),i]))
    end
    for i in 1:N
        prompt[i] = @views dot(
            code[start_sample+early_late_sample_shift:num_samples_left + start_sample + early_late_sample_shift - 1],(
            downconverted_signal_re[(start_sample:num_samples_left + start_sample - 1),i] +
            1im*downconverted_signal_im[(start_sample:num_samples_left + start_sample - 1),i]))
    end
    for i in 1:N
        early[i] = @views dot(
            code[start_sample+2*early_late_sample_shift:num_samples_left + start_sample + 2*early_late_sample_shift - 1],(
            downconverted_signal_re[(start_sample:num_samples_left + start_sample - 1),i] +
            1im*downconverted_signal_im[(start_sample:num_samples_left + start_sample - 1),i]))
    end
    EarlyPromptLateCorrelator(
        Tracking.get_early(correlator) .+ early,
        Tracking.get_prompt(correlator) .+ prompt,
        Tracking.get_late(correlator) .+ late
    )
end

# CuArray Ant=1
function gpu_correlate(
    correlator::EarlyPromptLateCorrelator,
    downconverted_signal::CuArray{ComplexF32},
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where NC
    late = zero(ComplexF32)
    prompt = zero(ComplexF32)
    early = zero(ComplexF32)
    @views late = downconverted_signal[start_sample:num_samples_left + start_sample - 1] ⋅ code[start_sample:num_samples_left + start_sample - 1]
    @views prompt = downconverted_signal[start_sample:num_samples_left + start_sample - 1] ⋅ code[early_late_sample_shift .+ (start_sample:num_samples_left + start_sample - 1)]
    @views early = downconverted_signal[start_sample:num_samples_left + start_sample - 1] ⋅ code[2*early_late_sample_shift .+ (start_sample:num_samples_left + start_sample - 1)]
    EarlyPromptLateCorrelator(
        Tracking.get_early(correlator) + early,
        Tracking.get_prompt(correlator) + prompt,
        Tracking.get_late(correlator) + late
    )
end

# Mulit Antenne CuArray
function gpu_correlate(
    correlator::EarlyPromptLateCorrelator{<: SVector{N}},
    downconverted_signal::CuArray{Complex{Float32},2},
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
)  where {N, NC}
    late = zero(MVector{N, ComplexF32})
    prompt = zero(MVector{N, ComplexF32})
    early = zero(MVector{N, ComplexF32})
    for i in 1:N
        late[i] = @views dot(code[start_sample:num_samples_left + start_sample - 1],
        downconverted_signal[(start_sample:num_samples_left + start_sample - 1),i])
    end
    for i in 1:N
        prompt[i] = @views dot(code[start_sample+early_late_sample_shift:num_samples_left + start_sample + early_late_sample_shift - 1],
        downconverted_signal[(start_sample:num_samples_left + start_sample - 1),i])
    end
    for i in 1:N
        early[i] = @views dot(code[start_sample+2*early_late_sample_shift:num_samples_left + start_sample + 2*early_late_sample_shift - 1],
        downconverted_signal[(start_sample:num_samples_left + start_sample - 1),i])
    end
    EarlyPromptLateCorrelator(
        Tracking.get_early(correlator) .+ early,
        Tracking.get_prompt(correlator) .+ prompt,
        Tracking.get_late(correlator) .+ late
    )
end

function gpu_correlate(
    correlator::EarlyPromptLateCorrelator,
    downconverted_signal::StructArray,
    code,
    early_late_sample_shift,
    start_sample,
    num_samples_left,
    agc_attenuation,
    agc_bits,
    carrier_bits::Val{NC}
) where NC
    gpu_correlate(
        correlator,
        downconverted_signal.re,
        downconverted_signal.im,
        code,
        early_late_sample_shift,
        start_sample,
        num_samples_left,
        agc_attenuation,
        agc_bits,
        carrier_bits
    )
end

ext: jl


text:
export iauPvxpv
"""
Outer (=vector=cross) product of two pv-vectors.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a        double[2][3]      first pv-vector
   b        double[2][3]      second pv-vector

Returned:
   axb      double[2][3]      a x b

Notes:

   1. If the position and velocity components of the two pv-vectors are
   ( ap, av ) and ( bp, bv ), the result, a x b, is the pair of
   vectors ( ap x bp, ap x bv + av x bp ).  The two vectors are the
   cross-product of the two p-vectors and its derivative.

   2. It is permissible to re-use the same array for any of the
   arguments.

Called:
   iauCpv       copy pv-vector
   iauPxp       vector product of two p-vectors
   iauPpp       p-vector plus p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvxpv(double a[2][3], double b[2][3], double axb[2][3])
function iauPvxpv(a::AbstractMatrix{<:Real}, b::AbstractMatrix{<:Real})
   # Preallocate return values
   axb = zeros(Float64, 3, 2)

   ccall((:iauPvxpv, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Array{Float64, 2}, a'),
          convert(Array{Float64, 2}, b'),
          axb)

   return SMatrix{2,3}(axb')
end


text:
@testset "boundary conditions" begin

    @testset "Circle Particle" begin

        ωs = [0.1,0.2,0.3]

        Nh = 8
        basis_order = Nh
        medium = Acoustic(1.,1.,2)

        # Choose particles
        soft_medium = Acoustic(0.,0.0 + 0.0im,2)
        p_soft = Particle(soft_medium, Sphere([1.0,2.0], .5))

        hard_medium = Acoustic(Inf,Inf + 0.0im,2)
        p_hard = Particle(hard_medium, Sphere([-3.0,-2.0], 0.3))

        sound = Acoustic(medium.ρ, 4. + 0.0im,2)
        p1 = Particle(sound, Sphere([-10.0,0.0], .2))

        particles = [p_soft, p_hard, p1]

        # Create two point sources
        source_position = SVector(0.0,0.2)
        amplitude = 1.0
        source1 = point_source(medium, source_position, amplitude);
        source2 = plane_source(medium, SVector(0.0,0.0), SVector(1.0,0.0), amplitude);
        # source2 = point_source(medium, -source_position, amplitude)
        source = 0.5*source1 + 0.5*source2;

        sim = FrequencySimulation(particles, source)
        sim_source = FrequencySimulation(source)

        result = run(sim_source, SVector(1.0,2.0), 0.1)

        pressure_results, displace_results =  boundary_data(
            shape(particles[1]), particles[1].medium, medium, sim, ωs;
            min_basis_order = 8, basis_order = 16
        )
        pressure_source_results, displace_source_results =  boundary_data(
            shape(particles[1]), particles[1].medium, medium, sim_source, ωs;
            min_basis_order = 6,
            basis_order = 10
        )

        # Zero pressure (Dirichlet) boundary condition
        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) < 1e-6 * mean(norm.(pressure_source_results[2].field))

        pressure_results, displace_results =  boundary_data(
            shape(particles[2]), particles[2].medium, medium, sim, ωs;
            min_basis_order = 8,
            basis_order = 16
        )
        pressure_source_results, displace_source_results =  boundary_data(
            shape(particles[2]), particles[2].medium, medium, sim_source, ωs;
            basis_order = 6,
            min_basis_order = 8
        );

        # Zero displacement (Neuman) boundary condition
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) < 5e-5 * mean(norm.(displace_source_results[2].field))

        pressure_results, displace_results =  boundary_data(shape(particles[3]), particles[3].medium, medium, sim, ωs; min_basis_order = 8, basis_order = 14, dr = 8e-6);
        pressure_source_results, displace_source_results =  boundary_data(shape(particles[3]), particles[3].medium, medium, sim_source, ωs; min_basis_order = 6, basis_order = 10, dr = 1e-6);

        # Continuous pressure and displacement accross particl boundary
        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) < 2e-6 * mean(norm.(pressure_source_results[2].field))
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) < 6e-5 * mean(norm.(displace_source_results[1].field))

        # The source pressure should always be continuous across any interface, however the displacement is only continuous because p1.medium.ρ == medium.ρ
        @test mean(norm.(pressure_source_results[1].field - pressure_source_results[2].field)) < 4e-8 * mean(norm.(pressure_source_results[2].field))
        @test mean(norm.(displace_source_results[1].field - displace_source_results[2].field)) < 5e-7 * mean(norm.(displace_source_results[1].field))
    end

    @testset "CapsuleParticle" begin

        concen_particles1 = [
            Particle(Acoustic(2.0,2.0,2),Sphere((0.0,0.0),1.0)),
            Particle(Acoustic(1.0,1.,2),Sphere((0.0,0.0),2.0))
        ]
        concen_particles2 = [
            Particle(Acoustic(2.0,1.0,2),Sphere((6.0,3.0),1.)),
            Particle(Acoustic(0.4,0.8,2),Sphere((6.0,3.0),1.4))
        ]
        particle = Particle(Acoustic(1.5,1.5,2),Sphere((-6.0,-3.0),0.8))

        ps = [CapsuleParticle(concen_particles2...), CapsuleParticle(concen_particles1...), particle]

        medium = Acoustic(0.8, 0.5 + 0.0im,2)
        source = plane_source(medium, SVector(0.0,0.0), SVector(1.0,0.0), 1.)
        sim = FrequencySimulation(ps, source)

        ωs = [0.01,0.2,0.3,1.]
        basis_vec = [8,16,16,24]

        pressure_results, displace_results =  boundary_data(shape(ps[1].outer), ps[1].outer.medium, medium, sim, ωs; basis_order_vec = basis_vec, dr = 1e9 * eps(Float64))

        # Continuous pressure and displacement accross particl boundary
        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) / mean(norm.(pressure_results[2].field)) < 1e-6
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) / mean(norm.(displace_results[2].field)) < 5e-6

        pressure_results, displace_results =  boundary_data(shape(ps[1].inner), ps[1].inner.medium, ps[1].outer.medium, sim, ωs;  basis_order_vec = basis_vec, dr = 1e9 * eps(Float64))
        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) / mean(norm.(pressure_results[2].field)) < 1e-6
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) / mean(norm.(displace_results[2].field)) < 5e-6

        pressure_results, displace_results =  boundary_data(shape(ps[2].outer), ps[2].outer.medium, medium, sim, ωs; basis_order_vec = basis_vec .+ 2, dr = 1e9 * eps(Float64))
        # Continuous pressure and displacement accross particl boundary
        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) / mean(norm.(pressure_results[2].field)) < 5e-6
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) / mean(norm.(displace_results[2].field)) < 1e-5

        pressure_results, displace_results =  boundary_data(shape(ps[3]), ps[3].medium, medium, sim, ωs; basis_order_vec = basis_vec, dr = 1e9 * eps(Float64))

        @test maximum(norm.(pressure_results[1].field - pressure_results[2].field)) / mean(norm.(pressure_results[2].field)) < 1e-6
        @test maximum(norm.(displace_results[1].field - displace_results[2].field)) / mean(norm.(displace_results[2].field)) < 5e-6
    end

    @testset "Spherical Particles" begin
        # ωs = [0.1,0.2,0.3]
        ωs = [0.3]
        medium = Acoustic(3,ρ=1.1,c=0.5)

        # Choose particles
        soft_medium = Acoustic(3; ρ = 0.2, c= 0.3 + 0.0im)
        p1 = Particle(soft_medium,Sphere([1.0,2.0,0.0], .5))

        hard_medium = Acoustic(3; ρ = 3.2, c= 6.3 + 0.0im)
        p2 = Particle(hard_medium,Sphere([-1.0,-2.0,0.0], 0.3))

        sound = Acoustic(3,ρ=medium.ρ,c=1.5)
        p3 = Particle(sound,Sphere([-4.0,0.0,0.0], .2))

        particles = [p1, p2, p3]

        os = origin.(particles)
        [norm(o-u) for o in os, u in os]

        # Create two point sources
        # source_position = SVector(0.0,0.2)
        amplitude = 1.0
        # source1 = point_source(medium, source_position, amplitude);
        source2 = plane_source(medium, SVector(0.0,0.0,0.0), SVector(0.0,0.0,1.0), amplitude);
        # source2 = point_source(medium, -source_position, amplitude)
        # source = 0.5*source1 + 0.5*source2;
        source = 1.5*source2;

        sim = FrequencySimulation(particles, source);
        sim_source = FrequencySimulation(source);

        # result = run(sim_source, SVector(1.0,2.0,1.0), ωs)

        dr = 1e-7
        map(particles) do p
            pressure_results, displace_results =  boundary_data(
                shape(p), p.medium, medium, sim, ωs;
                dr = dr, basis_order = 8 #,min_basis_order = 5
            );
            pressure_source_results, displace_source_results =  boundary_data(
                shape(p), p.medium, medium, sim_source, ωs;
                dr = dr, basis_order = 8 #,min_basis_order = 5
            );

            @test maximum(norm.(pressure_results[1].field - pressure_results[2].field) ./ norm.(pressure_results[2].field)) < 1e-7
            @test maximum(norm.(displace_results[1].field - displace_results[2].field) ./ norm.(displace_results[2].field)) < 4e-6
        end
    end
end


text:
using Turing, Distributions

#####
##### Sampling from posteriors
#####

"""
Generate test realizations of the gradient from the posterior of the system
    parameters. No integration is performed
"""
function sample_posterior_grad(chain, dat, sample_ind, vars, system_model;
                                num_samples=100, t = [0])
    param_samples = sample(chain, num_samples)
    all_vals = zeros(num_samples, length(sample_ind), size(dat,1))
    all_noise = zeros(num_samples)
    for i in 1:num_samples
        these_params = [param_samples[v].value[i] for v in vars]
        all_noise[i] = param_samples[:noise].value[i]
        for (i_save, i_dat) in enumerate(sample_ind)
            all_vals[i, i_save, :] = system_model(dat[:,i_dat],
                                    these_params, t)
        end
    end

    return (sample_trajectories=all_vals, sample_noise=all_noise)
end

"""
Generate test realizations of the gradient from the posterior of the system
    parameters. No integration is performed
"""
function sample_posterior_trajectories(chain, dat, sample_ind, vars, system_model,
                                tspan, ts;
                                num_samples=100)
    param_samples = sample(chain, num_samples)
    all_vals = []
    all_noise = zeros(num_samples)
    for i in 1:num_samples
        these_params = [param_samples[v].value[i] for v in vars]
        all_noise[i] = param_samples[:noise].value[i]
        prob = ODEProblem(system_model, dat[:,sample_ind[1]], tspan, these_params)
        sol = solve(prob, Tsit5(), saveat=ts);
        push!(all_vals, Array(sol))
    end

    return (sample_trajectories=all_vals, sample_noise=all_noise)
end

#####
##### Generating the posteriors via NUTS
#####

"""
Generate a chain with settings that have worked well for me
    Uses NUTS sampler

Input:
    dat - Each row is a variable
    numerical_grad - Each row is a variable
    turing_model - a model which takes in the data and predicts the gradient
        e.g:
        turing_model(grad, data)
"""
function generate_chain(dat, numerical_grad, turing_model;
                        opt=NUTS,
                        iterations=1000,
                        num_training_pts=500, start_ind=1,
                        train_ind=nothing)
    n_adapts = Int(minimum([iterations/5, 100]))
    j_max = 1.0
    # Try to predict the GRADIENT from data
    if train_ind == nothing
        train_ind = start_ind:num_training_pts+start_ind-1
    end
    y = numerical_grad[:,train_ind]
    x = dat[:,train_ind]
    # Turing 0.6.* syntax
    # chain = sample(turing_model(y, x),
    #                 NUTS(iterations, n_adapts, 0.6j_max));
    # Turing 0.7 syntax
    chain = sample(turing_model(y, x),
                    NUTS(n_adapts, 0.6j_max), iterations);
    return (chain=chain, train_ind=train_ind)
end

export sample_posterior_trajectories, sample_posterior_grad,
    generate_chain


text:
module TestCopy

include("../preamble.jl")
import Tables
using DataFrames: DataFrame, eachrow
using StructArrays: StructVector, StructArray
using TypedTables: Table

@testset "copy [base types]" begin
    @test copy(Map(x -> x^2), Vector, 1:3) == (1:3) .^ 2
    @test copy(Map(x -> x^2), [1:3;]) == (1:3) .^ 2
    @test copy(Map(x -> x^2), Set, 1:3) == Set((1:3) .^ 2)
    @test copy(Map(x -> x => x^2), Dict, 2:2) == Dict(2 => 4)
end

@testset "$copy" for copy in [copy, tcopy, dcopy]
    @testset "$copy(_, ::$(prettytypeof(src)))" for src in Any[
        DataFrame(a=[1:4;], b=[5:8;]),
        StructVector(a=[1:4;], b=[5:8;]),
        Table(a=[1:4;], b=[5:8;]),
    ]
        if copy in (tcopy, dcopy) && src isa DataFrame
            @test_broken copy(Map(identity), src) ==ₜ src
            @test_broken copy(Map(identity), src; basesize=1) ==ₜ src
            continue
        end
        @test copy(Map(identity), src) ==ₜ src
        @test copy(eduction(identity(x) for x in src)) ==ₜ src
        if copy in (tcopy, dcopy)
            @test copy(Map(identity), src; basesize=1) ==ₜ src
            @test copy(identity(x) for x in src) ==ₜ src
            @test copy((identity(x) for x in src); basesize=1) ==ₜ src
        end
    end
    @testset "$copy(_, $T, ::$(prettytypeof(src)))" for src in Any[
            StructVector(a = [1:4;], b = [5:8;]),
            Table(a = [1:4;], b = [5:8;]),
        ],
        T in Any[DataFrame, StructArray, Table]

        @test copy(Map(identity), T, src) ==ₜ T(src)
        @test copy(T, eduction(identity(x) for x in src)) ==ₜ T(src)
        if copy in (tcopy, dcopy)
            @test copy(Map(identity), T, src; basesize = 1) ==ₜ T(src)
            @test copy(T, identity(x) for x in src) ==ₜ T(src)
            @test copy(T, (identity(x) for x in src); basesize = 1) ==ₜ T(src)
        end
    end
    @testset "$copy(_, eachrow(df))" begin
        df = DataFrame(a=[1:4;], b=[5:8;])
        if Tables.materializer(eachrow(df)) == DataFrame
            # requires https://github.com/JuliaData/DataFrames.jl/pull/2055
            if copy in (tcopy, dcopy)
                @test copy(Map(identity), eachrow(df); basesize=1) ==ₜ df
                @test copy((identity(x) for x in eachrow(df)); basesize=1) ==ₜ df
                @test copy(eduction(identity(x) for x in eachrow(df)); basesize=1) ==ₜ df
            else
                @test copy(Map(identity), eachrow(df)) ==ₜ df
                @test copy(eduction(identity(x) for x in eachrow(df))) ==ₜ df
            end
        else
            @test_broken copy(Map(identity), eachrow(df)) ==ₜ df
        end
    end
end

@testset "$collect" for collect in [collect, tcollect, dcollect]
    @testset "$collect(_, ::$(prettytypeof(src)))" for src in Any[
        1:64,
        StructVector(a=[1:4;], b=[5:8;]),
     ]
        @test collect(Map(identity), src) ==ₜ Base.collect(src)
        @test collect(IdentityTransducer(), src) ==ₜ Base.collect(src)
        @test collect(eduction(identity(x) for x in src)) ==ₜ Base.collect(src)
        if collect in (tcollect, dcollect)
            @test collect(Map(identity), src; basesize=1) ==ₜ Base.collect(src)
            @test collect(identity(x) for x in src) ==ₜ Base.collect(src)
            @test collect((identity(x) for x in src); basesize=1) ==ₜ Base.collect(src)
        end
    end
end

end  # module


text:
eval(quote
    using .WriteVTK: VTKPolyhedron, vtk_grid

    export writevtk
    include("../Plotting/paraview.jl")
end)


text:
module Types

using Base.Random: UUID
using Base.Pkg.Types: VersionSet

export SHA1, VersionSpec, @vs_str

## ordering of UUIDs ##

Base.isless(a::UUID, b::UUID) = a.value < b.value

## SHA1 ##

struct SHA1
	bytes::Vector{UInt8}
	function SHA1(bytes::Vector{UInt8})
		length(bytes) == 20 ||
			throw(ArgumentError("wrong number of bytes for SHA1 hash: $(length(bytes))"))
		return new(bytes)
	end
end

Base.convert(::Type{SHA1}, s::String) = SHA1(hex2bytes(s))
Base.convert(::Type{Vector{UInt8}}, hash::SHA1) = hash.bytes
Base.convert(::Type{String}, hash::SHA1) = bytes2hex(Vector{UInt8}(hash))

Base.string(hash::SHA1) = String(hash)
Base.show(io::IO, hash::SHA1) = print(io, "SHA1(", String(hash), ")")
Base.isless(a::SHA1, b::SHA1) = lexless(a.bytes, b.bytes)
Base.hash(a::SHA1, h::UInt) = hash((SHA1, a.bytes), h)
Base.:(==)(a::SHA1, b::SHA1) = a.bytes == b.bytes

## VersionSpec ##

struct VersionSpec{n}
	spec::NTuple{n,Int}
	function VersionSpec{n}(spec::NTuple{n,Integer}) where n
		n <= 3 || throw(ArgumentError("VersionSpec: you can only specify major, minor and patch versions"))
		return new(spec)
	end
end
VersionSpec(spec::Integer...) = VersionSpec{length(spec)}(spec)

function Base.convert(::Type{VersionSpec}, s::AbstractString)
	isempty(s) || s == "*" ? VersionSpec() :
	VersionSpec(map(x->parse(Int, x), split(s, '.'))...)
end

macro vs_str(s::String); VersionSpec(s); end

function Base.show(io::IO, s::VersionSpec)
	print(io, "vs\"")
	join(io, s.spec, '.')
	print(io, '"')
end
Base.show(io::IO, ::VersionSpec{0}) = print(io, "vs\"*\"")
Base.getindex(s::VersionSpec, i::Int) = s.spec[i]

Base.in(v::VersionNumber, s::VersionSpec{0}) = true
Base.in(v::VersionNumber, s::VersionSpec{1}) = v.major == s.spec[1]
Base.in(v::VersionNumber, s::VersionSpec{2}) = v.major == s.spec[1] && v.minor == s.spec[2]
Base.in(v::VersionNumber, s::VersionSpec{3}) = v.major == s.spec[1] && v.minor == s.spec[2] && v.patch == s.spec[3]

Base.convert(::Type{VersionSet}, s::VersionSpec{0}) = VersionSet()
Base.convert(::Type{VersionSet}, s::VersionSpec{1}) =
	VersionSet(VersionNumber(s[1]), VersionNumber(s[1]+1))
Base.convert(::Type{VersionSet}, s::VersionSpec{2}) =
	VersionSet(VersionNumber(s[1], s[2]), VersionNumber(s[1], s[2]+1))
Base.convert(::Type{VersionSet}, s::VersionSpec{3}) =
	VersionSet(VersionNumber(s[1], s[2], s[3]), VersionNumber(s[1], s[2], s[3]+1))
Base.convert(::Type{VersionSet}, v::VersionNumber) =
	VersionSet(v, Base.nextpatch(v))

end # module


text:
const NOTEBOOKS = [
    ("初识 Julia",  "01-Introducing-Julia.jl"),
    ("概率入门",     "02-Basic-Probability.jl"),
    ("概率分布",     "03-Probability-Distributions.jl"),
    ("数据整理",     "04-Processing-and-Summarizing-Data.jl"),
    ("统计推论",     "05-Statistical-Inference-Concepts.jl"),
    ("置信区间",     "06-Confidence-Intervals.jl"),
    ("假设检验",     "07-Hypothesis-Testing.jl"),
    ("线性回归",     "08-Linear-Regression-and-Extensions.jl"),
    ("机器学习",     "09-Machine-Learning-Basics.jl"),
    ("动态模型",     "10-Simulation-of-Dynamic-Models.jl"),
    ("用 Julia 分析新冠病毒", "Covid-with-Julia.jl")
    ]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

Core.include(Main, "Base.jl")

using .Base

# Ensure this file is also tracked
pushfirst!(Base._included_files, (@__MODULE__, joinpath(@__DIR__, "Base.jl")))
pushfirst!(Base._included_files, (@__MODULE__, joinpath(@__DIR__, "sysimg.jl")))

# set up depot & load paths to be able to find stdlib packages
@eval Base creating_sysimg = true
Base.init_depot_path()
Base.init_load_path()

if Base.is_primary_base_module
# load some stdlib packages but don't put their names in Main
let
    # Stdlibs manually sorted in top down order
    stdlibs = [
            # No deps
            :Base64,
            :CRC32c,
            :SHA,
            :FileWatching,
            :Unicode,
            :Mmap,
            :Serialization,
            :Libdl,
            :Markdown,
            :LibGit2,
            :Logging,
            :Sockets,
            :Printf,
            :Profile,
            :Dates,
            :DelimitedFiles,
            :Random,
            :UUIDs,
            :Future,
            :LinearAlgebra,
            :SparseArrays,
            :SuiteSparse,
            :Distributed,
            :SharedArrays,
            :Pkg,
            :Test,
            :REPL,
            :Statistics,
        ]

    maxlen = maximum(textwidth.(string.(stdlibs)))

    print_time = (mod, t) -> (print(rpad(string(mod) * "  ", maxlen + 3, "─")); Base.time_print(t * 10^9); println())
    print_time(Base, (Base.end_base_include - Base.start_base_include) * 10^(-9))

    Base._track_dependencies[] = true
    Base.tot_time_stdlib[] = @elapsed for stdlib in stdlibs
        tt = @elapsed Base.require(Base, stdlib)
        print_time(stdlib, tt)
    end
    for dep in Base._require_dependencies
        dep[3] == 0.0 && continue
        push!(Base._included_files, dep[1:2])
    end
    empty!(Base._require_dependencies)
    Base._track_dependencies[] = false

    print_time("Stdlibs total", Base.tot_time_stdlib[])
end
end

# Clear global state
empty!(Core.ARGS)
empty!(Base.ARGS)
empty!(LOAD_PATH)
@eval Base creating_sysimg = false
Base.init_load_path() # want to be able to find external packages in userimg.jl

let
tot_time_userimg = @elapsed (Base.isfile("userimg.jl") && Base.include(Main, "userimg.jl"))

tot_time_base = (Base.end_base_include - Base.start_base_include) * 10.0^(-9)
tot_time = tot_time_base + Base.tot_time_stdlib[] + tot_time_userimg

println("Sysimage built. Summary:")
print("Total ─────── "); Base.time_print(tot_time               * 10^9); print(" \n");
print("Base: ─────── "); Base.time_print(tot_time_base          * 10^9); print(" "); show(IOContext(stdout, :compact=>true), (tot_time_base          / tot_time) * 100); println("%")
print("Stdlibs: ──── "); Base.time_print(Base.tot_time_stdlib[] * 10^9); print(" "); show(IOContext(stdout, :compact=>true), (Base.tot_time_stdlib[] / tot_time) * 100); println("%")
if isfile("userimg.jl")
print("Userimg: ──── "); Base.time_print(tot_time_userimg       * 10^9); print(" "); show(IOContext(stdout, :compact=>true), (tot_time_userimg       / tot_time) * 100); println("%")
end
end

empty!(LOAD_PATH)
empty!(DEPOT_PATH)

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
export MPO,
       randomMPO

struct MPO
  N_::Int
  A_::Vector{ITensor}

  MPO() = new(0,Vector{ITensor}())

  function MPO(N::Int, A::Vector{ITensor})
    new(N,A)
  end

  function MPO(sites)
    N = length(sites)
    v = Vector{ITensor}(undef, N)
    l = [Index(1, "Link,l=$ii") for ii ∈ 1:N-1]
    for ii in 1:N
      s = sites[ii]
      sp = prime(s)
      if ii == 1
        v[ii] = ITensor(s, sp, l[ii])
      elseif ii == N
        v[ii] = ITensor(l[ii-1], s, sp)
      else
        v[ii] = ITensor(l[ii-1], s, sp, l[ii])
      end
    end
    new(N,v)
  end

  function MPO(sites::SiteSet,
               ops::Vector{String})
    N = length(sites)
    its = Vector{ITensor}(undef, N)
    links = Vector{Index}(undef, N)
    for ii in 1:N
        si = sites[ii]
        spin_op = op(sites, ops[ii])
        links[ii] = Index(1, "Link,n=$ii")
        local this_it
        if ii == 1
            this_it = ITensor(links[ii], si, si')
            this_it[links[ii](1), s[:], s'[:]] = spin_op[si[:], si'[:]]
        elseif ii == N
            this_it = ITensor(links[ii-1], si, si')
            this_it[links[ii-1](1), si[:], si'[:]] = spin_op[si[:], si'[:]]
        else
            this_it = ITensor(links[ii-1], links[ii], si, si')
            this_it[links[ii-1](1), links[ii](1), si[:], si'[:]] = spin_op[si[:], si'[:]]
        end
        its[ii] = this_it
    end
    new(N,its)
  end

  function MPO(sites::SiteSet,
               ops::String)
    return MPO(sites, fill(ops, length(sites)))
  end

end

function randomMPO(sites,
                   m::Int=1)
  M = MPO(sites)
  for i=1:length(M)
    randn!(M[i])
    normalize!(M[i])
  end
  if m > 1
    error("randomMPS: currently only m==1 supported")
  end
  return M
end

length(m::MPO) = m.N_

getindex(m::MPO, n::Integer) = getindex(m.A_,n)
setindex!(m::MPO,T::ITensor,n::Integer) = setindex!(m.A_,T,n)

copy(m::MPO) = MPO(m.N_,copy(m.A_))

# TODO: optimize finding the index a little bit
# First do: scom = commonindex(A[j],x[j])
# Then do: uniqueindex(A[j],A[j-1],A[j+1],(scom,))
function siteindex(A::MPO,x::MPS,j::Integer)
  N = length(A)
  if j == 1
    si = uniqueindex(A[j],(A[j+1],x[j]))
  elseif j == N
    si = uniqueindex(A[j],(A[j-1],x[j]))
  else
    si = uniqueindex(A[j],(A[j-1],A[j+1],x[j]))
  end
  return si
end

function siteinds(A::MPO,x::MPS)
  N = length(A)
  is = IndexSet(N)
  for j in 1:N
    is[j] = siteindex(A,x,j)
  end
  return is
end

function prime!(M::T,vargs...) where {T <: Union{MPS,MPO}}
  N = length(M)
  for i ∈ 1:N
    prime!(M[i],vargs...)
  end
end

function primelinks!(M::T, plinc::Integer = 1) where {T <: Union{MPS,MPO}}
  N = length(M)
  for i ∈ 1:N-1
    l = linkindex(M,i)
    prime!(M[i],plinc,l)
    prime!(M[i+1],plinc,l)
  end
end

function simlinks!(M::T) where {T <: Union{MPS,MPO}}
  N = length(M)
  for i ∈ 1:N-1
    l = linkindex(M,i)
    l̃ = sim(l)
    #M[i] *= δ(l,l̃)
    replaceindex!(M[i],l,l̃)
    #M[i+1] *= δ(l,l̃)
    replaceindex!(M[i+1],l,l̃)
  end
end

function maxDim(M::T) where {T <: Union{MPS,MPO}}
  md = 1
  for b=1:length(M)-1
    md = max(md,dim(linkindex(M,b)))
  end
  return md
end

function show(io::IO,
              W::MPO)
  print(io,"MPO")
  (length(W) > 0) && print(io,"\n")
  for i=1:length(W)
    println(io,"$i  $(W[i])")
  end
end

"""
inner(y::MPS, A::MPO, x::MPS)

Compute <y|A|x>
"""
function inner(y::MPS,
               A::MPO,
               x::MPS)::Number
  N = length(A)
  if length(y) != N || length(x) != N
    error("inner: mismatched lengths $N and $(length(x)) or $(length(y))")
  end
  ydag = dag(y)
  simlinks!(ydag)
  sAx = siteinds(A,x)
  replacesites!(ydag,sAx)
  O = ydag[1]*A[1]*x[1]
  for j=2:N
    O = O*ydag[j]*A[j]*x[j]
  end
  return O[]
end


text:
module BattleSim

include("Stats.jl")
include("Ability.jl")
include("Pokemon.jl")

function __init__()

end

end

max_stars_repo_name: skeetcha/pta-battle-sim

max_issues_repo_name: skeetcha/pta-battle-sim

max_forks_repo_name: skeetcha/pta-battle-sim


text:
#=
parser2:
- Julia version:
- Author: bramb
- Date: 2019-02-28
=#
using Query,LibExpat,DataStructures,LightGraphs,MetaGraphs

@enum(TextState,TEXT,TAIL)

mutable struct Element
    name::String
    attributes::Dict{String,String}
    in_text_divergence::Bool
    parent::Element

    Element(name::String) = new(name,Dict{String,String}(),false)
end

mutable struct Triple
    element::Element
    text
    tail

    function Triple()
        element = Element("")
        text = IOBuffer()
        tail = IOBuffer();
        new(element,text,tail)
    end
end

mutable struct Context
    last_element_is_open::Bool
    triples::Array{Triple}
    text
    tail
    text_state::TextState
    open_elements::Array{Element}
    open_divergence_elements::Array{Element}

    function Context()
        last_element_is_open = false
        triples = []
        text = IOBuffer()
        tail = IOBuffer()
        text_state = TEXT;
        open_elements = []
        open_divergence_elements = []
        new(last_element_is_open,triples,text,tail,text_state,open_elements,open_divergence_elements)
    end
end

is_divergence_element(name::String) = name in ["subst", "choice", "app"]

function get_triples(xml::String)
    ctx = Context()

    cbs = XPCallbacks()
    cbs.start_element = function(h, name, attrs)
#         println("<$name>")
        triple = Triple()
        triple.element = Element(name)
        triple.element.in_text_divergence = !isempty(h.data.open_divergence_elements)
        push!(h.data.triples,triple)
        h.data.last_element_is_open = true
        h.data.text_state = TEXT
        if (!isempty(h.data.open_elements))
            parent = h.data.open_elements[1]
            triple.element.parent = parent
        end
        pushfirst!(h.data.open_elements,triple.element)
        if (is_divergence_element(triple.element.name))
            pushfirst!(h.data.open_divergence_elements,triple.element)
        end
    end
    cbs.end_element = function(h, name)
#         println("</$name>")
        h.data.text_state = TAIL
        popfirst!(h.data.open_elements)
        if (is_divergence_element(name))
            popfirst!(h.data.open_divergence_elements)
        end
    end
    cbs.character_data = function(h, txt)
#         println("\"$txt\"")
        # now we add to either text or tail, depending on text_state
        if (h.data.text_state == TEXT)
            print(h.data.triples[end].text,txt)
        else
            print(h.data.triples[end].tail,txt)
        end
    end

    parse(xml,cbs,data = ctx)
    return ctx.triples
end

parent(x::Triple) = isdefined(x.element, :parent) ? x.element.parent : "XML"

group_triples(triples) = triples |> @groupby(parent(_)) |> collect

serialize_text(t::Triple) = String(take!(t.text))

serialize_tail(t::Triple) = String(take!(t.tail))

function serialize_group(g)
    buf = IOBuffer()
    if (length(g) > 1)
        if (g[1].element.in_text_divergence)
            print(buf,"<|")
            for t in g
                print(buf,serialize_text(t),"|")
            end
            print(buf,">",serialize_tail(g[end]))
        else
            for t in g
                print(buf,serialize_text(t),serialize_tail(t))
            end
        end
    else
        print(buf,serialize_text(g[1]),serialize_tail(g[1]))
    end
    return String(take!(buf))
end

function serialize_grouped_triples(grouped_triples)
    serbuf = IOBuffer()
    for group in grouped_triples
        print(serbuf,serialize_group(group))
    end
    return String(take!(serbuf))
end

@enum(VertexType,TEXTNODE,DIVERGENCE,CONVERGENCE)

# function to_graph(grouped_triples::Array{Group{Triple},1})
#     mg = MetaGraph(SimpleGraph())
#     for group in grouped_triples
#         for triple in group
#             text = serialize_text(triple)
#             @show(text)
#             if (!isempty(text))
#                 add_vertices!(mg.graph,1)
#                 v = nv(mg.graph)
#                 set_props!(mg,v,Dict(:type => TEXTNODE,:text => text))
#                 v>1 && add_edge!(mg.graph, v-1, v)
#             end
#             tail = serialize_tail(triple)
#             @show(tail)
#             if (!isempty(tail))
#                 add_vertices!(mg.graph,1)
#                 v = nv(mg.graph)
#                 set_props!(mg,v,Dict(:type => TEXTNODE,:text => tail))
#                 v>1 && add_edge!(mg.graph, v-1, v)
#             end
#         end
#     end
#     return mg
# end


text:
type UBEV<:FiniteHorizonAgent
    policy::Matrix{UInt64}
    ns::Array{UInt64, 4}
    phat::Array{Float64, 4}
    n::Array{UInt64, 3}
    rewards::Array{Float64, 3}
    Vopt::Matrix{Float64}
    Qopt::Array{Float64, 3}
    δ::Float64
    rewards_known::Bool
    maxRet::Float64
    maxR::Float64
    explore_bonus::Float64
end

function UBEV(S, A, H, rewards::Array, δ, explore_bonus=1.)
    policy = zeros(UInt64, S, H)
    ns = zeros(UInt64, S, S, A, H)
    phat = zeros(S, S, A, H)
    phat[1, :, :, :] = 1.
    n = zeros(UInt64, S, A, H)
    maxRet = sum(maximum(reshape(rewards, S*A, H), 1))
    m = UBEV(policy, ns, phat, n, rewards, zeros(S, H+1), zeros(S, A, H+1), δ, true, maxRet, maximum(rewards[:]), explore_bonus)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end

function UBEV(S, A, H, δ::Real, explore_bonus=1.)
    policy = zeros(UInt64, S, H)
    ns = zeros(UInt64, S, S, A, H)
    phat = zeros(S, S, A, H)
    phat[1, :, :, :] = 1.
    n = zeros(UInt64, S, A, H)
    maxr = 1.
    maxRet = H * maxr
    rewards = zeros(S, A, H)
    m = UBEV(policy, ns, phat, n, rewards, zeros(S, H+1), zeros(S, A, H+1), δ, false, maxRet, maxr, explore_bonus)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end
maxV(m::UBEV) = m.maxRet
maxR(m::UBEV) = m.maxR

llnp(x) = max(0, log(max(0, log(max(0, x)))))
phi(n, S, A, H, δprime, range) = range * √(1./n * (2*llnp(n) + log(3 * S * A * H / δprime)))
confint(n, lC, valmax) =  valmax * √(1./n * (2*llnp(n) + lC))

function update_policy!(m::UBEV)
    S = nS(m)
    H = horizon(m)
    A = nA(m)
    δprime = m.δ
    if m.rewards_known
        δprime /= 7
    else
        δprime /= 9
    end
    Vmax = maxV(m)
    Rmax = maxR(m)
    Q = zeros(A)
    lC = log(3 * S * A * H / δprime)
    m.phat .= m.ns ./ reshape(m.n, 1, S, A, H)
    for t=H:-1:1
        V = m.Vopt[:, t+1]
        curmaxV = min(maximum(V), Vmax)
        curminV = min(minimum(V), Vmax)
        for s ∈ 1:S
            for a ∈ 1:A

                EV = curmaxV
                if m.n[s,a,t] > 0
                    EV = min(curmaxV, m.phat[:, s, a, t] ⋅ V + confint(m.n[s,a,t], lC, curmaxV - curminV) * m.explore_bonus)
                end

                r = min(m.rewards[s, a, t], Rmax)
                if !m.rewards_known
                    r = Rmax
                    if m.n[s,a,t] > 0
                        r = min(m.rewards[s, a, t] / m.n[s, a, t] + confint(m.n[s, a, t], lC, Rmax), Rmax)
                    end
                end

                Q[a] = r + EV
            end
            m.Qopt[s,:,t] = Q
            bestA = indmax(Q)
            m.policy[s, t] = bestA
            m.Vopt[s, t] = Q[bestA]
        end
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
ACASX grammar problem. Three runtypes are currently supported.\n
'nvn' or :nmacs_vs_nonnmacs (1 is nmac, 2 is non-nmac).\n
'nc' or :nmac_clusters (cluster labels provided as input, i.e., by manuals)\n
'nec' or :nonnmacs_extra_cluster.  same as 'nc' except nonnmacs are added as an
extra cluster to the end.  Default: 'nvn'
Two versions of the fitness function are provided, one with pruning (early stop) and one without
(always evaluates over entire dataset)
"""
module ACASXProblem

export ACASXClustering, create_grammar, get_fitness, eval_expr, apply_expr,
    get_members, get_grammar
export FMT_PRETTY, FMT_NATURAL
export entropy_metrics, gini_metrics

import Compat.view

using Datasets
using ExprSearch
using Devectorize
using RLESUtils, LogicUtils, Interpreter, TreeUtils
import ExprSearch: ExprProblem, get_fitness, get_grammar
import GBDTs.apply_expr

include("labeleddata.jl")
include("infogain.jl")
include("format.jl")

#weights for fitness function
const W_METRIC = 100.0 #entropy
const W_NUM_NODES = 0.1 #length of expression

typealias RealVec Union{DataArray{Float64,1}, Vector{Float64}}

type ACASXClustering{T} <: ExprProblem
  Dl::DFSetLabeled{T}
  w_metric::Float64
  w_num_nodes::Float64
  labelset::Vector{T}
  grammar::Grammar
end

"""
Example: problem = ACASXClustering(:nmacs_vs_nonnmacs, "dasc", "", "")
"""
function ACASXClustering(runtype::Symbol, dataname::AbstractString, manuals::AbstractString,
    clustername::AbstractString, w_metric::Float64=W_METRIC, w_num_nodes::Float64=W_NUM_NODES)
    if runtype == :nmacs_vs_nonnmacs
        out = ACASXClustering(dataname, w_metric, w_num_nodes)
    elseif runtype == :nmac_clusters
        out = ACASXClustering(dataname, manuals, clustername, w_metric, w_num_nodes;
            incl_nonnmacs=false)
    elseif runtype == :nonnmacs_extra_cluster
        out = ACASXClustering(dataname, manuals, clustername, w_metric, w_num_nodes;
            incl_nonnmacs=true)
    else
        error("Runtype not defined ($runtype)")
    end
    out
end

#nmacs vs non-nmacs
function ACASXClustering(dataname::AbstractString,
                         w_metric::Float64=W_METRIC, w_num_nodes::Float64=W_NUM_NODES)
    data = dataset(dataname)
    Dl = nmacs_vs_nonnmacs(data)
    labelset = unique(labels(Dl))
    grammar = create_grammar()
    return ACASXClustering(Dl, w_metric, w_num_nodes, labelset, grammar)
end

#with manual clustering
function ACASXClustering(dataname::AbstractString,
                         manuals::AbstractString,
                         clustername::AbstractString,
                         w_metric::Float64=W_METRIC, w_num_nodes::Float64=W_NUM_NODES;
                         incl_nonnmacs::Bool=true)

    data = dataset(dataname)
    clustering = dataset(manuals, clustername)

    if incl_nonnmacs
        Dl = nonnmacs_extra_cluster(clustering, data)
    else
        Dl = nmac_clusters(clustering, data)
    end
    labelset = unique(labels(Dl))
    return ACASXClustering(Dl, w_metric, w_num_nodes, labelset, grammar)
end

ExprSearch.get_grammar{T}(problem::ACASXClustering{T}) = problem.grammar
function ExprSearch.get_fitness{T}(problem::ACASXClustering{T}, derivtree::DerivationTree,
    userargs::SymbolTable)
    ids = userargs[:ids]

    expr = get_expr(derivtree)
    induced_labels = apply_expr(problem, ids, expr)
    num_nodes = count_nodes(derivtree.root)

    #_, _, metric = entropy_metrics(predicts, Dl.labels, Float64(problem.nlabels))
    _, _, metric = gini_metrics(induced_labels, view(problem.Dl.labels, ids))
    return problem.w_metric * metric + problem.w_num_nodes * num_nodes
end

function create_grammar()
    @grammar grammar begin
        start = bin

        bin = always | eventually | implies | count
        always = Expr(:call, :G, bin_vec) #global
        eventually = Expr(:call, :F, bin_vec) #future
        implies = Expr(:call, :Y, bin_vec, bin_vec)
        count = Expr(:call, :ctlt, bin_vec, timestep) | Expr(:call, :ctle, bin_vec, timestep) | Expr(:call, :ctgt, bin_vec, timestep) | Expr(:call, :ctge, bin_vec, timestep) | Expr(:call, :cteq, bin_vec, timestep)

        #produces a bin_vec
        bin_vec = bin_feat | and | or | not  | eq | lt | lte | abseq | abslt | abslte | diff_eq | diff_lt | diff_lte | sign | absdiff_eq | absdiff_lt | absdiff_lte
        and = Expr(:call, :&, bin_vec, bin_vec)
        or = Expr(:call, :|, bin_vec, bin_vec)
        not = Expr(:call, :!, bin_vec)

    #equal
        eq = vrate_eq | altdiff_eq | chi_angle_eq | psi_angle_eq | sr_eq | timer_eq | psid_eq | v_eq | alt_eq | abs_altdiff_eq
        vrate_eq = Expr(:call, :eq, vrate_feat, vrate_val) | Expr(:call, :eq, vrate_feat, vrate_feat)
        altdiff_eq = Expr(:call, :eq, altdiff_feat, altdiff_val)
        chi_angle_eq = Expr(:call, :eq, chi_angle_feat, angle_val) | Expr(:call, :eq, chi_angle_feat, chi_angle_feat)
        psi_angle_eq = Expr(:call, :eq, psi_angle_feat, psi_angle_feat)
        sr_eq = Expr(:call, :eq, sr_feat, sr_val)
        timer_eq = Expr(:call, :eq, timer_feat, timer_val) | Expr(:call, :eq, timer_feat, timer_feat)
        psid_eq = Expr(:call, :eq, psid_feat, psid_val) | Expr(:call, :eq, psid_feat, psid_feat)
        v_eq = Expr(:call, :eq, v_feat, v_val) | Expr(:call, :eq, v_feat, v_feat)
        alt_eq = Expr(:call, :eq, alt_feat, alt_val) | Expr(:call, :eq, alt_feat, alt_feat)
        abs_altdiff_eq = Expr(:call, :eq, abs_altdiff_feat, abs_altdiff_val)

        #absolute equals
        abseq = vrate_abseq | altdiff_abseq | chi_angle_abseq | psi_angle_abseq | sr_abseq | timer_abseq | psid_abseq | v_abseq | alt_abseq
        vrate_abseq = Expr(:call, :abeq, vrate_feat, vrate_val) | Expr(:call, :abeq, vrate_feat, vrate_feat)
        altdiff_abseq = Expr(:call, :abeq, altdiff_feat, altdiff_val)
        chi_angle_abseq = Expr(:call, :abeq, chi_angle_feat, angle_val) | Expr(:call, :abeq, chi_angle_feat, chi_angle_feat)
        psi_angle_abseq = Expr(:call, :abeq, psi_angle_feat, psi_angle_feat)
        sr_abseq = Expr(:call, :abeq, sr_feat, sr_val)
        timer_abseq = Expr(:call, :abeq, timer_feat, timer_val) | Expr(:call, :abeq, timer_feat, timer_feat)
        psid_abseq = Expr(:call, :abeq, psid_feat, psid_val) | Expr(:call, :abeq, psid_feat, psid_feat)
        v_abseq = Expr(:call, :abeq, v_feat, v_val) | Expr(:call, :abeq, v_feat, v_feat)
        alt_abseq = Expr(:call, :abeq, alt_feat, alt_val) | Expr(:call, :abeq, alt_feat, alt_feat)

        #less than
        lt = vrate_lt | altdiff_lt  | chi_angle_lt | psi_angle_lt | sr_lt | timer_lt | psid_lt | v_lt | alt_lt | abs_altdiff_lt
        vrate_lt = Expr(:call, :lt, vrate_feat, vrate_val) | Expr(:call, :lt, vrate_feat, vrate_feat)
        altdiff_lt = Expr(:call, :lt, altdiff_feat, altdiff_val)
        chi_angle_lt = Expr(:call, :lt, chi_angle_feat, angle_val) | Expr(:call, :lt, chi_angle_feat, chi_angle_feat)
        psi_angle_lt = Expr(:call, :lt, psi_angle_feat, psi_angle_feat)
        sr_lt = Expr(:call, :lt, sr_feat, sr_val)
        timer_lt = Expr(:call, :lt, timer_feat, timer_val) | Expr(:call, :lt, timer_feat, timer_feat)
        psid_lt = Expr(:call, :lt, psid_feat, psid_val) | Expr(:call, :lt, psid_feat, psid_feat)
        v_lt = Expr(:call, :lt, v_feat, v_val) | Expr(:call, :lt, v_feat, v_feat)
        alt_lt = Expr(:call, :lt, alt_feat, alt_val) | Expr(:call, :lt, alt_feat, alt_feat)
        abs_altdiff_lt = Expr(:call, :lt, abs_altdiff_feat, abs_altdiff_val)

        #absolute less than
        abslt = vrate_abslt | altdiff_abslt  | chi_angle_abslt | psi_angle_abslt | sr_abslt | timer_abslt | psid_abslt | v_abslt | alt_abslt
        vrate_abslt = Expr(:call, :ablt, vrate_feat, vrate_val) | Expr(:call, :ablt, vrate_feat, vrate_feat)
        altdiff_abslt = Expr(:call, :ablt, altdiff_feat, altdiff_val)
        chi_angle_abslt = Expr(:call, :ablt, chi_angle_feat, angle_val) | Expr(:call, :ablt, chi_angle_feat, chi_angle_feat)
        psi_angle_abslt = Expr(:call, :ablt, psi_angle_feat, psi_angle_feat)
        sr_abslt = Expr(:call, :ablt, sr_feat, sr_val)
        timer_abslt = Expr(:call, :ablt, timer_feat, timer_val) | Expr(:call, :ablt, timer_feat, timer_feat)
        psid_abslt = Expr(:call, :ablt, psid_feat, psid_val) | Expr(:call, :ablt, psid_feat, psid_feat)
        v_abslt = Expr(:call, :ablt, v_feat, v_val) | Expr(:call, :ablt, v_feat, v_feat)
        alt_abslt = Expr(:call, :ablt, alt_feat, alt_val) | Expr(:call, :ablt, alt_feat, alt_feat)

        #less then or equal
        lte = vrate_lte | altdiff_lte | chi_angle_lte | psi_angle_lte | sr_lte | timer_lte | psid_lte | v_lte | alt_lte | abs_altdiff_lte
        vrate_lte = Expr(:call, :lte, vrate_feat, vrate_val) | Expr(:call, :lte, vrate_feat, vrate_feat)
        altdiff_lte = Expr(:call, :lte, altdiff_feat, altdiff_val)
        chi_angle_lte = Expr(:call, :lte, chi_angle_feat, angle_val) | Expr(:call, :lte, chi_angle_feat, chi_angle_feat)
        psi_angle_lte = Expr(:call, :lte, psi_angle_feat, psi_angle_feat)
        sr_lte = Expr(:call, :lte, sr_feat, sr_val)
        timer_lte = Expr(:call, :lte, timer_feat, timer_val) | Expr(:call, :lte, timer_feat, timer_feat)
        psid_lte = Expr(:call, :lte, psid_feat, psid_val) | Expr(:call, :lte, psid_feat, psid_feat)
        v_lte = Expr(:call, :lte, v_feat, v_val) | Expr(:call, :lte, v_feat, v_feat)
        alt_lte = Expr(:call, :lte, alt_feat, alt_val) | Expr(:call, :lte, alt_feat, alt_feat)
        abs_altdiff_lte = Expr(:call, :lte, abs_altdiff_feat, abs_altdiff_val)

        #abs less then or equal
        abslte = vrate_abslte | altdiff_abslte | chi_angle_abslte | psi_angle_abslte | sr_abslte | timer_abslte | psid_abslte | v_abslte | alt_abslte
        vrate_abslte = Expr(:call, :able, vrate_feat, vrate_val) | Expr(:call, :able, vrate_feat, vrate_feat)
        altdiff_abslte = Expr(:call, :able, altdiff_feat, altdiff_val)
        chi_angle_abslte = Expr(:call, :able, chi_angle_feat, angle_val) | Expr(:call, :able, chi_angle_feat, chi_angle_feat)
        psi_angle_abslte = Expr(:call, :able, psi_angle_feat, psi_angle_feat)
        sr_abslte = Expr(:call, :able, sr_feat, sr_val)
        timer_abslte = Expr(:call, :able, timer_feat, timer_val) | Expr(:call, :able, timer_feat, timer_feat)
        psid_abslte = Expr(:call, :able, psid_feat, psid_val) | Expr(:call, :able, psid_feat, psid_feat)
        v_abslte = Expr(:call, :able, v_feat, v_val) | Expr(:call, :able, v_feat, v_feat)
        alt_abslte = Expr(:call, :able, alt_feat, alt_val) | Expr(:call, :able, alt_feat, alt_feat)

        #sign
        sign = vrate_sign | chi_angle_sign | psid_sign
        vrate_sign = Expr(:call, :sn, vrate_feat, vrate_feat)
        chi_angle_sign = Expr(:call, :sn, chi_angle_feat, chi_angle_feat)
        psid_sign = Expr(:call, :sn, psid_feat, psid_feat)

        #difference is equal
        diff_eq = vrate_diff_eq | chi_angle_diff_eq | psi_angle_diff_eq | timer_diff_eq | psid_diff_eq | v_diff_eq
        vrate_diff_eq = Expr(:call, :dfeq, vrate_feat, vrate_feat, vrate_val)
        chi_angle_diff_eq = Expr(:call, :dfeq, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_diff_eq = Expr(:call, :dfeq, psi_angle_feat, psi_angle_feat, angle_val)
        timer_diff_eq = Expr(:call, :dfeq, timer_feat, timer_feat, timer_val)
        psid_diff_eq = Expr(:call, :dfeq, psid_feat, psid_feat, psid_val)
        v_diff_eq = Expr(:call, :dfeq, v_feat, v_feat, v_val)

        #abs difference is equal
        absdiff_eq = vrate_absdiff_eq | chi_angle_absdiff_eq | psi_angle_absdiff_eq | timer_absdiff_eq | psid_absdiff_eq | v_absdiff_eq
        vrate_absdiff_eq = Expr(:call, :adfeq, vrate_feat, vrate_feat, vrate_val)
        chi_angle_absdiff_eq = Expr(:call, :adfeq, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_absdiff_eq = Expr(:call, :adfeq, psi_angle_feat, psi_angle_feat, angle_val)
        timer_absdiff_eq = Expr(:call, :adfeq, timer_feat, timer_feat, timer_val)
        psid_absdiff_eq = Expr(:call, :adfeq, psid_feat, psid_feat, psid_val)
        v_absdiff_eq = Expr(:call, :adfeq, v_feat, v_feat, v_val)

        #difference is less than some value.  Values may not be appropriate after taking difference
        diff_lt = vrate_diff_lt | chi_angle_diff_lt | psi_angle_diff_lt | timer_diff_lt | psid_diff_lt | v_diff_lt
        vrate_diff_lt = Expr(:call, :dflt, vrate_feat, vrate_feat, vrate_val)
        chi_angle_diff_lt = Expr(:call, :dflt, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_diff_lt = Expr(:call, :dflt, psi_angle_feat, psi_angle_feat, angle_val)
        timer_diff_lt = Expr(:call, :dflt, timer_feat, timer_feat, timer_val)
        psid_diff_lt = Expr(:call, :dflt, psid_feat, psid_feat, psid_val)
        v_diff_lt = Expr(:call, :dflt, v_feat, v_feat, v_val)

        #abs difference is less than some value.  Values may not be appropriate after taking difference
        absdiff_lt = vrate_absdiff_lt | chi_angle_absdiff_lt | psi_angle_absdiff_lt | timer_absdiff_lt | psid_absdiff_lt | v_absdiff_lt
        vrate_absdiff_lt = Expr(:call, :adflt, vrate_feat, vrate_feat, vrate_val)
        chi_angle_absdiff_lt = Expr(:call, :adflt, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_absdiff_lt = Expr(:call, :adflt, psi_angle_feat, psi_angle_feat, angle_val)
        timer_absdiff_lt = Expr(:call, :adflt, timer_feat, timer_feat, timer_val)
        psid_absdiff_lt = Expr(:call, :adflt, psid_feat, psid_feat, psid_val)
        v_absdiff_lt = Expr(:call, :adflt, v_feat, v_feat, v_val)

        #difference is less than or equal to some value.  Values may not be appropriate after taking difference
        diff_lte = vrate_diff_lte | chi_angle_diff_lte | psi_angle_diff_lte | timer_diff_lte | psid_diff_lte | v_diff_lte
        vrate_diff_lte = Expr(:call, :dfle, vrate_feat, vrate_feat, vrate_val)
        chi_angle_diff_lte = Expr(:call, :dfle, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_diff_lte = Expr(:call, :dfle, psi_angle_feat, psi_angle_feat, angle_val)
        timer_diff_lte = Expr(:call, :dfle, timer_feat, timer_feat, timer_val)
        psid_diff_lte = Expr(:call, :dfle, psid_feat, psid_feat, psid_val)
        v_diff_lte = Expr(:call, :dfle, v_feat, v_feat, v_val)

        #absolute difference is less than or equal to some value.
        absdiff_lte = vrate_absdiff_lte | chi_angle_absdiff_lte | psi_angle_absdiff_lte | timer_absdiff_lte | psid_absdiff_lte | v_absdiff_lte
        vrate_absdiff_lte = Expr(:call, :adfle, vrate_feat, vrate_feat, vrate_val)
        chi_angle_absdiff_lte = Expr(:call, :adfle, chi_angle_feat, chi_angle_feat, angle_val)
        psi_angle_absdiff_lte = Expr(:call, :adfle, psi_angle_feat, psi_angle_feat, angle_val)
        timer_absdiff_lte = Expr(:call, :adfle, timer_feat, timer_feat, timer_val)
        psid_absdiff_lte = Expr(:call, :adfle, psid_feat, psid_feat, psid_val)
        v_absdiff_lte = Expr(:call, :adfle, v_feat, v_feat, v_val)

        #read features
        bin_feat = Expr(:call, :g, :D, bin_feat_id)
        vrate_feat = Expr(:call, :g, :D, vrate_feat_id)
        altdiff_feat = Expr(:call, :g, :D, altdiff_feat_id)
        abs_altdiff_feat = Expr(:call, :g, :D, abs_altdiff_feat_id)
        angle_feat = Expr(:call, :g, :D, angle_feat_id)
        psi_angle_feat = Expr(:call, :g, :D, psi_angle_feat_id)
        chi_angle_feat = Expr(:call, :g, :D, chi_angle_feat_id)
        sr_feat = Expr(:call, :g, :D, sr_feat_id)
        timer_feat = Expr(:call, :g, :D, timer_feat_id)
        psid_feat = Expr(:call, :g, :D, psid_feat_id)
        v_feat = Expr(:call, :g, :D, v_feat_id)
        alt_feat = Expr(:call, :g, :D, alt_feat_id)

        #indices of each type
        #include bit flags
        #bin_feat_id = 1 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 23 | 24 | 25 | 26 | 27 | 28 | 30 | 31 | 32 | 38 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 60 | 61 | 62 | 63 | 64 | 65 | 67 | 68 | 69 | 75
        #no bit flags
        bin_feat_id = 1 | 23 | 24 | 25 | 30 | 31 | 32 | 38 | 60 | 61 | 62 | 67 | 68 | 69 | 75

        vrate_feat_id = 2 | 22 | 34 | 39 | 59 | 71
        altdiff_feat_id = 3 | 40
        abs_altdiff_feat_id = 76
        angle_feat_id = psi_angle_feat_id | chi_angle_feat_id
        psi_angle_feat_id = 4 | 41
        chi_angle_feat_id = 6 | 43
        sr_feat_id = 5 | 42 | 77
        timer_feat_id = 33 | 70
        psid_feat_id = 35 | 72
        v_feat_id = 36 | 73
        alt_feat_id = 37 | 74

        #values of each type
        vrate_val = -50 | -40 | -30 | -20 | -10 | -1 | 0 | 1 | 10 | 20 | 30 | 40 | 50
        altdiff_val = -2000 | -1500 | -1000 | -500 | -250 | -100 | -50 | -25 | -10 | -5 | -1 | 0 | 1 | 5 | 10 | 25 | 50 | 100 | 250 | 500 | 1000 | 1500 | 2000
        abs_altdiff_val = 0 | 1 | 5 | 10 | 25 | 50 | 100 | 250 | 500 | 1000 | 1500 | 2000
        angle_val = -180 | -135 | -90 | -45 | 0 | 45 | 90 | 135 | 180
        sr_val = 30000 | 25000 | 20000 | 15000 | 10000 | 7500 | 5000 | 2500 | 1000 | 500 | 250 | 100 | 50 | 25 | 10 | 1 | 0
        timer_val = 0:5
        psid_val = -10:10
        v_val = 50 | 75 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 450 | 500
        alt_val = 1000 | 2000 | 3000 | 5000 | 10000 | 15000 | 18000 | 20000 | 25000 | 29000
        timestep = 0:50
    end

    #automatically determine real vs bool columns from DataFrame
    #bin_ids, _ = feat_type_ids(D)
    #bin_terms = map(GrammaticalEvolution.Terminal, bin_ids)
    #grammar.rules[:bin_feat_id] = OrRule("bin_feat_id", bin_terms, nothing)

    grammar
end

function eq(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] == v2[i]
    end
    r
end
function eq(v1::RealVec, b::Real)
    v1 = convert(Array, v1)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] == b
    end
    r
end
function lt(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] < v2[i]
    end
    r
end
function lt(v1::RealVec, b::Real)
    v1 = convert(Array, v1)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] < b
    end
    r
end
function lte(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] <= v2[i]
    end
    r
end

function lte(v1::RealVec, b::Real)
    v1 = convert(Array, v1)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] <= b
    end
    r
end

function abs_eq(v::RealVec, b::Real)
    v = convert(Array, v)
    r = Array(Bool, length(v))
    @inbounds for i = 1:length(v)
        r[i] = abs(v[i]) == b
    end
    r
end
function abs_lte(v::RealVec, b::Real)
    v = convert(Array, v)
    r = Array(Bool, length(v))
    @inbounds for i = 1:length(v)
        r[i] = abs(v[i]) <= b
    end
    r
end
function abs_lt(v::RealVec, b::Real)
    v = convert(Array, v)
    r = Array(Bool, length(v))
    @inbounds for i = 1:length(v)
        r[i] = abs(v[i]) < b
    end
    r
end

function abs_eq(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i]) == abs(v2[i])
    end
    r
end
function abs_lte(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i]) <= abs(v2[i])
    end
    r
end
function abs_lt(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i]) < abs(v2[i])
    end
    r
end

function diff_eq(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = (v1[i] - v2[i]) == b
    end
    r
end
function diff_lte(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = (v1[i] - v2[i]) <= b
    end
    r
end
function diff_lt(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = (v1[i] - v2[i]) < b
    end
    r
end

function abs_diff_eq(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i] - v2[i]) == b
    end
    r
end
function abs_diff_lte(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i] - v2[i]) <= b
    end
    r
end
function abs_diff_lt(v1::RealVec, v2::RealVec, b::Real)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = abs(v1[i] - v2[i]) < b
    end
    r
end

function eventually(v::AbstractVector{Bool})
    v = convert(Array, v)
    any(v)
end
function globally(v::AbstractVector{Bool})
    v = convert(Array, v)
    all(v)
end
function implies(v1::AbstractVector{Bool}, v2::AbstractVector{Bool})
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v2[i] || !v1[i]
    end
    all(r)
end

function sign_(v1::RealVec, v2::RealVec)
    v1 = convert(Array, v1)
    v2 = convert(Array, v2)
    r = Array(Bool, length(v1))
    @inbounds for i = 1:length(v1)
        r[i] = v1[i] * v2[i] >= 0.0
    end
    r
end

function count_eq(v::AbstractVector{Bool}, b::Real)
    v = convert(Array, v)
    count(identity, v) == b
end
function count_lt(v::AbstractVector{Bool}, b::Real)
    v = convert(Array, v)
    count(identity, v) < b
end
function count_lte(v::AbstractVector{Bool}, b::Real)
    v = convert(Array, v)
    count(identity, v) <= b
end
function count_gt(v::AbstractVector{Bool}, b::Real)
    v = convert(Array, v)
    count(identity, v) > b
end
function count_gte(v::AbstractVector{Bool}, b::Real)
    v = convert(Array, v)
    count(identity, v) >= b
end

#get
get_ref(D, id) = D[:, id]

#shorthands used in grammar to reduce impact on code length
#g = get_ref
#abeq = abs_eq
#able = abs_lte
#ablt = abs_lt
#dfeq = diff_eq
#dfle = diff_lte
#dflt = diff_lt
#adfeq = abs_diff_eq
#adfle = abs_diff_lte
#adflt = abs_diff_lt
#F = eventually
#G = globally
#Y = implies
#sn = sign_ #avoid conflict with Base.sign
#ctlt = count_lt
#ctle = count_lte
#ctgt = count_gt
#ctge = count_gte
#cteq = count_eq

const SYMTABLE = SymbolTable(
    :g => get_ref,
    :abeq => abs_eq,
    :able => abs_lte,
    :ablt => abs_lt,
    :dfeq => diff_eq,
    :dfle => diff_lte,
    :dflt => diff_lt,
    :adfeq => abs_diff_eq,
    :adfle => abs_diff_lte,
    :adflt => abs_diff_lt,
    :F => eventually,
    :G => globally,
    :Y => implies,
    :sn => sign_, #avoid conflict with Base.sign
    :ctlt => count_lt,
    :ctle => count_lte,
    :ctgt => count_gt,
    :ctge => count_gte,
    :cteq => count_eq,
    :eq => eq,
    :lt => lt,
    :lte => lte,
    :| => |,
    :& => &,
    :! => !
    )

function eval_expr(problem::ACASXClustering, expr, D)
    SYMTABLE[:D] = D
    return interpret(SYMTABLE, expr)
end

function apply_expr{T}(problem::ACASXClustering{T}, ids::Vector{Int64}, expr)
    expr_labels = map(i->eval_expr(problem, expr, getrecords(problem.Dl, i)), ids)
    expr_labels::AbstractVector{Bool}
end

function get_members{T}(problem::ACASXClustering{T}, ids::Vector{Int64},
    induced_labels::Vector{Bool})
    meta = getmeta(problem.Dl)
    encounter_ids = meta[:encounter_id]
    true_ids = encounter_ids[ids[find(induced_labels)]]
    false_ids = encounter_ids[ids[find(!induced_labels)]]
    (true_ids, false_ids)
end

function get_members{T}(problem::ACASXClustering{T}, ids::Vector{Int64}, expr)
    expr_labels = apply_expr(problem, ids, expr)
    members = get_members(problem, ids, expr_labels)
    members
end

type CountTracker{T}
    c_true::Dict{T,Int64} #counts for the true cluster
    c_false::Dict{T,Int64} #counts for the false cluster
    N_true::Int64 #number of items in true cluster
    N_false::Int64 #number of items in false cluster
end

function CountTracker{T}(problem::ACASXClustering{T})
    c_true = Dict{T,Int64}()
    c_false = Dict{T,Int64}()
    N_true = N_false = 0
    for l in problem.labelset
        c_true[l] = 0
        c_false[l] = 0
    end
    CountTracker(c_true, c_false, 0, 0)
end

function increment!{T}(tracker::CountTracker{T}, ::Type{Val{true}}, label::T)
    tracker.c_true[label] += 1
    tracker.N_true += 1
end

function increment!{T}(tracker::CountTracker{T}, ::Type{Val{false}}, label::T)
    tracker.c_false[label] += 1
    tracker.N_false += 1
end

function gini_optimistic{T}(tracker::CountTracker{T}, N_remaining::Int64)
    c_true = collect(values(tracker.c_true))
    c_false = collect(values(tracker.c_false))

    #add to true
    i = indmax(c_true)
    c_true[i] += N_remaining
    g_true = gini_from_counts(c_true, c_false)

    #add to false
    c_true = collect(values(tracker.c_true)) #restore
    i = indmax(c_false) #mode
    c_false[i] += N_remaining
    g_false = gini_from_counts(c_true, c_false)

    min(g_true, g_false)
end

end #module

#Code that might be useful later
#=
function get_fitness{T}(problem::ACASXClustering{T}, expr,
                                   thresh::Float64, default::Float64,
                                   earlystop_div::Int64)

  Dl = problem.Dl
  codelen = length(string(expr))
  f = to_function(problem, expr)

  metric_thresh = (thresh - problem.w_len * codelen) / problem.w_metric #translate thresh to bound on metric
  c_tracker = CountTracker(problem)

  earlystop_int = round(Int64, length(Dl) / earlystop_div)

  predicts = Array(Bool, length(Dl))
  for i = 1:length(Dl)
    predicts[i] = f(getrecords(Dl, i))

    increment!(c_tracker, Val{predicts[i]}, Dl.labels[i])

    #evaluate early exit
    if rem(i, earlystop_int) == 0 #on interval
      optim = gini_optimistic(c_tracker, length(Dl) - i)
      if optim > metric_thresh #if most optimistic case still doesn't meet thresh, early exit
        #@show metric_thresh
        #@show optim
        #@show i
        return default
      end
    end
  end

  #do full calc
  _, _, metric = gini_metrics(predicts, Dl.labels)
  return problem.w_metric * metric + problem.w_len * codelen
end
=#


text:
# In our second example, we will compare the `Arctan` function
# to various partial sums of its power series. To do this,
# we will make a slider to compare the phase portraits,
# and allow the user to interactively choose the size of the sum.
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("sans-serif")

# Build our slider
function get_param_n()
    value_min, value_max = 1, 19
    # Our slider position only moves in integer-length steps
    slider_pos = CV_TranslateByOffset(Int)
    # Our starting slider position is 4
    slider_pos.value = 4
    set_slider_value_func = z -> begin
        slider_pos.value = round(UInt, max(min(z, value_max), value_min))
        return CV_Response(; redraw_flag=true)
    end
    # The next 4 lines are used to style the slider
    label_style = cv_black → fontface → cv_fontsize(20)
    app = CV_TickLabelAppearance(; label_style)
    app_notick = CV_TickLabelAppearance(; tick_length=0, gap=15, label_style)
    # Odd values between 1 and 19 are labeled
    rulers=(app ↦ ("%.0f" ⇒ value_min:2:value_max),
            # Even values are not labeled           # Set the label 'n' to the slider left
            app ↦ ("" ⇒ value_min+1:2:value_max), app_notick ↦ ("n:" ⇒ 0.5,))

    return CV_ParamWithSlider(;
        slider_pos, value_min=value_min-0.1, value_max=value_max+0.1,
        rulers, set_slider_value_func)
end

# creates the slider object. The slider position
# can be gotten by `param_n.slider_pos.value`
param_n = get_param_n()

# Power series for the Arctan function
function trafo2(z)
    n = param_n.slider_pos.value
    result = 0.0
    for k in 1:(n÷2)
      result += (isodd(k) ? +1 : -1)*z^(2*k-1)
    end
    return result
end

trafo1 = z -> log((1+z*1im)/(1-z*1im))/(2im)

codomain1 = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)
codomain2 = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)

# Create custom axis ruler for the (codomain) portraits
label_style = cv_color(0,0,0) → fontface → cv_fontsize(20)
rulers=(CV_Ruler(cv_format_ticks("%.0f", -2.0:1.0:2.0...),
    CV_TickLabelAppearance(; label_style)),)

circle = (cv_black → cv_linewidth(3)) ↦ CV_LinePainter(
    cv_arc_lines(0.0, 2π, (1.0,)))

# Now we build the scene. We use `cv_scene_comp_codomains_std`
# to create a 'standard' left-right layout with slider(s) above
scene = cv_scene_comp_codomains_std((param_n, ), trafo1, trafo2,
    codomain1, codomain2; codomain1_re_rulers=rulers,
    painter1=CV_PortraitPainter(trafo1) → circle,
    painter2=CV_PortraitPainter(trafo2) → circle)
cv_get_redraw_func(scene)()

cv_save_image(cv_get_can_layout(scene), "./ArctanSeries.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4: