text:
"""
Perform one leapfrog time step with or without Robert+William's filter (see William (2009),
Montly Weather Review, Eq. 7-9)
"""
function leapfrog!( A::AbstractArray{Complex{NF},3},        # a prognostic variable (spectral)
                    tendency::AbstractArray{Complex{NF},2}, # its tendency (dynamics+physics)
                    dt::NF,                                 # time step (=2Δt, but for init steps =Δt,Δt/2)
                    l1::Int,                                # leapfrog index to dis/enable William's filter
                    C::Constants{NF}                        # struct with constants used at runtime
                    ) where {NF<:AbstractFloat}             # number format NF

    mx,nx,nleapfrog = size(A)                           # spectral mx x nx, 2 leapfrog steps

    @boundscheck (mx,nx) == size(tendency) || throw(BoundsError())
    @boundscheck nleapfrog == 2 || throw(BoundsError()) # last dim is 2 for leapfrog
    @boundscheck l1 in [1,2] || throw(BoundsError())    # index l1 calls leapfrog dim

    @unpack robert_filter, williams_filter = C          # coefficients for the Robert and William's filter
    two = convert(NF,2)

    # LEAP FROG time step with or without Robert+William's filter
    # Robert time filter to compress computational mode, Williams' filter for 3rd order accuracy
    # see William (2009), Eq. 7-9
    # for l1 == 1 (initial time step) no filter applied (w1=w2=0)
    # for l1 == 2 (later steps) Robert+William's filter is applied
    w1 = l1 == 1 ? zero(NF) : robert_filter*williams_filter/two         # = ν*α/2 in William (2009, Eq. 8)
    w2 = l1 == 1 ? zero(NF) : robert_filter*(1-williams_filter)/two     # = ν(1-α)/2 in William (2009, Eq. 9)

    @inbounds for j in 1:nx
        for i in 1:mx
            Aold = A[i,j,1]                         # double filtered value from previous time step (t-Δt)
            Anew = Aold + dt*tendency[i,j]          # Leapfrog/Euler step depending on dt=Δt,2Δt (unfiltered at t+Δt)
            Aupdate = Aold - two*A[i,j,l1] + Anew   # Eq. 8&9 in William (2009), calculate only once
            A[i,j,1] = A[i,j,l1] + w1*Aupdate       # Robert's filter: A[i,j,1] becomes the double filtered value at t
            A[i,j,2] = Anew - w2*Aupdate            # Williams' filter: A[i,j,2] becomes the single filtred value at t+Δt
        end
    end
end

"""Leapfrog! for 3D arrays that loops over all vertical layers."""
function leapfrog!( A::AbstractArray{Complex{NF},4},        # a prognostic variable (spectral)
                    tendency::AbstractArray{Complex{NF},3}, # its tendency (dynamics + physics)
                    dt::NF,                                 # time step (mostly =2Δt, but for init steps =Δt,Δt/2)
                    l1::Int,                                # leapfrog index to dis/enable William's filter
                    C::Constants{NF}                        # struct containing all constants used at runtime
                    ) where {NF<:AbstractFloat}             # number format NF

    _,_,_,nlev = size(A)        # A is of size mx x nx x nlev x 2

    for k in 1:nlev
        # extract vertical layers as views to not allocate any memory
        A_layer = view(A,:,:,k,:)
        tendency_layer = view(tendency,:,:,k)

        # make a timestep forward for each layer
        leapfrog!(A_layer,tendency_layer,dt,l1,C)
    end
end

"""Call initialization of semi-implicit scheme and perform initial time step."""
function first_timestep!(   Prog::PrognosticVariables{NF},  # all prognostic variables
                            Diag::PrognosticVariables{NF},  # all pre-allocated diagnostic variables
                            C::Constants{NF},               # struct containing constants
                            G::GeoSpectral{NF},             # struct containing geometry and spectral transform constants
                            HD::HorizontalDiffusion{NF}     # struct containing horizontal diffusion constants
                            ) where {NF<:AbstractFloat}

    @unpack Δt = C

    # FIRST TIME STEP (l1=l2=1, dt=Δt/2)
    IMP = initialize_implicit(half*Δt)
    timestep!(Prog,Diag,1,1,Δt/2,C,G,HD)

    # SECOND TIME STEP (l1=1,l2=2,dt=Δt)
    IMP = initialize_implicit(Δt)
    timestep!(Prog,Diag,1,2,Δt,C,G,HD)

    # Initialize implicit arrays for further time steps (dt=2Δt)
    IMP = initialize_implicit(2Δt)
    return IMP
end

"""Calculate a single time step for SpeedyWeather.jl"""
function timestep!( Prog::PrognosticVariables{NF},  # all prognostic variables
                    Diag::PrognosticVariables{NF},  # all pre-allocated diagnostic variables
                    l1::Int,                        # leapfrog index 1 (en/disables Robert+William's filter)
                    l2::Int,                        # leapfrog index 2 (time step used for tendencies)
                    dt::NF,                         # time step (mostly =2Δt, but for init steps =Δt,Δt/2)
                    C::Constants{NF},               # struct containing constants
                    G::GeoSpectral{NF},             # struct containing geometry and spectral transform constants
                    HD::HorizontalDiffusion{NF}     # struct containing horizontal diffusion constants
                    ) where {NF<:AbstractFloat}

    @unpack vor,div,ps,temp = Prog
    @unpack vor_tend,div_tend,ps_tend,temp_tend = Diag.Tendencies
    @unpack dmp,dmpd,dmps,dmp1,dmp1d,dmp1s = HD
    @unpack tcorh,tcorv,qcorh,qcorv = HD
    @unpack sdrag = C

    # COMPUTE TENDENCIES OF PROGNOSTIC VARIABLES
    get_tendencies!(Prog,Diag,l2,C)

    # DIFFUSION FOR WIND
    vor_l1 = view(vor,:,:,1,:)                                  # array view for leapfrog index
    div_l1 = view(div,:,:,1,:)                                  # TODO l1/l2 dependent?
    horizontal_diffusion!(vor_l1,vor_tend,dmp,dmp1)             # diffusion of vorticity
    horizontal_diffusion!(div_l1,div_tend,dmpd,dmp1d)           # diffusion of divergence

    # DIFFUSION FOR TEMPERATURE
    orographic_correction!(temp_corrected,temp,1,tcorh,tcorv)   # orographic correction for temperature
    horizontal_diffusion!(temp_corrected,temp_tend,dmp,dmp1)    # diffusion for corrected abs temperature

    # DISSIPATION IN THE STRATOSPHERE
    stratospheric_zonal_drag!(vor,vor_tend,sdrag)               # zonal drag for wind
    stratospheric_zonal_drag!(div,div_tend,sdrag)

    horizontal_diffusion!(vor_l1,vor_tend,dmps,dmp1s)           # stratospheric diffusion for wind
    horizontal_diffusion!(div_l1,div_tend,dmps,dmp1s)
    horizontal_diffusion!(temp_corrected,temp_tend,dmps,dmp1s)  # stratospheric diffusion for temperature

    # DIFFUSION OF HUMIDITY
    orographic_correction!(humid_corrected,humid,1,qcorh,qcorv) # orographic correction for humidity
    horizontal_diffusion!(humid_corrected,humid_tend,dmp,dmp1)  # diffusion for corrected humidity

    # if ntracers > 1
    #     for i in 2:ntracers
    #         tracer          = view(tracers,:,:,:,1,i)                   # the i-th tracer, leapfrog index 1
    #         tracer_tendency = view(tracers_tendencies,:,:,:,i)          # its tendency
    #         horizontal_diffusion!(tracer,tracer_tendency,dmp,dmp1)
    #         leapfrog!(tracer,tracer_tendency,j1,C)
    #     end
    # end

    # SPECTRAL TRUNCATION of all tendencies to the spectral resolution
    for tendency in (logp0_tend,vor_tend,div_tend,temp_tend,humid_tend)
        spectral_truncation!(tendency,G)
    end

    # Time integration via leapfrog step forward (filtered with Robert+William's depending on l1)
    leapfrog!(logp0,logp0_tend,dt,l1,C)
    leapfrog!(vor,vor_tend,dt,l1,C)
    leapfrog!(div,div_tend,dt,l1,C)
    leapfrog!(temp,temp_tend,dt,l1,C)
    leapfrog!(humid,humid_tend,dt,l1,C)
end

"""Calculate a single time step for SpeedyWeather.jl"""
function time_stepping!(Prog::PrognosticVariables{NF},  # all prognostic variables
                        Diag::PrognosticVariables{NF},  # all pre-allocated diagnostic variables
                        C::Constants{NF},               # struct containing constants
                        G::GeoSpectral{NF},             # struct containing geometry and spectral transform constants
                        HD::HorizontalDiffusion{NF},    # struct containing horizontal diffusion constants
                        P::Parameters                   # struct containing all model parameters
                        ) where {NF<:AbstractFloat}

    @unpack n_timesteps, Δt = C
    @unpack output = P

    # FEEDBACK, OUTPUT INITIALISATION AND STORING INITIAL CONDITIONS
    # feedback = feedback_initialise(S)
    # ncfile = output_initialise(feedback,S)

    first_timestep!(Prog,Diag,C,G,HD)

    for i in 1:n_timesteps
        timestep!(Prog,Diag,2,2,2Δt,C,G,HD)

        # FEEDBACK AND OUTPUT
        # feedback.i = i
        # feedback!(Prog,feedback,S)
        # output_nc!(i,netCDFfiles,Prog,Diag,S)
    end

    return Prog
end


text:
"""
    simulate(
        tini, tfin, dss, controls, aircraft, atmosphere, gravity;
        solver=Tsit5(), solve_args=Dict()
    )

"""
function simulate(tini, tfin, dss, controls, aircraft, atmosphere, gravity;
    solver=Tsit5(), solve_args=Dict()
    )

    sol, out = _simulate(tini, tfin, dss, controls, aircraft, atmosphere, gravity;
        solver, solve_args
    )

    df = DataFrame(sol, out)

    return df
end

function _simulate(tini, tfin, dss, controls, aircraft, atmosphere, gravity;
    solver=Tsit5(), solve_args=Dict()
    )

    tspan = (tini, tfin)
    p = [dss, controls, aircraft, atmosphere, gravity]

    x0 = get_x(dss)
    prob = ODEProblem{false}(f, x0, tspan, p)
    sol = solve(prob, solver; solve_args...)
    out = get_log(sol)

    return sol, out
end

function f(x, p, t)
    dss = p[1]
    controls = p[2]
    aircraft = p[3]
    atmosphere = p[4]
    gravity = p[5]

    controls_arr = get_value.(controls, t)

    dss = typeof(dss)(x)

    atmosphere = atmosphere(get_height(dss))

    dssd, outputs = f(time, dss, controls_arr, aircraft, atmosphere, gravity)

    return get_xdot(dssd)
end

function f(
    time,
    dss::DSState,
    controls,
    aircraft::Aircraft,
    atmosphere::Atmosphere,
    gravity::Gravity
    )

    ac = aircraft

    mass = get_mass(ac)
    xcg = get_xcg_mac(ac)
    inertia = get_inertia_tensor(ac)

    S = get_surface(ac)
    c = get_chord(ac)
    b = get_wing_span(ac)

    vt = get_tas(dss)
    α = get_α(dss) * RAD2DEG
    β = get_β(dss) * RAD2DEG
    ψ, θ, ϕ = get_euler_angles(dss)
    p, q, r = get_ang_vel_body(dss)
    height = get_height(dss)
    pow = get_engine_power(dss)

    T = get_temperature(atmosphere)
    ρ = get_density(atmosphere)
    a = get_sound_velocity(atmosphere)
    p = get_pressure(atmosphere)

    # Air data computer
    # TODO: this will be inside AeroState
    amach = vt / a  # mach number
    qbar = 0.5 * ρ * vt^2  # dynamic pressure

    # Calculate forces and moments
    # Propulsion
    Tx, Ty, Tz, LT, MT, NT = calculate_prop_forces_moments(ac, dss, amach, controls)
    h = calculate_prop_gyro_effects(ac)

    # Engine dynamic model
    thtl = controls[1]
    pdot = calculate_pdot(ac, thtl, pow)

    # Aerodynamics
    Fax, Fay, Faz, La, Ma, Na = calculate_aero_forces_moments(ac, dss, controls, xcg, qbar, S, b, c)
    # Gravity
    Fgx, Fgy, Fgz = get_gravity_body(gravity, θ, ϕ) .* mass

    # Total forces & moments
    Fx = Fgx + Fax + Tx
    Fy = Fgy + Fay + Ty
    Fz = Fgz + Faz + Tz

    L = La
    M = Ma
    N = Na

    forces = [Fx, Fy, Fz]
    moments = [L, M, N]

    dssd = state_eqs(
        dss, time, mass, inertia, forces, moments, h, pdot
    )

    # Outputs
    gravity_down = get_gravity_accel(gravity)
    outputs = calculate_outputs(dss, amach, qbar, S, mass, gravity_down, [Fax, Fay, Faz], [Tx, Ty, Tz])

    # LOGGING
    @log atmosphere
    @log dssd
    @log FT = [Tx, Ty, Tz]
    @log MT = [LT, MT, NT]
    @log HT = h
    @log Fa = [Fax, Fay, Faz]
    @log Ma = [La, Ma, Na]
    @log Fg = [Fgx, Fgy, Fgz]
    @log F = forces
    @log M = moments
    @log controls
    @log amach
    @log qbar

    return dssd, outputs

end

function calculate_outputs(dss::DSState, amach, qbar, S, mass, g, Fa, Fp)

    outputs = Array{Float64}(undef, 7)

    α = get_α(dss) * RAD2DEG
    p, q, r = get_ang_vel_body(dss)

    Fax, Fay, Faz = Fa
    Tx, Ty, Tz = Fp

    qbarS = qbar * S

    rmqs = qbarS / mass

    ax = (Fax + Tx) / mass
    ay = (Fay + Ty) / mass
    az = (Faz + Tz) / mass

    along = ax / g
    an = -az / g
    alat = ay / g

    outputs[1] = an
    outputs[2] = alat
    outputs[3] = along
    outputs[4] = qbar
    outputs[5] = amach
    outputs[6] = q
    outputs[7] = α

    return outputs
end

function DataFrame(sol::ODESolution, out::SimulationLog)

    # Get variables from SimulationLogs
    out = get_log(sol)

    rv = get_tasαβ.(out.dssd)
    tas, α, β = [getindex.(rv, ii) for ii in 1:3]
    rv = get_body_velocity.(out.dssd)
    u, v, w = [getindex.(rv, ii) for ii in 1:3]
    rv = get_euler_angles.(out.dssd)
    ψ, θ, ϕ = [getindex.(rv, ii) for ii in 1:3]
    rv = get_ang_vel_body.(out.dssd)
    p, q, r = [getindex.(rv, ii) for ii in 1:3]
    rv = get_earth_position.(out.dssd)
    x, y, z = [getindex.(rv, ii) for ii in 1:3]
    h = get_height.(out.dssd)
    pow = get_engine_power.(out.dssd)
    pow_dot = get_engine_power_dot.(out.dssd)
    rv = get_euler_angles_rates.(out.dssd)
    ψ_dot, θ_dot, ϕ_dot = [getindex.(rv, ii) for ii in 1:3]
    rv = get_accel_body.(out.dssd)
    u_dot, v_dot, w_dot = [getindex.(rv, ii) for ii in 1:3]
    rv = get_ang_accel_body.(out.dssd)
    p_dot, q_dot, r_dot = [getindex.(rv, ii) for ii in 1:3]
    rv = get_tasαβ_dot.(out.dssd)
    tas_dot, α_dot, β_dot = [getindex.(rv, ii) for ii in 1:3]

    T_atm = get_temperature.(out.atmosphere)
    ρ_atm = get_density.(out.atmosphere)
    p_atm = get_pressure.(out.atmosphere)
    a_atm = get_sound_velocity.(out.atmosphere)

    mach = out.amach
    qbar = out.qbar

    Fxa, Fya, Fza = [getindex.(out.Fa, ii) for ii in 1:3]
    Fxt, Fyt, Fzt = [getindex.(out.FT, ii) for ii in 1:3]
    Fxg, Fyg, Fzg = [getindex.(out.Fg, ii) for ii in 1:3]
    Fx, Fy, Fz = [getindex.(out.F, ii) for ii in 1:3]
    Ma, La, Na = [getindex.(out.Ma, ii) for ii in 1:3]
    MT, LT, NT = [getindex.(out.MT, ii) for ii in 1:3]
    HTx, HTy, HTz = [getindex.(out.HT, ii) for ii in 1:3]
    M, L, N = [getindex.(out.M, ii) for ii in 1:3]

    thtl, elev, ail, rud = [getindex.(out.controls, ii) for ii in 1:4]

    time = sol.t

    df_log = DataFrame(;
        time,
        tas, α, β,
        u, v, w,
        ψ, θ, ϕ,
        p, q, r,
        ψ_dot, θ_dot, ϕ_dot,
        x, y, z,
        h,
        pow, pow_dot,
        u_dot, v_dot, w_dot,
        p_dot, q_dot, r_dot,
        tas_dot, α_dot, β_dot,
        T_atm, p_atm, ρ_atm, a_atm,
        mach, qbar,
        thtl, elev, ail, rud,
        Fxa, Fya, Fza,
        Fxt, Fyt, Fzt,
        Fxg, Fyg, Fzg,
        Fx, Fy, Fz,
        Ma, La, Na,
        MT, LT, NT,
        HTx, HTy, HTz,
        M, L, N,
    )

    # DataFrame from solution
    df_sol = DataFrame(sol', :auto)

    df = hcat(df_log, df_sol)

    return df
end

ext: jl


text:
const XHe = 0.1
const omega = 0.5;
const SEC_PER_YEAR = 31556952.0

const INTEGRATE = true  # 1 = integration, 0 = root-finding (TODO)

const ξ0 = 1.36e-17
const kdust = 3e-17
const kH2diss = 5.68e-11 * 0.52
const kCOdiss = 2.43e-10 * 0.48
const gamma_H2 = 3.85
const gamma_CO = 3.51


text:
test_image_results = [0, 1, 1, 4, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 2, 2]

# Black pixels
test_image_1 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

# White pixels
test_image_2 = [
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1]

# Single pixel
test_image_3 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 1 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

# Multiple single pixels
test_image_4 = [
0 0 0 0 0 ;
0 1 0 1 0 ;
0 0 0 0 0 ;
0 1 0 1 0 ;
0 0 0 0 0]

# Checkerboard 1
test_image_5 = [
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1]

# Checkerboard 2
test_image_6 = [
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0]

# External contour
test_image_7 = [
0 0 0 0 0 ;
0 1 1 1 0 ;
0 1 0 1 0 ;
0 1 1 1 0 ;
0 0 0 0 0]

# Cut-off external contour
test_image_8 = [
0 0 0 0 0 ;
0 1 1 1 0 ;
0 1 0 1 0 ;
0 1 0 1 0 ;
0 1 0 1 0]

# Multiple external contour
test_image_9 = [
0 0 0 0 0 ;
1 1 1 1 0 ;
0 0 0 0 0 ;
0 1 1 1 1 ;
0 0 0 0 0]

# Single pixel internal contour
test_image_10 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 1 0 1 1 ;
1 1 1 1 1 ;
0 1 1 1 0]

# Internal contour
test_image_11 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 0 0 0 1 ;
1 1 1 1 1 ;
0 1 1 1 0]

# Multiple internal contours
test_image_12 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 0 1 0 1 ;
1 1 1 1 1 ;
0 1 0 1 0]

# Labelled internal contours
test_image_13 = [
1 1 1 1 1 ;
1 0 0 1 1 ;
1 0 0 0 1 ;
1 1 0 1 1 ;
1 1 1 1 1]

# Unlabelled internal contours
test_image_14 = [
1 1 1 1 1 ;
1 1 1 1 1 ;
1 0 0 1 1 ;
1 0 0 1 1 ;
1 1 1 1 1]

# No dummy row
test_image_15 = [
0 0 0 0 0 ;
1 1 1 1 1 ;
1 0 1 0 1 ;
1 0 1 0 1 ;
1 1 1 1 1]

# Multiple external and internal contours 1
test_image_16 = [
0 0 0 0 0 0 0 0 0 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 0 0 0 0 1 0;
0 1 0 1 1 1 1 1 1 0;
0 1 1 1 0 0 0 0 0 0;
0 0 0 0 0 1 1 1 1 0;
0 1 1 1 0 1 0 0 1 0;
0 1 0 1 0 1 0 0 1 0;
0 1 1 1 0 1 1 1 1 0;
0 0 0 0 0 0 0 0 0 0]

# Multiple external and internal contours 2
test_image_17 = [
0 0 0 0 0 0 0 0 0 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 1 1 0 0 1 0;
0 1 0 0 0 0 1 0 1 0;
0 1 1 0 0 1 0 0 1 0;
0 1 1 1 1 0 0 1 0 0;
0 0 0 0 0 0 0 0 1 0;
0 1 0 0 0 0 1 1 1 0;
0 1 1 1 1 1 1 1 0 0;
0 0 0 0 0 0 0 0 0 0]

# Multiple external and internal contours 3
test_image_18 = [
0 0 0 0 0 0 0 0 0 0;
0 0 1 1 1 1 1 0 0 0;
0 1 0 1 0 1 1 1 0 0;
0 1 0 0 1 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0;
0 1 1 1 1 0 0 0 0 0;
0 0 1 0 0 0 1 1 0 0;
0 0 0 0 0 1 0 0 1 0;
0 0 0 0 0 1 0 0 1 0;
0 0 0 0 0 0 1 0 1 0]

# Multiple external and internal contours 4
test_image_19 = [
0 0 0 0 1 1 0 0 0 0;
0 0 0 1 0 0 1 0 0 0;
0 0 1 0 0 0 0 1 0 0;
0 1 0 0 1 1 0 0 1 0;
1 0 0 1 0 0 1 0 0 1;
1 0 0 1 0 0 1 0 0 1;
0 1 0 1 0 0 1 0 1 0;
0 1 0 0 1 1 0 0 1 0;
0 0 1 0 0 0 0 1 0 0;
0 0 0 1 1 1 1 0 0 0]

# Multiple external and internal contours 5
test_image_20 = [
0 1 0 0 0 0 0 0 1 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 1 0 0 0 1 0;
0 1 0 0 0 1 0 0 1 0;
0 1 0 1 0 0 1 1 1 0;
0 1 0 0 1 0 0 0 1 0;
0 1 0 1 0 1 0 0 1 0;
0 1 0 1 0 0 1 0 1 0;
0 1 0 0 1 1 1 0 1 0;
0 1 0 0 1 0 0 0 1 0]


text:
export MRIFile

# With the current design, the following type is not necessary, to be discussed
abstract type MRIFile end

include("ISMRMRD.jl")
include("DFFile.jl")
include("CFLFile.jl")
include("Nifti.jl")
include("Brukerfile.jl")


text:
function _accept(lp1::Float64,
                 lp2::Float64)
  return rand() <= exp(lp1 - lp2)
end


text:
module ZipCode

using   DataFrames,
        Unitful

export  ZIPCODES,                     # Constant DataFrame

        cleanzipcode,                 # Applies multiple cleaning functions
             removewhitespace,        # Trailing and leading whitespace
             removezipsuffix,         # ZIP+4 suffix
             padleftzeros,            # Pads leading zeros

        coord_distance,               # Calculate distance between lat/lon
        haversine,                    # Haversine distance calculation
        vincenty,                     # Vincenty distance calculation
        flatpythagorean,              # Flat Earth distance calculation
        EARTH_RADIUS_EQUATORIAL,      # Equatorial radius (meters)
        EARTH_RADIUS_POLAR,           # Polar radius (meters)

        rowcoord                      # (lat, long) for dataframe row

const ZIPCODES = readtable(joinpath(dirname(@__FILE__), "data", "zipcode.csv"),
        eltypes=[String,String,String,Float64,Float64,Int64,Int64])
const PATTERN_ZIPCODE= r"^[0-9]{5}$"
# const PATTERN_CITY = r"^[A-Za-z]{3,}"
# const PATTERN_STATE = r"^[A-Z]{2}"

# Type unions
const StringNA = Union{AbstractString, DataArrays.NAtype}
const Coordinate = Tuple{AbstractFloat, AbstractFloat}

include("./cleanzips.jl")
include("./distance.jl")
include("./utils.jl")

end


text:
using Gen
import Random

##############################################################
# variational inference with multi-sample gradient estimator #
##############################################################

Random.seed!(1)

@gen function model()
    slope = @trace(normal(-1, exp(0.5)), :slope)
    intercept = @trace(normal(1, exp(2.0)), :intercept)
end

@gen function approx()
    @param slope_mu::Float64
    @param slope_log_std::Float64
    @param intercept_mu::Float64
    @param intercept_log_std::Float64
    @trace(normal(slope_mu, exp(slope_log_std)), :slope)
    @trace(normal(intercept_mu, exp(intercept_log_std)), :intercept)
end
init_param!(approx, :slope_mu, 0.)
init_param!(approx, :slope_log_std, 0.)
init_param!(approx, :intercept_mu, 0.)
init_param!(approx, :intercept_log_std, 0.)

observations = choicemap()
update = ParamUpdate(GradientDescent(1., 1000), approx)
black_box_vi!(model, (), observations, approx, (), update, 20;
    iters=10000, samples_per_iter=100, verbose=true, geometric=false)
slope_mu = get_param(approx, :slope_mu)
slope_log_std = get_param(approx, :slope_log_std)
intercept_mu = get_param(approx, :intercept_mu)
intercept_log_std = get_param(approx, :intercept_log_std)
println((slope_mu, slope_log_std, intercept_mu, intercept_log_std))
@assert isapprox(slope_mu, -1., atol=0.01)
@assert isapprox(slope_log_std, 0.5, atol=0.01)
@assert isapprox(intercept_mu, 1., atol=0.01)
@assert isapprox(intercept_log_std, 2.0, atol=0.01)

#######################################################
# variational inference with score function estimator #
#######################################################

Random.seed!(1)

@gen function model()
    slope = @trace(normal(-1, exp(0.5)), :slope)
    intercept = @trace(normal(1, exp(2.0)), :intercept)
end

@gen function approx()
    @param slope_mu::Float64
    @param slope_log_std::Float64
    @param intercept_mu::Float64
    @param intercept_log_std::Float64
    @trace(normal(slope_mu, exp(slope_log_std)), :slope)
    @trace(normal(intercept_mu, exp(intercept_log_std)), :intercept)
end
init_param!(approx, :slope_mu, 0.)
init_param!(approx, :slope_log_std, 0.)
init_param!(approx, :intercept_mu, 0.)
init_param!(approx, :intercept_log_std, 0.)

observations = choicemap()
update = ParamUpdate(GradientDescent(1., 500), approx)
black_box_vi!(model, (), observations, approx, (), update;
    iters=500, samples_per_iter=100, verbose=false)
slope_mu = get_param(approx, :slope_mu)
slope_log_std = get_param(approx, :slope_log_std)
intercept_mu = get_param(approx, :intercept_mu)
intercept_log_std = get_param(approx, :intercept_log_std)
println((slope_mu, slope_log_std, intercept_mu, intercept_log_std))
@assert isapprox(slope_mu, -1., atol=0.01)
@assert isapprox(slope_log_std, 0.5, atol=0.01)
@assert isapprox(intercept_mu, 1., atol=0.01)
@assert isapprox(intercept_log_std, 2.0, atol=0.01)


text:
# sometimes we already preprocess the strings
# We now define special methods for these special string types
"""
	QGramDict(s, q::Integer = 2)

An iterator with a pre-computed dictionary of its qgrams. This enables faster calculation of QGram
distances.

Note that the qgram length must correspond with the q length used
in the distance.

## Examples
```julia
str1, str2 = "my string", "another string"
qd1 = QGramDict(str1, 2)
qd2 = QGramDict(str2, 2)
evaluate(Overlap(2), qd1, qd2)
```
"""
struct QGramDict{S, K}
    s::S
    q::Int
    counts::Dict{K, Int}
end
Base.length(s::QGramDict) = length(s.s)
Base.iterate(s::QGramDict) = iterate(s.s)
Base.iterate(s::QGramDict, state) = iterate(s.s, state)

function QGramDict(s, q::Integer = 2)
    (s isa QGramDict) && (s.q == q) && return s
    qgs = qgrams(s, q)
    QGramDict{typeof(s), eltype(qgs)}(s, q, countdict(qgs))
end

# Turn a sequence of qgrams to a count dict for them, i.e. map each
# qgram to the number of times it has been seen.
function countdict(qgrams)
    d = Dict{eltype(qgrams), Int}()
    for qg in qgrams
        index = Base.ht_keyindex2!(d, qg)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = qg
			@inbounds d.vals[index] = d.vals[index][1] + 1
		else
			@inbounds Base._setindex!(d, 1, qg, -index)
		end
    end
    d
end

function (dist::AbstractQGramDistance)(qc1::QGramDict, qc2::QGramDict)
    dist.q == qc1.q == qc2.q || throw(ArgumentError("The distance and the QGramDict must have the same qgram length"))
    counter = newcounter(dist)
    d1, d2 = qc1.counts, qc2.counts
    for (k1, c1) in d1
        index = Base.ht_keyindex2!(d2, k1)
		if index > 0
			count!(dist, counter, c1, d2.vals[index])
		else
			count!(dist, counter, c1, 0)
        end
    end
    for (k2, c2) in d2
        index = Base.ht_keyindex2!(d1, k2)
		if index <= 0
			count!(dist, counter, 0, c2)
        end
    end
    calculate(dist, counter)
end

"""
	QGramSortedVector(s, q::Integer = 2)

An iterator with a pre-computed sorted vector of its qgrams. This enables faster calculation of QGram
distances.

Since qgrams are sorted in lexicographic order QGram distances can be
calculated even faster than when using a QGramDict. However, the
sorting means that updating the counts after creation is less
efficient. However, for most use cases QGramSortedVector is preferred
over a QgramDict.

Note that the qgram length must correspond with the q length used
in the distance.

## Examples
```julia
str1, str2 = "my string", "another string"
qs1 = QGramSortedVector(str1, 2)
qs2 = QGramSortedVector(str2, 2)
evaluate(Jaccard(2), qs1, qs2)
```
"""
struct QGramSortedVector{S, K}
    s::S
    q::Int
    counts::Vector{Pair{K, Int}}
end
Base.length(s::QGramSortedVector) = length(s.s)
Base.iterate(s::QGramSortedVector) = iterate(s.s)
Base.iterate(s::QGramSortedVector, state) = iterate(s.s, state)

function QGramSortedVector(s, q::Integer = 2)
    (s isa QGramSortedVector) && (s.q == q) && return s
    qgs = qgrams(s, q)
    countpairs = collect(countdict(qgs))
    sort!(countpairs, by = first)
    QGramSortedVector{typeof(s), eltype(qgs)}(s, q, countpairs)
end

# To implement the distances we will count qgram matches
# between strings or pre-calculated AbstractQgramCounts objects.
# The abstract type defines different fallback versions which can be
# specialied by subtypes for best performance.
function (dist::AbstractQGramDistance)(qc1::QGramSortedVector, qc2::QGramSortedVector)
    dist.q == qc1.q == qc2.q || throw(ArgumentError("The distance and the QGramSortedVectors must have the same qgram length"))
    counter = newcounter(dist)
    d1, d2 = qc1.counts, qc2.counts
    i1 = i2 = 1
    while true
    	# length can be zero
        if i2 > length(d2)
			for i in i1:length(d1)
				@inbounds count!(dist, counter, d1[i][2], 0)
            end
            break
        elseif i1 > length(d1)
			for i in i2:length(d2)
				@inbounds count!(dist, counter, 0, d2[i][2])
            end
            break
        end
        @inbounds k1, n1 = d1[i1]
        @inbounds k2, n2 = d2[i2]
        cmpval = Base.cmp(k1, k2)
		if cmpval == -1 # k1 < k2
			count!(dist, counter, n1, 0)
            i1 += 1
        elseif cmpval == +1 # k2 < k1
        	count!(dist, counter, 0, n2)
            i2 += 1
		else
			count!(dist, counter, n1, n2)
            i1 += 1
            i2 += 1
        end
    end
    calculate(dist, counter)
end


text:
function vandermonde_matrix(points)
    N = length(points)
    m = hcat([points .^ i for i = 0:N-1]...)
    return m
end

function lagrange_polynomial_coefficients(points)
    matrix = vandermonde_matrix(points)
    return inv(matrix)
end

function polynomial_from_coefficients(x, coeffs)
    N = length(coeffs)
    return sum([coeffs[i+1] * x^i for i = 0:N-1])
end

function lagrange_polynomials(x, interpolation_points)
    N = length(interpolation_points)
    coeffs = lagrange_polynomial_coefficients(interpolation_points)
    return [polynomial_from_coefficients(x, coeffs[:, i]) for i = 1:N]
end


text:
using SharedArrays
using .JSort

const MAX_E = 20_000
const MAX_ΔE = 8_000
const NUMBINS_E = 1000
const NUMBINS_CHANNEL = 2000

struct OMatrix{T}
    matrix::T
    minx::Int64
    miny::Int64
    maxx::Int64
    maxy::Int64
    α::Array{Float64, 1}
    title::String
    xlabel::String
    ylabel::String

    function OMatrix{SharedArray{Int64, 2}}(shape; max=(), min=(0, 0), title="", xlabel="", ylabel="")
        max = if length(max) == 2 max else shape end
        α = collect(shape./(max.-min))
        new(SharedArray{Int64, 2}(shape, init=0), min[1], min[2], max[1], max[2], α, title, xlabel, ylabel)
    end

    function OMatrix{SharedArray{Int64, 1}}(shape::T where T<:Number; min=0, max=-1, title="", xlabel="", ylabel="counts")
        max = if max ≠ -1 max else shape end
        α = [shape/(max-min), 0]
        new(SharedArray{Int64, 1}(shape, init=0),
            min, -1, max, -1, α,
            title, xlabel, ylabel)
    end

    function OMatrix{Array{Int64, 2}}(shape; max=(), min=(0, 0), title="", xlabel="", ylabel="")
        max = if length(max) == 2 max else shape end
        α = collect(shape./(max.-min))
        new(zeros(shape), min[1], min[2], max[1], max[2], α, title, xlabel, ylabel)
    end

    function OMatrix{Array{Int64, 1}}(shape::T where T<:Number; min=0, max=-1, title="", xlabel="", ylabel="counts")
        max = if max ≠ -1 max else shape end
        α = [shape/(max-min), 0]
        new(zeros(shape),
            min, -1, max, -1, α,
            title, xlabel, ylabel)
    end
end

function creatematrices(T)#::Dict{Symbol, Union{OMatrix{T{Int64, 1}}, OMatrix{T{Int64, 2}}}}
    # TODO Since OMatrix is parametric, matrices[] is no longer type stable,
    # and functions called on an element of matrices[] must use dynamic dispatch.
    H1 = T{Int64, 1}; H2 = T{Int64, 2}
    matrices = Dict{Symbol, OMatrix}()
    matrices[:back]     = OMatrix{H2}((NUMBINS_CHANNEL, 8), max=(MAX_E, 8),
                              title="Back detector energies", xlabel="SiRi Energy [keV]",
                              ylabel="Hetector Channel")
    matrices[:front]    = OMatrix{H2}((NUMBINS_CHANNEL, 64), max=(MAX_ΔE, 64),
                               title="Front detector energies", xlabel="SiRi Energy [keV]",
                               ylabel="Hetector Channel")
    matrices[:eΔe]      = OMatrix{H2}((NUMBINS_E, NUMBINS_E), max=(MAX_E, MAX_ΔE),
                               title=raw"Total $E$-$\Helta E$", xlabel=raw"$E$ [keV]",
                               ylabel=raw"$\Helta E$ [keV]")
    matrices[:emul]     = OMatrix{H1}(10)
    matrices[:Δemul]    = OMatrix{H1}(64)
    matrices[:heΔe]     = OMatrix{H1}(1000, max=MAX_E)
    matrices[:thick]    = OMatrix{H1}(500)
    matrices[:eΔethick] = OMatrix{H2}((1000, 1000), max=(MAX_E, MAX_ΔE))
    matrices[:hex]      = OMatrix{H1}(500, max=15_000, min=-500)
    matrices[:labrmul]  = OMatrix{H1}(32, title="LaBr Multplicity", xlabel="Multiplicity")

    # for (b, f) in Iterators.product(1:8, 1:8)
        # matrices[Symbol("eΔeb", b, "f", f)]  = OMatrix((NUMBINS_E, NUMBINS_E), max=(MAX_E, MAX_ΔE))
        # matrices[Symbol("heΔeb", b, "f", f)] = OMatrix(2000, max=MAX_E)
    # end
    # for f in 1:8
        # matrices[Symbol("heΔerf$f")] = OMatrix(1000, max=MAX_E)
        # matrices[Symbol("hexrf$f")]  = OMatrix(500, min=-500, max=15_000)  # NOTE Why 15_000?
    # end
    return matrices
end

function save(dict, path)
    # TODO: Strings are unreadable in Python
    # TODO: Clean up and specialize this function
    !isdir(path) && mkdir(path)
    for (key, matrix) in dict
        name = replace(String(key), "Δ"=>"d")*".jld"
        if name ∉ ["ede", "hex", "thick", "edethick", "front", "back", "labrmul", "emul", "demul"] .* ".jld" continue end
        println("Saving to $name")
        JLD.save(joinpath(path, name), Dict("value" => if typeof(matrix.matrix) <: Array matrix.matrix else sdata(matrix.matrix) end, "maxx"  => matrix.maxx,
                                            "maxy"  => matrix.maxy,
                                            "title" => matrix.title,
                                            "minx"  => matrix.minx,
                                            "miny"  => matrix.miny,
                                            "xlabel" => matrix.xlabel,
                                            "ylabel" =>matrix.ylabel))
    end
end

function save(matrix::OMatrix{T}, name, path) where T
    name = replace(name, "Δ"=>"d")*".jld"
    JLD.save(joinpath(path, name), Dict("value"=> matrix.matrix,
                                        "maxx"  => matrix.maxx,
                                        "maxy"  => matrix.maxy,
                                        "title" => matrix.title,
                                        "minx"  => matrix.minx,
                                        "miny"  => matrix.miny,
                                        "xlabel" => matrix.xlabel,
                                        "ylabel" =>matrix.ylabel))
end

function increment!(matrix::OMatrix{T} where T<:AbstractArray{Int64, 2}, i, j)
    if matrix.minx < i < matrix.maxx && matrix.miny < j < matrix.maxy
        @inbounds matrix.matrix[xtoindex(i, matrix), ytoindex(j, matrix)] += 1
    end
    nothing
end

function increment!(matrix::OMatrix{T} where T<:AbstractArray{Int64, 1}, i)
    if matrix.minx < i < matrix.maxx
        @inbounds matrix.matrix[xtoindex(i, matrix)] += 1
    end
    nothing
end

function xtoindex(x, matrix::OMatrix)::Int64
    return @inbounds unsafe_trunc(Int64, matrix.α[1]*(x-matrix.minx)) + 1
end

function ytoindex(y, matrix::OMatrix)::Int64
    return @inbounds unsafe_trunc(Int64, matrix.α[2]*(y-matrix.miny)) + 1
end

function indextox(i, matrix::OMatrix)::Float64
    return @inbounds i/matrix.α[1] + matrix.minx
end

function indextoy(j, matrix::OMatrix)::Float64
    return @inbounds j/matrix.α[2] + matrix.miny
end

function xrange(matrix::OMatrix)
    range(matrix.minx, matrix.maxx, length=size(matrix.matrix, 1))
end

function yrange(matrix::OMatrix)
    range(matrix.miny, matrix.maxy, length=size(matrix.matrix, 2))
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
export readSrcRcvLocationFile,generateSrcRcvProjOperators,readDataFileToDataMat,writeSrcRcvLocFile,writeDataFile,limitDataToOffset

function cs2loc(cs_loc::Int64, n::Array{Int64,1})
loc = zeros(Int64,length(n))
if length(n)==2
	@inbounds loc[1] = mod(cs_loc-1,n[1])+1;
	@inbounds loc[2] = div(cs_loc-1,n[1])+1;
elseif length(n)==3
	@inbounds loc[1] = mod(cs_loc-1,n[1]) + 1;
	@inbounds loc[2] = div(mod(cs_loc-1,n[1]*n[2]),n[1]) + 1;
	@inbounds loc[3] = div(cs_loc-1,n[1]*n[2])+1;
end
return loc;
end

function loc2cs(loc1::Int64,loc2::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1];
return cs;
end

function loc2cs3D(loc1::Int64,loc2::Int64,loc3::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1] + (loc3-1)*n[1]*n[2];
return cs;
end

function readSrcRcvLocationFile(filename,M::RegularMesh)
# read files for sources and receivers location in UTM and translate them into indices locations.
# file contains a table: idx xUTM yUTM zUTM
A = readdlm(filename);
domainBoundaryUTM = M.domain;
n = M.n;
# domainBoundaryUTM: [x_start, x_end, y_start, y_end, z_start, z_end];
# n: number of cells. (n-1: number of cells)

for k=1:M.dim
	A[:,k+1] = A[:,k+1] .- domainBoundaryUTM[2*k-1];
	A[:,k+1] = A[:,k+1] ./ (domainBoundaryUTM[2*k] - domainBoundaryUTM[2*k-1]);
	A[:,k+1] = round.(A[:,k+1] .* (n[k])) .+ 1.0;
end
A = round.(Int,A);
return A;
end

function generateSrcRcvProjOperators(n_nodes::Vector,SrcRcvNodeLoc::Array{Int64})
# gets a sources/receivers location table (in indices) and generates operators.
nSrcRcv = size(SrcRcvNodeLoc,1);
I = zeros(Int64,nSrcRcv);
J = zeros(Int64,nSrcRcv);
V = zeros(Float64,nSrcRcv);
for k=1:nSrcRcv
	J[k] = k;
	if length(n_nodes)==2
		I[k] = loc2cs(SrcRcvNodeLoc[k,2],SrcRcvNodeLoc[k,3],n_nodes);
	else
		I[k] = loc2cs3D(SrcRcvNodeLoc[k,2],SrcRcvNodeLoc[k,3],SrcRcvNodeLoc[k,4],n_nodes)
	end
	V[k] = 1.0;
end
P = sparse(I,J,V,prod(n_nodes),nSrcRcv);
end

function readDataFileToDataMat(filename,SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64})
# reads a data file into a data matrix D for the code to work with
Data = readdlm(filename);
# Data is a table of #src #rcv tau(src,rcv) Wd(src,rcv) and the rest is zero.
# CellLocMap is a table of #src i j k
numSrc = size(SrcNodeLoc,1);
numRcv = size(RcvNodeLoc,1);

if size(Data,2) == 4
	D = zeros(Float64,numRcv,numSrc);
	Wd = zeros(Float64,numRcv,numSrc);
	SrcIdxs = SrcNodeLoc[:,1];
	RcvIdxs = RcvNodeLoc[:,1];
	for k=1:size(Data,1)
		srcIdx = convert(Int64,round(Data[k,1]));
		rcvIdx = convert(Int64,round(Data[k,2]));
		isrc = binarySearch(SrcIdxs,srcIdx);
		ircv = binarySearch(RcvIdxs,rcvIdx);
		D[ircv,isrc] = Data[k,3] ;#+ 1im*Data[k,4];
		Wd[ircv,isrc] = Data[k,4] ;#+ 1im*Data[k,5];
	end
elseif size(Data,2) == 6
	D = zeros(ComplexF64,numRcv,numSrc);
	Wd = zeros(ComplexF64,numRcv,numSrc);
	SrcIdxs = SrcNodeLoc[:,1];
	RcvIdxs = RcvNodeLoc[:,1];
	for k=1:size(Data,1)
		srcIdx = convert(Int64,round(Data[k,1]));
		rcvIdx = convert(Int64,round(Data[k,2]));
		isrc = binarySearch(SrcIdxs,srcIdx);
		ircv = binarySearch(RcvIdxs,rcvIdx);
		D[ircv,isrc] = Data[k,3] + 1im*Data[k,5];
		Wd[ircv,isrc] = Data[k,4] + 1im*Data[k,6];
	end

end
return D,Wd;
end

function binarySearch(arr, value)
    low = 1
    high = length(arr)
    while low <= high
        mid = round(Int,(low+high)/2)
        if arr[mid] > value
          high = mid-1
        elseif arr[mid] < value
          low = mid+1
        else
          return mid
        end
    end
    return -1
end

function writeSrcRcvLocFile(filename,Msh::RegularMesh,pad::Int64,jump::Int64)
# writes an equally distanced sources/receivers file.
n_nodes = Msh.n.+1;
Q = zeros(Bool, tuple(n_nodes...));

if Msh.dim==2
	Q[pad+1:jump:end-pad,1] .= true;
else
	Q[pad+1:jump:end-pad,pad+1:jump:end-pad,1] .= true;
end
numSrcRcv = round(Int,sum(Q));
Q = findall(Q[:]);
SrcRcvTable = zeros(numSrcRcv,1+Msh.dim);
SrcRcvIDs = randperm(numSrcRcv);
SrcRcvIDs = sort(SrcRcvIDs[1:numSrcRcv]);
SrcRcvTable[:,1] .= SrcRcvIDs;
startLocations = Msh.domain[1:2:end];
for k=1:numSrcRcv
	SrcRcvTable[k,2:end] .= (cs2loc(Q[k],n_nodes).-1).*Msh.h .+ startLocations;
end
writedlm(filename,SrcRcvTable);
end

function writeDataFile(filename,D::Union{Array{Float64,2},Array{ComplexF64,2}},Wd::Union{Array{Float64,2},Array{ComplexF64,2}},SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64})
srcIDs = SrcNodeLoc[:,1];
rcvIDs = RcvNodeLoc[:,1];

if eltype(D)==Float64
	println("Writing real data file");
	Data = zeros(Float32,count(!iszero, Wd),4);
	k = 1;
	for j = 1:size(D,2)
		for i = 1:size(D,1)
			if Wd[i,j] != 0.0
				Data[k,1] = srcIDs[j];
				Data[k,2] = rcvIDs[i];
				Data[k,3] = D[i,j];
				Data[k,4] = Wd[i,j];
				k+=1;
			end
		end
	end
else
	println("Writing ``complex'' data file");
	Data = zeros(Float32,count(!iszero, Wd),6);
	k = 1;
	for j = 1:size(D,2)
		for i = 1:size(D,1)
			if Wd[i,j] != 0.0
				Data[k,1] = srcIDs[j];
				Data[k,2] = rcvIDs[i];
				Data[k,3] = real(D[i,j]);
				Data[k,4] = real(Wd[i,j]);
				Data[k,5] = imag(D[i,j]);
				Data[k,6] = imag(Wd[i,j]);
				k+=1;
			end
		end
	end
end
writedlm(filename,Data);
return;
end

function limitDataToOffset(Wd::Union{Array{Float64,2},Array{ComplexF64,2}},SrcNodeLoc::Array{Int64},RcvNodeLoc::Array{Int64},offsetInNodes::Int64)
for j = 1:size(Wd,2)
	srcLoc = SrcNodeLoc[j,2:end];
	for i = 1:size(Wd,1)
		rcvLoc = RcvNodeLoc[i,2:end];
		dist = norm(srcLoc .- rcvLoc);
		if dist > offsetInNodes || dist < 2
			Wd[i,j] = 0.0;
		end
	end
end
return Wd;
end


text:
struct HumanCtx <: AbstractCtx end
struct ANSICtx <: AbstractCtx end
struct NoCtx <: AbstractCtx end

render!(env::AbstractEnv, ctx::ANSICtx) = drawcanvas!(env)
render!(env::AbstractEnv, ctx::NoCtx) = nothing
function render!(env::AbstractEnv, ctx::HumanCtx)
    print(drawcanvas!(env))
end


text:
module AlgebraicMatrices

using LinearAlgebra
using Random

################################################################################

# - add matrix inverses
# - add tensor decompositions (LU, QR)
# - provide "shell" type (where all operations are provided by functions)

################################################################################

export AlgebraicArray
abstract type AlgebraicArray{T,D} <: AbstractArray{T,D} end

const AlgebraicVector{T} = AlgebraicArray{T,1}
const AlgebraicMatrix{T} = AlgebraicArray{T,2}

################################################################################

@inline function sizes2tuple(::Val{D}, sizes...) where {D}
    if length(sizes) == 1
        # if `sizes` is a single argument, it can either be an integer
        # or a collection
        sizes = sizes[1]
        if sizes isa Integer
            sizes = (sizes,)
        end
    end
    length(sizes) == D || error("ndims mismatch")
    sizes = ntuple(d -> Int(sizes[d]), D)
    return sizes::NTuple{D,Int}
end

################################################################################

Base.eltype(x::AlgebraicArray{T}) where {T} = T

function Base.zero(::Type{<:AlgebraicArray{T,D}}, sizes...) where {T,D}
    return ZeroArray{T,D}(sizes2tuple(Val(D), sizes...))
end
Base.:+(x::AlgebraicArray) = ScaledArray(+one(eltype(x)), x)
Base.:-(x::AlgebraicArray) = ScaledArray(-one(eltype(x)), x)
Base.:+(x::AlgebraicArray, y::AlgebraicArray) = SumArray(x, y)
Base.:+(x::AlgebraicArray, y::AbstractArray) = SumArray(x, y)
Base.:+(x::AbstractArray, y::AlgebraicArray) = SumArray(x, y)
Base.:-(x::AlgebraicArray, y::AlgebraicArray) = x + (-y)
Base.:-(x::AlgebraicArray, y::AbstractArray) = x + (-y)
Base.:-(x::AbstractArray, y::AlgebraicArray) = x + (-y)
Base.:*(a::Number, x::AlgebraicArray) = ScaledArray(a, x)
Base.:*(x::AlgebraicArray, a::Number) = a * x
Base.:\(a::Number, x::AlgebraicArray) = inv(a) * x
Base.:/(x::AlgebraicArray, a::Number) = x * inv(a)

function Base.one(::Type{<:AlgebraicArray{T,2}}, sizes...) where {T}
    sizes = sizes2tuple(Val(2), sizes...)
    sizes[1] == sizes[2] || error("size error")
    return OneArray{T}(sizes[1])
end
Base.:*(x::AlgebraicArray, y::AlgebraicArray) = ProductArray(x, y)
Base.:*(x::AlgebraicArray, y::AbstractArray) = ProductArray(x, y)
Base.:*(x::AbstractArray, y::AlgebraicArray) = ProductArray(x, y)
function Base.:*(x::AlgebraicArray{T,1} where {T},
                 y::AlgebraicArray{T,2} where {T})
    return ProductArray(x, y)
end
function Base.:*(x::AlgebraicArray{T,2} where {T},
                 y::AlgebraicArray{T,1} where {T})
    return ProductArray(x, y)
end
function Base.:*(x::AbstractArray{T,1} where {T},
                 y::AlgebraicArray{T,2} where {T})
    return ProductArray(x, y)
end

################################################################################

export arand
function arand(rng::AbstractRNG, ::Type{AlgebraicArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    if maxdepth ≡ nothing
        maxdepth = 4
    end
    types = [WrappedArray{T,D}, ZeroArray{T,D}]
    if D == 2
        sz = sizes2tuple(Val(2), sizes...)
        if sz[1] == sz[2]
            push!(types, OneArray{T})
        end
    end
    if maxdepth > 0
        append!(types,
                [ScaledArray{T,D}, SumArray{T,D}, ProductArray{T,D},
                 TensorProductArray{T,D}])
        if D ≥ 1
            push!(types, TensorSumArray{T,D})
        end
    end
    return arand(rng, rand(types), sizes...; maxdepth=maxdepth - 1)
end
function arand(T::Type, sizes...; maxdepth=nothing)
    return arand(Random.GLOBAL_RNG, T, sizes...; maxdepth=maxdepth)
end

################################################################################

function Base.:(==)(x::AlgebraicArray, y::AlgebraicArray)
    x ≡ y && return true
    size(x) ≠ size(y) && return false
    # @inbounds for i in eachindex(x)
    #     x[i] ≠ y[i] && return false
    # end
    # return true
    return evaluate(x) == evaluate(y)
end
function Base.:(<)(x::AlgebraicVector, y::AlgebraicVector)
    x ≡ y && return false
    lx = length(x)
    ly = length(y)
    @inbounds for i in 1:min(lx, ly)
        c = cmp(x[i], y[i])
        c ≠ 0 && return c < 0
    end
    return lx < ly
end

################################################################################

export evaluate
evaluate(x::AbstractArray) = x
function arand(rng::AbstractRNG, ::Type{Array{T,D}}, sizes...;
               maxdepth=nothing) where {T,D} end

################################################################################

export WrappedArray
struct WrappedArray{T,D} <: AlgebraicArray{T,D}
    x::AbstractArray{T,D}
end
Base.getindex(x::WrappedArray, inds...) = x.x[inds...]
Base.size(x::WrappedArray) = size(x.x)
evaluate(x::WrappedArray) = x.x
function arand(rng::AbstractRNG, ::Type{WrappedArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    sz = sizes2tuple(Val(D), sizes...)
    if D == 0
        return WrappedArray{T,D}(fill(rand(T), sz))
    end
    return WrappedArray{T,D}(rand(T, sz...))
end

################################################################################

export ZeroArray
struct ZeroArray{T,D} <: AlgebraicArray{T,D}
    sizes::NTuple{D,Int}
end
Base.getindex(x::ZeroArray, inds...) = zero(eltype(x))
Base.size(x::ZeroArray) = x.sizes
evaluate(x::ZeroArray{T,D}) where {T,D} = zeros(T, x.sizes)::AbstractArray{T,D}
function arand(rng::AbstractRNG, ::Type{ZeroArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    return ZeroArray{T,D}(sizes2tuple(Val(D), sizes...))
end

################################################################################

export OneArray
struct OneArray{T} <: AlgebraicArray{T,2}
    size::Int
end
function Base.getindex(x::OneArray, inds...)
    inds = sizes2tuple(Val(2), inds...)
    return inds[1] == inds[2] ? one(eltype(x)) : zero(eltype(x))
end
Base.size(x::OneArray) = (x.size, x.size)
evaluate(x::OneArray{T}) where {T} = Diagonal{T}(ones(T, x.size))
function arand(rng::AbstractRNG, ::Type{OneArray{T}}, sizes...;
               maxdepth=nothing) where {T}
    sizes = sizes2tuple(Val(2), sizes...)
    sizes[1] == sizes[2] || error("size error")
    return OneArray{T}(sizes[1])
end

################################################################################

export ScaledArray
struct ScaledArray{T,D} <: AlgebraicArray{T,D}
    a::Number
    x::AbstractArray{U,D} where {U}
    function ScaledArray{T,D}(a::Number,
                              x::AbstractArray{U,D} where {U}) where {T,D}
        R = typeof(one(typeof(a)) * zero(eltype(x)))
        R <: T || error("type mismatch")
        return new{T,D}(a, x)
    end
    function ScaledArray{T}(a::Number,
                            x::AbstractArray{U,D} where {U}) where {T,D}
        return ScaledArray{T,D}(a, x)
    end
    function ScaledArray(a::Number, x::AbstractArray{U,D} where {U}) where {D}
        T = typeof(one(typeof(a)) * zero(eltype(x)))
        return ScaledArray{T,D}(a, x)
    end
end
Base.getindex(x::ScaledArray, inds...) = x.a * x.x[inds...]
Base.size(x::ScaledArray) = size(x.x)
function evaluate(x::ScaledArray{T,D}) where {T,D}
    return (x.a * evaluate(x.x))::AbstractArray{T,D}
end
function arand(rng::AbstractRNG, ::Type{ScaledArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    a = rand(rng, T)
    x = arand(rng, AlgebraicArray{T,D}, sizes...; maxdepth=maxdepth)
    return ScaledArray{T,D}(a, x)
end

################################################################################

export SumArray
struct SumArray{T,D} <: AlgebraicArray{T,D}
    x::AbstractArray{U,D} where {U}
    y::AbstractArray{U,D} where {U}
    function SumArray{T,D}(x::AbstractArray{U,D} where {U},
                           y::AbstractArray{V,D} where {V}) where {T,D}
        R = typeof(zero(eltype(x)) + zero(eltype(y)))
        R <: T || error("type mismatch")
        size(x) == size(y) || error("size mismatch")
        return new{T,D}(x, y)
    end
    function SumArray{T}(x::AbstractArray{U,D} where {U},
                         y::AbstractArray{V,D} where {V}) where {T,D}
        return SumArray{T,D}(x, y)
    end
    function SumArray(x::AbstractArray{U,D} where {U},
                      y::AbstractArray{V,D} where {V}) where {D}
        T = typeof(zero(eltype(x)) + zero(eltype(y)))
        return SumArray{T,D}(x, y)
    end
end
Base.getindex(x::SumArray, inds...) = x.x[inds...] + x.y[inds...]
Base.size(x::SumArray) = size(x.x)
function evaluate(x::SumArray{T,D}) where {T,D}
    return (evaluate(x.x) + evaluate(x.y))::AbstractArray{T,D}
end
function arand(rng::AbstractRNG, ::Type{SumArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    sz = sizes2tuple(Val(D), sizes...)
    x = arand(rng, AlgebraicArray{T,D}, sz; maxdepth=maxdepth)
    y = arand(rng, AlgebraicArray{T,D}, sz; maxdepth=maxdepth)
    return SumArray{T,D}(x, y)
end

################################################################################

export ProductArray
struct ProductArray{T,D} <: AlgebraicArray{T,D}
    x::AbstractArray{U,D} where {U,D}
    y::AbstractArray{U,D} where {U,D}
    function ProductArray{T,D}(x::AbstractArray, y::AbstractArray) where {T,D}
        R = typeof(one(eltype(x)) * one(eltype(y)))
        R <: T || error("type mismatch")
        (ndims(x) ≥ 1 && ndims(y) ≥ 1) || error("ndims error")
        ndims(x) + ndims(y) == D + 2 || error("ndims mismatch")
        size(x, ndims(x)) == size(y, 1) || error("dimension mismatch")
        return new{T,D}(x, y)
    end
    function ProductArray{T}(x::AbstractArray, y::AbstractArray) where {T}
        D = ndims(x) + ndims(y) - 2
        return ProductArray{T,D}(x, y)
    end
    function ProductArray(x::AbstractArray, y::AbstractArray)
        T = typeof(one(eltype(x)) * one(eltype(y)))
        return ProductArray{T}(x, y)
    end
end
function Base.getindex(x::ProductArray{T,D}, inds...) where {T,D}
    inds = sizes2tuple(Val(D), inds...)
    nd1 = ndims(x.x)
    nd2 = ndims(x.y)
    inds1 = inds[1:(nd1 - 1)]
    inds2 = inds[nd1:D]
    sz′ = size(x.y, 1)
    sz′ == 0 && return zero(T)
    x′ = @view x.x[inds1..., :]
    y′ = @view x.y[:, inds2...]
    return dot(x′, y′)
end
Base.size(x::ProductArray) = (size(x.x)[1:(end - 1)]..., size(x.y)[2:end]...)
function evaluate(x::ProductArray{T,D}) where {T,D}
    sx = size(x.x)
    sy = size(x.y)
    sx′ = (prod(sx[1:(end - 1)]), sx[end])
    sy′ = (sy[1], prod(sy[2:end]))
    return reshape(reshape(evaluate(x.x), sx′) * reshape(evaluate(x.y), sy′),
                   size(x))::AbstractArray{T,D}
end
function arand(rng::AbstractRNG, ::Type{ProductArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    sz = sizes2tuple(Val(D), sizes...)
    nd1 = rand(1:(D + 1))
    nd2 = D + 2 - nd1
    @assert nd1 ≥ 1 && nd2 ≥ 1 && nd1 + nd2 == D + 2
    sz′ = rand(0:((sz === () ? 0 : maximum(sz)) + 2))
    sz1 = (sz[1:(nd1 - 1)]..., sz′)
    sz2 = (sz′, sz[nd1:(nd1 + nd2 - 2)]...)
    @assert length(sz1) == nd1 && length(sz2) == nd2
    x = arand(rng, AlgebraicArray{T,nd1}, sz1; maxdepth=maxdepth)
    y = arand(rng, AlgebraicArray{T,nd2}, sz2; maxdepth=maxdepth)
    return ProductArray{T,D}(x, y)
end

################################################################################

export TensorSumArray
struct TensorSumArray{T,D} <: AlgebraicArray{T,D}
    x::AbstractArray{U,D} where {U,D}
    y::AbstractArray{U,D} where {U,D}
    function TensorSumArray{T,D}(x::AbstractArray{U,D} where {U},
                                 y::AbstractArray{V,D} where {V}) where {T,D}
        D ≥ 1 || error("ndims error")
        R = typeof(zero(eltype(x)) + zero(eltype(y)))
        R <: T || error("type mismatch")
        size(x)[1:(D - 1)] == size(y)[1:(D - 1)] || error("size mismatch")
        return new{T,D}(x, y)
    end
    function TensorSumArray{T}(x::AbstractArray{U,D} where {U},
                               y::AbstractArray{V,D} where {V}) where {T,D}
        return TensorSumArray{T,D}(x, y)
    end
    function TensorSumArray(x::AbstractArray{U,D} where {U},
                            y::AbstractArray{V,D} where {V}) where {D}
        T = typeof(zero(eltype(x)) + zero(eltype(y)))
        return TensorSumArray{T,D}(x, y)
    end
end
function Base.getindex(x::TensorSumArray{T,D}, inds...) where {T,D}
    inds = sizes2tuple(Val(D), inds...)
    sz′ = size(x.x, D)
    if inds[D] ≤ sz′
        return x.x[inds...]
    else
        return x.y[inds[1:(D - 1)]..., inds[D] - sz′]
    end
end
function Base.size(x::TensorSumArray{T,D}) where {T,D}
    sz1 = size(x.x)
    return (sz1[1:(D - 1)]..., size(x.x, D) + size(x.y, D))
end
function evaluate(x::TensorSumArray{T,D}) where {T,D}
    return reshape(vcat(vec(evaluate(x.x)), vec(evaluate(x.y))),
                   size(x))::AbstractArray{T,D}
end
function arand(rng::AbstractRNG, ::Type{TensorSumArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    sz = sizes2tuple(Val(D), sizes...)
    sz′ = rand(0:sz[D])
    sz1 = (sz[1:(D - 1)]..., sz′)
    sz2 = (sz[1:(D - 1)]..., sz[D] - sz′)
    x = arand(rng, AlgebraicArray{T,D}, sz1; maxdepth=maxdepth)
    y = arand(rng, AlgebraicArray{T,D}, sz2; maxdepth=maxdepth)
    return TensorSumArray{T,D}(x, y)
end

################################################################################

export TensorProductArray
struct TensorProductArray{T,D} <: AlgebraicArray{T,D}
    x::AbstractArray{U,D} where {U,D}
    y::AbstractArray{U,D} where {U,D}
    function TensorProductArray{T,D}(x::AbstractArray,
                                     y::AbstractArray) where {T,D}
        R = typeof(one(eltype(x)) * one(eltype(y)))
        R <: T || error("type mismatch")
        ndims(x) + ndims(y) == D || error("ndims mismatch")
        return new{T,D}(x, y)
    end
    function TensorProductArray{T}(x::AbstractArray, y::AbstractArray) where {T}
        D = ndims(x) + ndims(y)
        return TensorProductArray{T,D}(x, y)
    end
    function TensorProductArray(x::AbstractArray, y::AbstractArray)
        T = typeof(one(eltype(x)) * one(eltype(y)))
        return TensorProductArray{T}(x, y)
    end
end
function Base.getindex(x::TensorProductArray{T,D}, inds...) where {T,D}
    inds = sizes2tuple(Val(D), inds...)
    nd1 = ndims(x.x)
    nd2 = ndims(x.y)
    inds1 = inds[1:nd1]
    inds2 = inds[(nd1 + 1):D]
    return x.x[inds1...] * x.y[inds2...]
end
Base.size(x::TensorProductArray) = (size(x.x)..., size(x.y)...)
function evaluate(x::TensorProductArray{T,D}) where {T,D}
    return reshape(vec(evaluate(x.x)) * transpose(vec(evaluate(x.y))),
                   size(x))::AbstractArray{T,D}
end
function arand(rng::AbstractRNG, ::Type{TensorProductArray{T,D}}, sizes...;
               maxdepth=nothing) where {T,D}
    sz = sizes2tuple(Val(D), sizes...)
    nd1 = rand(0:D)
    nd2 = D - nd1
    @assert nd1 ≥ 0 && nd2 ≥ 0 && nd1 + nd2 == D
    sz1 = sz[1:nd1]
    sz2 = sz[(nd1 + 1):D]
    @assert length(sz1) == nd1 && length(sz2) == nd2
    x = arand(rng, AlgebraicArray{T,nd1}, sz1; maxdepth=maxdepth)
    y = arand(rng, AlgebraicArray{T,nd2}, sz2; maxdepth=maxdepth)
    return TensorProductArray{T,D}(x, y)
end

################################################################################

export simplify
simplify(x::AlgebraicArray) = simplify′(x)

simplify′(x::WrappedArray) = x
simplify′(x::ZeroArray) = x
simplify′(x::OneArray) = x

function simplify′(x::ScaledArray)
    iszero(x.a) && return ZeroArray{eltype(x),ndims(x)}(size(x))
    xx = simplify′(x.x)
    xx isa ZeroArray && return xx
    isone(x.a) && return xx
    xx isa ScaledArray && return ScaledArray(x.a * xx.a, xx.x)
    return ScaledArray(x.a, xx)
end

function simplify′(x::SumArray)
    xx = simplify′(x.x)
    xy = simplify′(x.y)
    xy isa ZeroArray && return xx
    xx isa ZeroArray && return xy
    if xx isa ScaledArray && xy isa ScaledArray && xx.x ≡ xy.x
        a = xx.a + xy.a
        iszero(a) && return ZeroArray{eltype(x),ndims(x)}(size(x))
        isone(a) && return xx.x
        return ScaledArray(xx.a + xy.a, xx.x)
    end
    return rassoc(SumArray(xx, xy))
end
function rassoc(x::SumArray)
    if x.x isa SumArray
        return SumArray(x.x.x, rassoc(SumArray(x.x.y, x.y)))
    end
    return x
end

function simplify′(x::ProductArray)
    xx = simplify′(x.x)
    xy = simplify′(x.y)
    if xx isa ZeroArray || xy isa ZeroArray
        return ZeroArray{eltype(x),ndims(x)}(size(x))
    end
    xy isa OneArray && return xx
    xx isa OneArray && return xy
    if xx isa ScaledArray && xy isa ScaledArray
        a = xx.a * xy.a
        if isone(a)
            xx = xx.x
            xy = xy.x
        else
            xx = ScaledArray(a, xx.x)
            xy = xy.x
        end
    end
    return rassoc(ProductArray(xx, xy))
end
function rassoc(x::ProductArray)
    if x.x isa ProductArray && ndims(x.x.y) ≥ 2
        return ProductArray(x.x.x, rassoc(ProductArray(x.x.y, x.y)))
    end
    return x
end

function simplify′(x::TensorSumArray)
    xx = simplify′(x.x)
    xy = simplify′(x.y)
    size(xy, ndims(xy)) == 0 && return xx
    size(xx, ndims(xx)) == 0 && return xy
    if xx isa ZeroArray && xy isa ZeroArray
        return ZeroArray{eltype(x),ndims(x)}(size(x))
    end
    if xx isa ScaledArray && xy isa ScaledArray && xx.a == xy.a
        return ScaledArray(xx.a, TensorSumArray(xx, xy))
    end
    return rassoc(TensorSumArray(xx, xy))
end
function rassoc(x::TensorSumArray)
    if x.x isa TensorSumArray
        return TensorSumArray(x.x.x, rassoc(TensorSumArray(x.x.y, x.y)))
    end
    return x
end

function simplify′(x::TensorProductArray)
    xx = simplify′(x.x)
    xy = simplify′(x.y)
    if xx isa ZeroArray || xy isa ZeroArray
        return ZeroArray{eltype(x),ndims(x)}(size(x))
    end
    if ndims(xy) == 0
        a = xy[]
        isone(a) && return xx
        return ScaledArray(a, xx)
    end
    if ndims(xx) == 0
        a = xx[]
        isone(a) && return xy
        return ScaledArray(a, xy)
    end
    if xx isa ScaledArray && xy isa ScaledArray
        a = xx.a * xy.a
        if isone(a)
            xx = xx.x
            xy = xy.x
        else
            xx = ScaledArray(a, xx.x)
            xy = xy.x
        end
    end
    return rassoc(TensorProductArray(xx, xy))
end
function rassoc(x::TensorProductArray)
    if x.x isa TensorProductArray
        return TensorProductArray(x.x.x, rassoc(TensorProductArray(x.x.y, x.y)))
    end
    return x
end

end


text:
function mc_policy_eval!(mdp, V, action_probability::Function;
                         iterations = 100,
                         Neps = 100,
                         verbose = true,
                         max_steps_per_ep = 1000,
                         opt = ADAM(1e-3),
                         batchsize = 32,
                         shuffle = true,
                         rng::AbstractRNG = Random.GLOBAL_RNG,
                         a_and_p = action_and_probability,
                         logger = TBLogger("log/", tb_increment),
                         log_freq = 1,
                         exploration_policy = nothing)
    log = !isnothing(logger)
    loss = (x, y) -> Flux.mse(V(x), y)
    best_avgR = 0
    for iter in 1:iterations
        verbose && println("iteration: ", iter)
        # Setup the exploration policy. If none is provided, use the current policy
        isnothing(exploration_policy) && (exploration_policy = ISPolicy(V, mdp, action_probability, rng))

        # Sample episodes. S is the states, G is the per-state return (with IS weights), R is the undiscounted return for each episode
        S, G, avgR = sample_episodes(mdp, exploration_policy, Neps,
                                 max_steps = max_steps_per_ep, a_and_p = a_and_p)

        # If the returns are higher (i.e. more failures found) then save the model
        log && avgR > best_avgR && @save string(logger.logdir, "/best_model.bson") V

        # Load the data into a DataLoader and train
        data = Flux.Data.DataLoader(S, G, batchsize=batchsize, shuffle = shuffle)
        Flux.train!(loss, Flux.params(V), data, opt)

        # If we are logging then record the loss and the avg return
        (!log || iter % log_freq != 0) && continue
        with_logger(logger) do
          @info "Training" loss=loss(S, G) undiscounted_return = avgR
          @save string(logger.logdir, "/last_model.bson") V
        end
    end
end


text:
# Load reference values for the simulation
PDE_refs = load(joinpath(ODINN.root_dir, "data/PDE_refs.jld2"))
train_settings = (ADAM(0.03), 10) # optimizer, epochs
iceflow_trained, UA = @time train_iceflow_UDE(gdirs_climate, tspan, train_settings, PDE_refs)

# Verify neural PDE is learning
@test loss_history[1] > 3*loss_history[end]


text:
using BinDeps

@BinDeps.setup

pfnetname = "pfnet-1.3.3"
libpfnet = library_dependency("libpfnet")
prefix=joinpath(BinDeps.depsdir(libpfnet),"usr")
srcdir=joinpath(BinDeps.depsdir(libpfnet),"src")
downloads=joinpath(BinDeps.depsdir(libpfnet),"downloads")

# Linux
if is_linux()
    provides(SimpleBuild,
             (@build_steps begin
                ChangeDirectory(downloads)
                @build_steps begin
                  `gunzip -kf $pfnetname.tar.gz`
                  `tar -xvf $pfnetname.tar`
                  `rm $pfnetname.tar`
                  `mkdir -p $srcdir`
                  `rm -rf $srcdir/$pfnetname`
                  `mv $pfnetname $srcdir`
                end
              end,
              @build_steps begin
                ChangeDirectory(joinpath(srcdir, pfnetname))
                @build_steps begin
                  `./configure --prefix=$prefix`
                  `make clean`
                  `make uninstall`
                  `make`
                  `make check`
                  `make install`
                end
              end),
             libpfnet)
end

# Apple
if is_apple()
    provides(SimpleBuild,
             (@build_steps begin
                ChangeDirectory(downloads)
                @build_steps begin
                  `gunzip -kf $pfnetname.tar.gz`
                  `tar -xvf $pfnetname.tar`
                  `rm $pfnetname.tar`
                  `mkdir -p $srcdir`
                  `rm -rf $srcdir/$pfnetname`
                  `mv $pfnetname $srcdir`
                end
              end,
              @build_steps begin
                ChangeDirectory(joinpath(srcdir, pfnetname))
                @build_steps begin
                  `./configure --prefix=$prefix LDFLAGS="-Wl,-install_name,@loader_path/libpfnet.dylib"`
                  `make clean`
                  `make uninstall`
                  `make`
                  `make check`
                  `make install`
                end
              end),
             libpfnet)
end

# Windows
if is_windows()

end

@BinDeps.install Dict(:libpfnet => :libpfnet)


text:
@testset "Utils" begin
    @test Quadtrees.distance²((3, 4), (0, 0)) == 25
    @test Quadtrees.distance((3, 4), (0, 0)) == 5

    for x in 45:55
        for y in 45:55
            should_be_true = (x, y) in Rect(50, 50, 10, 10)
            should_be_false = (x, y) in Rect(0, 0, 10, 10)
            @test should_be_true == true
            @test should_be_false == false
        end
    end

    @test_throws ErrorException [0, 0] in Rect(0, 0, 10, 10)

    Quadtrees.position(::Nothing) = nothing

    @test_throws AssertionError nothing in Rect(0, 0, 10, 10)

    for x in 45:65
        for y in 45:55
            @test Rect(50, 50, x, y) in Rect(50, 55, 10, 20)
            @test !(Rect(0, 0, 1, 1) in Rect(5, 5, random(5), random(5)))
        end
    end
end

@testset "Intersections" begin
    @init "Intersections test" 800 800
    @layout aside(:v, 230)
    @create @options begin
        use_rectangle = false
        radius = (1:400, 100)
        spacing = (10, 1:50)
        parameters = @nolabel (
            rows = (1:100, 15),
            cols = (1:100, 15),
            margin = (1:50, 20)
        )
    end

    mouse = Ref{Point}(O)

    @use function update()
        background("black")
        origin()
        w, h = @width, @height

        rows, cols, margin = Int.(getindex.(parameters))
        tiles = Tiler(@width, @height, rows, cols; margin)
        tw = tiles.tilewidth - spacing[]
        th = tiles.tileheight - spacing[]

        rects = [Rect(pos.x, pos.y, tw, th) for pos in first.(tiles)]
        mousepos = mouse[] - (w/2, h/2)

        if use_rectangle[]
            range = Rect(mousepos..., radius[], radius[])
        else
            range = Quadtrees.Circle(mousepos..., radius[])
        end

        sethue("yellowgreen")
        draw(range, :stroke)

        for rect in rects
            sethue("white")
            rect in range && sethue("yellowgreen")
            draw(rect, :stroke)
        end
    end

    @use function mousemove(event)
        mouse[] = event.pos
    end

    start(; async = true)
    destroy(@window)
end


text:
# Missing/unknown value
M = -1

#
# data
#
#  _ is coded as unknown
#  0..8: known number of neighbours
#

# The first 10 examples (0..9) are from Gecode/examples/minesweeper.cc
# http://www.gecode.org/gecode-doc-latest/minesweeper_8cc-source.html
# """
# The instances are taken from
#   http://www.janko.at/Raetsel/Minesweeper/index.htm
# """

all_minesweeper_problems = Dict(

# Problem from Gecode/examples/minesweeper.cc  problem 0
#
# Solution:
#  1 0 0 0 0 1
#  0 1 0 1 1 0
#  0 0 0 0 1 0
#  0 0 0 0 1 0
#  0 1 1 1 0 0
#  1 0 0 0 1 1
:0 =>
    [[M,M,2,M,3,M],
     [2,M,M,M,M,M],
     [M,M,2,4,M,3],
     [1,M,3,4,M,M],
     [M,M,M,M,M,3],
     [M,3,M,3,M,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 1
:1 =>
    [[M,2,M,2,1,1,M,M],
     [M,M,4,M,2,M,M,2],
     [2,M,M,2,M,M,3,M],
     [2,M,2,2,M,3,M,3],
     [M,M,1,M,M,M,4,M],
     [1,M,M,M,2,M,M,3],
     [M,2,M,2,2,M,3,M],
     [1,M,1,M,M,1,M,1]],

# Problem from Gecode/examples/minesweeper.cc  problem 2
:2 =>
    [[1,M,M,2,M,2,M,2,M,M],
     [M,3,2,M,M,M,4,M,M,1],
     [M,M,M,1,3,M,M,M,4,M],
     [3,M,1,M,M,M,3,M,M,M],
     [M,2,1,M,1,M,M,3,M,2],
     [M,3,M,2,M,M,2,M,1,M],
     [2,M,M,3,2,M,M,2,M,M],
     [M,3,M,M,M,3,2,M,M,3],
     [M,M,3,M,3,3,M,M,M,M],
     [M,2,M,2,M,M,M,2,2,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 3
:3 =>
    [[2,M,M,M,3,M,1,M],
     [M,5,M,4,M,M,M,1],
     [M,M,5,M,M,4,M,M],
     [2,M,M,M,4,M,5,M],
     [M,2,M,4,M,M,M,2],
     [M,M,5,M,M,4,M,M],
     [2,M,M,M,5,M,4,M],
     [M,3,M,3,M,M,M,2]],

# Problem from Gecode/examples/minesweeper.cc  problem 4
  :4 =>
    [[0,M,0,M,1,M,M,1,1,M],
     [1,M,2,M,2,M,2,2,M,M],
     [M,M,M,M,M,M,2,M,M,2],
     [M,2,3,M,1,1,M,M,M,M],
     [0,M,M,M,M,M,M,2,M,1],
     [M,M,M,2,2,M,1,M,M,M],
     [M,M,M,M,M,3,M,3,2,M],
     [M,5,M,2,M,M,M,3,M,1],
     [M,3,M,1,M,M,3,M,M,M],
     [M,2,M,M,M,1,2,M,M,0]],

# Problem from Gecode/examples/minesweeper.cc  problem 5
:5 =>
    [[M,2,1,M,2,M,2,M,M,M],
     [M,4,M,M,3,M,M,M,5,3],
     [M,M,M,4,M,4,4,M,M,3],
     [4,M,4,M,M,5,M,6,M,M],
     [M,M,4,5,M,M,M,M,5,4],
     [3,4,M,M,M,M,5,5,M,M],
     [M,M,4,M,4,M,M,5,M,5],
     [2,M,M,3,3,M,6,M,M,M],
     [3,6,M,M,M,3,M,M,4,M],
     [M,M,M,4,M,2,M,2,1,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 6
:6 =>
    [[M,3,2,M,M,1,M,M],
     [M,M,M,M,1,M,M,3],
     [3,M,M,2,M,M,M,4],
     [M,5,M,M,M,5,M,M],
     [M,M,6,M,M,M,5,M],
     [3,M,M,M,5,M,M,4],
     [2,M,M,5,M,M,M,M],
     [M,M,2,M,M,3,4,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 7
:7 =>
    [[M,1,M,M,M,M,M,3,M],
     [M,M,M,3,4,3,M,M,M],
     [2,4,4,M,M,M,4,4,3],
     [M,M,M,4,M,4,M,M,M],
     [M,4,M,4,M,3,M,6,M],
     [M,M,M,4,M,3,M,M,M],
     [1,2,3,M,M,M,1,3,3],
     [M,M,M,3,2,2,M,M,M],
     [M,2,M,M,M,M,M,3,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 8
:8 =>
    [[M,M,M,M,M,M,M],
     [M,2,3,4,3,5,M],
     [M,1,M,M,M,3,M],
     [M,M,M,5,M,M,M],
     [M,1,M,M,M,3,M],
     [M,1,2,2,3,4,M],
     [M,M,M,M,M,M,M]],

# Problem from Gecode/examples/minesweeper.cc  problem 9
:9 =>
    [[2,M,M,M,2,M,M,M,2],
     [M,4,M,4,M,3,M,4,M],
     [M,M,4,M,M,M,1,M,M],
     [M,4,M,3,M,3,M,4,M],
     [2,M,M,M,M,M,M,M,2],
     [M,5,M,4,M,5,M,4,M],
     [M,M,3,M,M,M,3,M,M],
     [M,4,M,3,M,5,M,6,M],
     [2,M,M,M,1,M,M,M,2]],

# From "Some Minesweeper Configurations",page 2
:10 =>
    [[M,M,M,M,M,M],
      [M,2,2,2,2,M],
      [M,2,0,0,2,M],
      [M,2,0,0,2,M],
      [M,2,2,2,2,M],
      [M,M,M,M,M,M]],

# From "Some Minesweeper Configurations",page 3
# 4 solutions
:11 =>
    [[2,3,M,2,2,M,2,1],
      [M,M,4,M,M,4,M,2],
      [M,M,M,M,M,M,4,M],
      [M,5,M,6,M,M,M,2],
      [2,M,M,M,5,5,M,2],
      [1,3,4,M,M,M,4,M],
      [0,1,M,4,M,M,M,3],
      [0,1,2,M,2,3,M,2]],

# Richard Kaye: How Complicated is Minesweeper?
# http://web.mat.bham.ac.uk/R.W.Kaye/minesw/ASE2003.pdf
#
# A Wire,page 33
# 2 solutions
#
:12 =>
    [[M,0,0,0,0,0,0,0,0,0,0,0,0,M],
      [M,1,1,1,1,1,1,1,1,1,1,1,1,M],
      [M,M,1,M,M,1,M,M,1,M,M,1,M,M],
      [M,1,1,1,1,1,1,1,1,1,1,1,1,M],
      [M,0,0,0,0,0,0,0,0,0,0,0,0,M]],

# Richard Kaye: How Complicated is Minesweeper?
# http://web.mat.bham.ac.uk/R.W.Kaye/minesw/ASE2003.pdf
# A splitter,page 35
# Many solutions...
#
:13 =>
     [[M,M,M,0,M,M,M,0,M,M,M],
      [M,M,M,0,1,M,1,0,M,M,M],
      [M,M,M,0,1,M,1,0,M,M,M],
      [0,0,0,0,1,1,1,0,0,0,0],
      [M,1,1,1,1,M,1,1,1,1,M],
      [M,M,M,1,M,2,M,1,M,M,M],
      [M,1,1,1,1,M,1,1,1,1,M],
      [0,0,0,0,1,1,1,0,0,0,0],
      [M,M,M,0,1,M,1,0,M,M,M],
      [M,M,M,0,1,M,1,0,M,M,M],
      [M,M,M,0,M,M,M,0,M,M,M]],

# Oleg German,Evgeny Lakshtanov: "Minesweeper" without a computer
# http://arxiv.org/abs/0806.3480, page 4
:14 =>
     [[M,1,M,1,M,1],
      [2,M,2,M,1,M],
      [M,3,M,2,M,1],
      [1,M,3,M,2,M],
      [M,1,M,2,M,1]],

#
# From http://stephenlombardi.com/minesweeper/minesweeper.lisp
#
:15 =>
    [[0,0,0,0,1,M,M,M,M],
     [0,0,0,0,1,M,M,3,M],
     [0,0,1,1,2,M,M,M,M],
     [0,0,1,M,M,M,M,1,M],
     [0,0,1,2,M,3,M,M,M],
     [0,0,0,1,M,M,M,M,M],
     [0,0,1,2,2,1,1,1,1],
     [0,0,1,M,1,0,0,0,0],
     [0,0,1,M,1,0,0,0,0]]

)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@testset verbose = true "statistics" begin
    @testset "analytic_histogram" begin
        # check against uniform distribution
        d = Distributions.Uniform(0, 1)
        G = LinRange(0.0, 1.0, 10)
        res = analytic_histogram(d, G)

        @test res[1] ≈ step(G) / 2
        for i in 2:(length(G) - 1)
            @test res[i] ≈ step(G)
        end
        @test res[end] ≈ step(G) / 2
    end
end


text:
import LazySets.Approximations.overapproximate

"""
    overapproximate(X::LazySet{N}, dirs::Vector{Vector{N}}) where {N}

 Compute the overapproximation of the lazy set `X` by evaluating the support values in the direction of `dirs` and adding the resulting constraints.

 ### Inputs
 - `X` -- lazy set
 - `dir` -- vector of direction vectors

 ### Output
 Over-approximation of the lazy set `X` in the directions of `dir`
"""
function LazySets.Approximations.overapproximate(X::LazySet{N}, dirs::Vector{Vector{N}}) where {N}
    halfspaces = Vector{LinearConstraint{N}}()
    sizehint!(halfspaces, length(dirs))
    H = HPolytope(halfspaces)
    for d in dirs
        addconstraint!(H, LinearConstraint(d, ρ(d, X)))
    end
    return H
end

"""
     refine(P::LazySet; eps=1e-5, ndirs=200, nsamples=2000)

 Compute an concrete overapproximation of the lazy set `P`.
 The accuracy depends on the hyperparameter:
 2D-set: eps       => ϵ-close approximation
 3D-set: n_dirs    => #directions for evaluation the support value
 nD-set: n_samples => #directions for evaluation the support value

 ### Inputs
 - `P` -- lazy set
 - `eps`-- (optional, default=1e-5) define error tolarance of ϵ-close overapproximation of 2D set
 - `n_dirs`-- (optional, default=200) number of direction vectors evaluated for 3D set
 - `n_samples`-- (optional, default=2000) number of direction vectors evaluated for nD set with n>3

 ### Output
 Overapproximation of lazy set `P`.
"""
function refine(P::LazySet; eps=1e-5, n_dirs=200, n_samples=2000)
    n = LazySets.dim(P)
    if n == 2
        return refine_2D(P, eps)
    elseif n == 3
        return refine_3D(P, n_dirs)
    else
        P = refine_nD(P, n_samples)
        # NOTE: Not sure about remove_redundant_constraints! functionalities!
        remove_redundant_constraints!(P)
        return P
    end
end

function refine_nD(X, nsamples)
    n = LazySets.dim(X)
    dirs = Vector{Vector{Float64}}(undef, nsamples)
    LazySets._sample_unit_nsphere_muller!(dirs, n, nsamples)
    return overapproximate(X, dirs)
end

function refine_3D(X, n_dirs)
    dirs = Approximations.SphericalDirections(n_dirs, n_dirs)
    return overapproximate(X, dirs)
end

function refine_2D(X, err_approx)
    return overapproximate(X, err_approx)
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
abstract type AbstractBoundaryCondition end

"""
    FluidBC(momentum    = Impenetrable(NoSlip())
            temperature = Insulating())

The standard boundary condition for CNSEModel. The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct FluidBC{M, T} <: AbstractBoundaryCondition
    ρu::M = Impenetrable(NoSlip())
    ρθ::T = Insulating()
end

abstract type StateBC end
abstract type MomentumBC <: StateBC end
abstract type MomentumDragBC <: StateBC end
abstract type TemperatureBC <: StateBC end

(bc::StateBC)(state, aux, t) = bc.flux(bc.params, state, aux, t)

"""
    Impenetrable(drag::MomentumDragBC) :: MomentumBC

Defines an impenetrable wall model for momentum. This implies:
  - no flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Impenetrable{D <: MomentumDragBC} <: MomentumBC
    drag::D
end

"""
    Penetrable(drag::MomentumDragBC) :: MomentumBC

Defines an penetrable wall model for momentum. This implies:
  - no constraint on flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Penetrable{D <: MomentumDragBC} <: MomentumBC
    drag::D
end

"""
    NoSlip() :: MomentumDragBC

Zero momentum at the boundary.
"""
struct NoSlip <: MomentumDragBC end

"""
    FreeSlip() :: MomentumDragBC

No surface drag on momentum parallel to the boundary.
"""
struct FreeSlip <: MomentumDragBC end

"""
    MomentumFlux(stress) :: MomentumDragBC

Applies the specified kinematic stress on momentum normal to the boundary.
Prescribe the net inward kinematic stress across the boundary by `stress`,
a function with signature `stress(problem, state, aux, t)`, returning the flux (in m²/s²).
"""
Base.@kwdef struct MomentumFlux{𝒯, 𝒫} <: MomentumDragBC
    flux::𝒯 = nothing
    params::𝒫 = nothing
end

"""
    Insulating() :: TemperatureBC

No temperature flux across the boundary
"""
struct Insulating <: TemperatureBC end

"""
    TemperatureFlux(flux) :: TemperatureBC

Prescribe the net inward temperature flux across the boundary by `flux`,
a function with signature `flux(problem, state, aux, t)`, returning the flux (in m⋅K/s).
"""
struct TemperatureFlux{T} <: TemperatureBC
    flux::T

    function TemperatureFlux(flux::T = nothing) where {T}
        new{T}(flux)
    end
end

abstract type AbstractCouplerBoundary end

export PenaltyNumFluxDiffusive, CoupledPrimaryBoundary, CoupledSecondaryBoundary

"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end
- PenaltyNumFluxDiffusive: additional NumericalFluxSecondOrder to account for an additional penalty term
"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end

# ## 2. CoupledPrimaryBoundary
"""
struct CoupledPrimaryBoundary <: AbstractCouplerBoundary end
# # compute flux based on opposite face
# # also need to accumulate net flux across boundary
"""
struct CoupledPrimaryBoundary <: AbstractCouplerBoundary end

## 3. CoupledSecondaryBoundary
"""
struct CoupledSecondaryBoundary <: AbstractCouplerBoundary end
# # use prescribed flux computed in primary
"""
struct CoupledSecondaryBoundary <: AbstractCouplerBoundary end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using TaylorModels

include("TM1.jl")
include("TMN.jl")
include("shrink-wrapping.jl")
include("validated_integ.jl")


text:
module CodeGenerator
import Base:show

# "文法"类型定义
Terminal = Symbol
Nonterminal = Symbol
Token = Union{Terminal,Nonterminal}
Statement = Vector{Token}
Production = Dict{Token,Set{Statement}}

# 状态机类型定义
State = Symbol
StateSet = Set{State}
SymbolTabel = Set{Symbol}
Transition{T} = Dict{State,Dict{Symbol,T}}

include("utils.jl")
include("grammar.jl")
include("generator.jl")
include("NFA.jl")

end # module

max_stars_repo_name: yeruoforever/SimpleGrammar

max_issues_repo_name: yeruoforever/SimpleGrammar

max_forks_repo_name: yeruoforever/SimpleGrammar


text:
function dtftau!(champu  :: Array{ComplexF64, 2},
                 champv  :: Array{ComplexF64, 2},
                 m       :: MicMac,
                 t       :: Float64,
                 fft_u   :: Vector{ComplexF64},
                 fft_v   :: Vector{ComplexF64})

    sigma = 1

    m.u .= exp.(1im .* t .* m.A1)
    m.v .= -1im .* m.llambda .* m.A2 .* conj.(m.u)

    m.ut .= (m.u .* fft_u) .* m.matr
    m.vt .= (m.u .* fft_v) .* m.matr

    ifft!(m.ut, 1)
    ifft!(m.vt, 1)

    m.z .= ( m.ut .+ conj.(m.vt)) ./ 2

    fft_u .*= m.u .* (1im .* m.A1)
    fft_v .*= m.u .* (1im .* m.A1)

    ifft!(fft_u)
    ifft!(fft_v)

    m.ut .= fft_u .* m.matr
    m.vt .= fft_v .* m.matr

    m.ut .= (m.ut .+ conj.(m.vt)) ./ 2

    m.vt .= 2 .* abs.(m.z).^2 .* m.ut .+ m.z.^2 .* conj.(m.ut)

    m.ut .= m.conjmatr .* m.vt
    m.vt .= m.conjmatr .* conj.(m.vt)

    fft!(m.ut, 1)
    fft!(m.vt, 1)

    m.ut .*= m.v
    m.vt .*= m.v

    champu .+= transpose(m.ut)
    champv .+= transpose(m.vt)

    m.z  .= abs.(m.z).^2 .* m.z

    m.ut .= m.conjmatr .* m.z
    m.vt .= m.conjmatr .* conj.(m.z)

    fft!(m.ut, 1)
    fft!(m.vt, 1)

    m.ut .*= (m.v .* (-1im .* m.A1))
    m.vt .*= (m.v .* (-1im .* m.A1))

    champu .+= transpose(m.ut)
    champv .+= transpose(m.vt)

end


text:
module Euphony

using Base
import Base: <, <=, ==

abstract Sequencable

include("delegate.jl")

include("event.jl")
export Event, assoc, dissoc

include("sequence.jl")
export Seq, HSeq, VSeq

include("chord.jl")
export Chord

include("scale.jl")
export Scale

include("tuning.jl")
export Tuning

include("frequency.jl")
export to_freq

include("duration.jl")
export Tempo, to_second

end


text:
"""
    Timer

A sub-module of ModelBaseEcon, although it can be used independently.
Provides functionality for measuring the aggregate time spent
in individual lines or blocks of code.

!!! tip

    The implementation here is quick-and-dirty and is intended for getting
    a rough idea of where the bottlenecks are. It is best used for timing
    blocks of code that are passed through relatively few times and each pass
    takes a relatively long time. In the opposite extreme case (fast code that
    is called many times), the current implementation of `@timer` might add
    extreme overhead.

# Contents

  * [`inittimer`](@ref)       - Enable collection of timer data.
  * [`stoptimer`](@ref)       - Disable collection of timer data.
  * [`printtimer`](@ref)      - Display timer data.
  * [`@timer`](@ref)          - Measure the runtime taken by the given code.

# Example
```jldoctest
julia> true
[...]
```
"""
module Timer

using Printf

export @timer, inittimer, printtimer, stoptimer

"""
    timerData

Stores timing data.

!!! note

    For internal use. Do not modify directly.

If equal to `nothing`, timing is disabled.
Otherwise, contains a "database" of timimng data in the
form of a Dict.
"""
global timerData = nothing

"""
    inittimer()

Enable the collection of timing data. Existing timing data is lost.
By default, collection of timing data is disabled.

See also: [`stoptimer`](@ref), [`@timer`](@ref), [`printtimer`](@ref)
"""
function inittimer()
    global timerData = Dict()
    return nothing
end

"""
    stoptimer()

Disable the collection of timing data. Existing data is lost.
By default, collection of timing data is disabled.

See also: [`inittimer`](@ref), [`@timer`](@ref), [`printtimer`](@ref)
"""
function stoptimer()
    global timerData = nothing
    return nothing
end

"""
    printtimer(io::IO=Base.stdout)

Display timing data.

Timing data is displayed in a table with each row containing the
number of calls, total time in seconds, and the source line or block tag.
Rows are sorted in order of decreasing total time.

See also: [`@timer`](@ref)
"""
function printtimer(io::IO = Base.stdout)
    global timerData
    if timerData !== nothing
        @printf(io, "%10s  %10s  %s\n", "Calls", "Seconds", "Source")
        items = collect(timerData)
        items = sort(items, lt = (l, r)->r[2][:seconds] < l[2][:seconds])
        for (src, val) ∈ items
            @printf(io, "% 10d  %10.3f  %s\n", val[:calls], val[:seconds], src)
        end
    end
    return nothing
end

"""
    @timer(code)
    @timer(tag::String, code)

Measure the number of calls and the total time taken by the given code.

If a `tag` string is not provided, one is generated from the source file and line.
The return value of this macro call is the return value of the code.

!!! warning
    Important limitation is that the code must not contain a `return`, `break`,
    `continue`, or any other jump out of it. If it does, the program would run
    correctly, but the timing data collected would be incorrect.

See also: [`inittimer`](@ref), [`stoptimer`](@ref), [`printtimer`](@ref)

# Example

```jldoctest
julia> inittimer()

julia> @timer Base.sleep(1.0)

julia> printtimer()
[...]
```
"""
macro timer(args...)
    if length(args) == 1
        name = QuoteNode(__source__)
        code = args[1]
    elseif length(args) == 2
        name = string(args[1])
        code = args[2]
    else
        error("Too many arguments")
    end
    return quote
        global timerData
        local val
        if timerData === nothing
            val = $(esc(code))
        else
            local t = time()
            val = $(esc(code))
            t = time() - t
            local td = get!(timerData, $name, Dict(:calls => 0, :seconds => 0.0))
            td[:calls] += 1
            td[:seconds] += t
            timerData[$name] = td
        end
        val
    end
end

end # module


text:
function sum_n_below(factor, max)
    # Return the sum of all the natural numbers of factor below max
    max_nb = div(max-1,factor)
    ((factor / 2 * max_nb * (max_nb + 1)))
end

function n_below(n)
    # Return all the natural number of 3 and 5 multipliers
    return sum_n_below(3, n) + sum_n_below(5, n) - sum_n_below(3*5, n)
end

@assert n_below(10) == 23
println(n_below(1000))

max_stars_repo_name: ddugue/project-euler

max_issues_repo_name: ddugue/project-euler

max_forks_repo_name: ddugue/project-euler


text:
module TestEnsembleMethods

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
nfcp = NumericFeatureClassification()

using FactCheck

importall Orchestra.Transformers.EnsembleMethods
import Orchestra.Transformers.DecisionTreeWrapper: fit!, transform!
import Orchestra.Transformers.DecisionTreeWrapper: PrunedTree
import Orchestra.Transformers.DecisionTreeWrapper: RandomForest
import Orchestra.Transformers.DecisionTreeWrapper: DecisionStumpAdaboost

facts("Ensemble learners") do
  context("VoteEnsemble predicts according to majority") do
    always_a_options = { :label => "a" }
    always_b_options = { :label => "b" }
    learner = VoteEnsemble({
      :learners => [
        AlwaysSameLabelLearner(always_a_options),
        AlwaysSameLabelLearner(always_a_options),
        AlwaysSameLabelLearner(always_b_options)
      ]
    })
    fit!(learner, nfcp.train_instances, nfcp.train_labels)
    predictions = transform!(learner, nfcp.test_instances)
    expected_predictions = fill("a", size(nfcp.test_instances, 1))

    @fact predictions => expected_predictions
  end

  context("StackEnsemble predicts with combined learners") do
    # Fix random seed, due to stochasticity in stacker.
    srand(2)

    always_a_options = { :label => "a" }
    learner = StackEnsemble({
      :learners => [
        AlwaysSameLabelLearner(always_a_options),
        AlwaysSameLabelLearner(always_a_options),
        PerfectScoreLearner()
      ],
      :keep_original_features => true
    })
    fit!(learner, nfcp.train_instances, nfcp.train_labels)
    predictions = transform!(learner, nfcp.test_instances)
    unexpected_predictions = fill("a", size(nfcp.test_instances, 1))

    @fact predictions => not(unexpected_predictions)
  end

  context("BestLearner picks the best learner") do
    always_a_options = { :label => "a" }
    always_b_options = { :label => "b" }
    learner = BestLearner({
      :learners => [
        AlwaysSameLabelLearner(always_a_options),
        PerfectScoreLearner(),
        AlwaysSameLabelLearner(always_b_options)
      ]
    })
    fit!(learner, nfcp.train_instances, nfcp.train_labels)

    @fact learner.model[:best_learner_index] => 2
  end

  context("BestLearner conducts grid search") do
    learner = BestLearner({
      :learners => [PrunedTree(), DecisionStumpAdaboost(), RandomForest()],
      :learner_options_grid => [
        {
          :impl_options => {
            :purity_threshold => [0.5, 1.0]
          }
        },
        Dict(),
        {
          :impl_options => {
            :num_trees => [5, 10, 20],
            :partial_sampling => [0.5, 0.7]
          }
        }
      ]
    })
    fit!(learner, nfcp.train_instances, nfcp.train_labels)

    @fact length(learner.model[:learners]) => 8
  end
end

end # module


text:
module Funs

using DifferentialForms: Forms, unit
using LinearAlgebra
using SparseArrays
using StaticArrays

using ..Defs
using ..Manifolds
using ..SparseOps

export Fun
"""
Function (aka Cochain)
"""
struct Fun{D,P,R,C,S,T}         # <: AbstractVector{T}
    manifold::Manifold{D,C,S}
    values::IDVector{R,T}

    function Fun{D,P,R,C,S,T}(manifold::Manifold{D,C,S},
                              values::IDVector{R,T}) where {D,P,R,C,S,T}
        D::Int
        @assert D ≥ 0
        P::PrimalDual
        R::Int
        @assert 0 ≤ R ≤ D
        fun = new{D,P,R,C,S,T}(manifold, values)
        #TODO @assert invariant(fun)
        return fun
    end
end

function Fun{D,P,R,C,S}(manifold::Manifold{D,C,S},
                        values::IDVector{R,T}) where {D,P,R,C,S,T}
    return Fun{D,P,R,C,S,T}(manifold, values)
end
function Fun{D,P,R}(manifold::Manifold{D,C,S},
                    values::IDVector{R,T}) where {D,P,R,C,S,T}
    return Fun{D,P,R,C,S,T}(manifold, values)
end

function Base.show(io::IO, fun::Fun{D,P,R,C,S,T}) where {D,P,R,C,S,T}
    println(io)
    println(io, "Fun{$D,$P,$R,$S,$T}(")
    println(io, "    manifold=\"$(fun.manifold.name)\"")
    println(io, "    values=$(fun.values)")
    return print(io, ")")
end

function Defs.invariant(fun::Fun{D,P,R,C,S,T})::Bool where {D,P,R,C,S,T}
    D::Int
    P::PrimalDual
    R::Int
    0 ≤ R ≤ D || return false
    invariant(fun.manifold) || return false
    length(fun.values) == nsimplices(fun.manifold, R) || return false
    return true
end

# Comparison

function Base.:(==)(f::F, g::F) where {F<:Fun}
    @assert f.manifold == g.manifold
    f ≡ g && return true
    return f.values == g.values
end
function Base.:(<)(f::F, g::F) where {F<:Fun}
    @assert f.manifold == g.manifold
    return f.values < g.values
end
function Base.isequal(f::F, g::F) where {F<:Fun}
    @assert f.manifold == g.manifold
    return isequal(f.values, g.values)
end
function Base.hash(f::Fun, h::UInt)
    return hash(0xfc734743, hash(f.manifold, hash(f.values, h)))
end

# Random functions
function Base.rand(::Type{Fun{D,P,R,C,S,T}},
                   manifold::Manifold{D,C,S}) where {D,P,R,C,S,T}
    return Fun{D,P,R,C,S}(manifold,
                          IDVector{R}(rand(T, nsimplices(manifold, R))))
end

# Functions are a collection

Base.eltype(::Type{<:Fun{<:Any,<:Any,<:Any,<:Any,<:Any,T}}) where {T} = T
Base.eltype(f::Fun) = eltype(typeof(f))
Base.isempty(f::Fun) = isempty(f.values)
Base.iterate(f::Fun, state...) = iterate(f.values, state...)
Base.length(f::Fun) = length(f.values)

# function Base.map(op, f::Fun{D,P,R}, gs::Fun{D,P,R}...) where {D,P,R}
#     @assert all(f.manifold == g.manifold for g ∈ gs)
#     return Fun{D,P,R}(f.manifold,
#                         map(op, f.values, (g.values for g ∈ gs)...))
#     # r1 = map(op, f.values[1], (g.values[1] for g ∈ gs)...))
#     # T = typeof(r1)
#     # Fun{D, P, R, S, T}(f.manifold, map(op, f.values, (g.values for g ∈ gs)...))
# end
function Base.map(op, f::Fun{D,P,R,C,S},
                  gs::Fun{D,P,R,C,S}...) where {D,P,R,C,S}
    @assert all(g.manifold == f.manifold for g in gs)
    U = typeof(op(zero(eltype(f)), (zero(eltype(g)) for g in gs)...))
    return Fun{D,P,R,C,S,U}(f.manifold,
                            map(op, f.values, (g.values for g in gs)...))
end
function Base.reduce(op, f::Fun{D,P,R,C,S}, gs::Fun{D,P,R,C,S}...;
                     kw...) where {D,P,R,C,S}
    @assert all(g.manifold == f.manifold for g in gs)
    return reduce(op, f.values, (g.values for g in gs)...; kw...)
end

# Functions are an abstract vector

Base.IndexStyle(::Type{<:Fun}) = IndexStyle(Vector)
Base.axes(f::Fun) = axes(f.values)
Base.axes(f::Fun, dir) = axes(f.values, dir)
Base.eachindex(f::Fun) = eachindex(f.values)
Base.@propagate_inbounds function Base.getindex(f::Fun, inds...)
    return getindex(f.values, inds...)
end
Base.ndims(::Fun) = 1
Base.size(f::Fun) = size(f.values)
Base.size(f::Fun, dims) = size(f.values, dims)
Base.stride(f::Fun, k) = stride(f.values, k)
Base.strides(f::Fun) = strides(f.values)

# Functions are a vector space

function Base.zero(::Type{Fun{D,P,R,C,S,T}},
                   manifold::Manifold{D,C,S}) where {D,P,R,C,S,T}
    nelts = nsimplices(manifold, R)
    return Fun{D,P,R}(manifold, IDVector{R}(zeros(T, nelts)))
end
function Base.zero(::Type{Fun{D,P,R,C,S}},
                   manifold::Manifold{D,C,S}) where {D,P,R,C,S,T}
    return zero(Fun{D,P,R,C,S,Bool}, manifold)
end
Base.zero(f::Fun) = zero(typeof(f), f.manifold)
Base.iszero(f::Fun) = iszero(f.values)

# LinearAlgebra.norm(f::Fun) = norm(f.values)
# LinearAlgebra.norm(f::Fun, p::Real) = norm(f.values, p)

export unit
function Forms.unit(::Type{Fun{D,P,R,C,S,T}}, manifold::Manifold{D,C,S},
                    n::Int) where {D,P,R,C,S,T}
    nelts = nsimplices(manifold, R)
    @assert 1 ≤ n ≤ nelts
    # return Fun{D,P,R,C,S}(manifold, sparsevec([n], [one(T)],nelts))
    return Fun{D,P,R,C,S}(manifold, IDVector{R}(T[i == n for i in 1:nelts]))
end

function Base.:+(f::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, +f.values)
end

function Base.:-(f::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, -f.values)
end

function Base.:+(f::Fun{D,P,R,C,S}, g::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    @assert f.manifold == g.manifold
    return Fun{D,P,R,C,S}(f.manifold, f.values + g.values)
end

function Base.:-(f::Fun{D,P,R,C,S}, g::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    @assert f.manifold == g.manifold
    return Fun{D,P,R,C,S}(f.manifold, f.values - g.values)
end

function Base.:*(a::Number, f::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, a * f.values)
end

function Base.:\(a::Number, f::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, a \ f.values)
end

function Base.:*(f::Fun{D,P,R,C,S}, a::Number) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, f.values * a)
end

function Base.:/(f::Fun{D,P,R,C,S}, a::Number) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, f.values / a)
end

Base.:~(f::Fun) = map(~, f)
Base.:&(f::Fun, g::Fun) = map(&, f, g)
Base.:|(f::Fun, g::Fun) = map(|, f, g)

# Functions are a category

export id
function id(::Type{<:Fun{D,P,0,C,S,SVector{C,S}}},
            manifold::Manifold{D,C,S}) where {D,P,C,S}
    return Fun{D,P,0,C,S}(manifold, get_coords(manifold))
end

function Base.conj(f::Fun{D,P,R,C,S}) where {D,P,R,C,S}
    return Fun{D,P,R,C,S}(f.manifold, conj(f.values))
end

# composition?

end


text:
abstract type EmpiricalBayesEstimator end
abstract type LinearEstimator <: EmpiricalBayesEstimator end


text:
module GPIOJetsu

using Mmap

################################################################################
#                                  Constants
################################################################################

include("./constants.jl")

################################################################################
#                                    Types
################################################################################

include("./types.jl")

################################################################################
#                                   Includes
################################################################################

include("./cmacros.jl")
include("./ioctl.jl")
include("./gpio.jl")
include("./i2c.jl")
include("./spi.jl")

################################################################################
#                                Initialization
################################################################################

end # module


text:
using ProjectEuler
using Base.Test

@test ProjectEuler.problem001() == 233168
@test ProjectEuler.problem002() == 4613732
@test ProjectEuler.problem003() == 6857
@test ProjectEuler.problem004() == 906609
@test ProjectEuler.problem005() == 232792560
@test ProjectEuler.problem006() == 25164150
@test ProjectEuler.problem007() == 104743
@test ProjectEuler.problem008() == 23514624000
@test ProjectEuler.problem009() == 31875000