text:
import BITEModel
const BM=BITEModel
using Plots, StatPlots, Measures, NamedTuples

#include("itmix-comp-fit-sigma.jl")
# include("../rgi-setup.jl")
# region = [3, 4, 14, 11, 7][end]
# runtyp = "test"
# update_cache = true
# glnrs1, glnrs2, glnrs_all = get_glaciers(region)
# sigma_of_model = @NT(sigma_h_model = 10.0,
#                      sigma_iv_model = 10.0)
# fit_sigma = false
# nr = glnrs2[end]

# # number of divisions in elevation band variables
# fit_vars = @NT(n_1d_btilde=3,
#                n_1d_fsl=3,
#                n_1d_temp=1)
# const progress_bar = true
# fit_target = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.length][1]
# out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(1), dir_output="", retsol=true, fit_sigma=fit_sigma)

fl = [# not working anymore because of type updates:
    # "figs/f2-5-34-42-h.jls",
      # "figs/f2-5-34-42-h_iv-abserr.jls",
      # "figs/f2-5-34-42-h-abserr.jls",
      # "figs/f2-5-34-42-h_iv-abserr-sigma20.jls",
      # "figs/f2-5-34-42-h_iv-abserr-thin3.jls",
      # "figs/f2-5-34-42-h-abserr-sigmah20-thin3.jls", # looks sweet
      # "figs/f2-5-34-42-h-abserr-thin3.jls",
      # "figs/f2-5-34-42-h-abserr-sigma30-thin1-lengths1e-3.jls",
      # "figs/f2-5-34-42-h-sqrerr-sigma20-thin1-lengths1e-3.jls",
      # "figs/f2-mar20-5-34-42-h-abserr-sigma30-thin1.jls",
      # "figs/f2-mar20-5-34-42-h-sqrerr-sigma40-thin1.jls",
      # "figs/f2-mar21-5-34-42-h-sqrerr-sigma50-thin1.jls",
      # "figs/f2-mar21-5-34-42-h-sqrerr-sigma60-thin1.jls",
      # "figs/f2-mar21-5-34-42-h-sqrerr-sigma70-thin1.jls",
      # "figs/f2-mar21-5-34-42-h-abserr-sigma20-thin1.jls",
      # "figs/f2-mar22-5-34-42-h-abserr-sigma20-thin1.jls",
      "figs/f2-apr01-5-34-42-h-sqrerr-sigma20-thin1.jls",
      # fitting sigma too
      "figs/f2-with_sigma-apr03-5-34-42-h-sqrerr.jls", # with Gaussian sigma prior
      "figs/f2-with_sigma-apr04-5-34-42-h_iv-sqrerr.jls", # Uniform
      "figs/f2-with_sigma-apr04-5-34-42-h-sqrerr.jls", # Uniform: very similar to Gaussian
      "figs/f2-with_sigma-apr04-5-34-42-iv-sqrerr.jls", # Uniform: doesn't actually look that bad...
      # noveau
      "figs/f2-with_sigma-better-apr16-5-34-42-h_iv-sqrerr.jls", # with sigma_btilde = 0.3
      "figs/f2-large-btilde-sigma-apr17-5-34-42-h_iv-sqrerr.jls", # with sigma_btilde = 2.8
      "figs/f2-sigma-btilde-1.4-sigma-apr17-5-34-42-h_iv-sqrerr.jls",  # with sigma_btilde = 1.4 (production value)
      #"figs/f2-sigma-btilde-1.4_sigma-dhdt-26-sigma-apr17-5-34-42-h_iv-sqrerr.jls",  # WRONG with sigma_btilde = 1.4 (production value), =26 (production value)
      "figs/f2-production-sigma-apr17-5-34-42-h_iv-sqrerr.jls",
      "figs/f2-production-length97-sigma-apr17-5-34-42-h_iv-sqrerr.jls", # <--- in publication
      "figs/f2-production-length98-sigma-apr17-5-34-42-h_iv-sqrerr.jls", # BAD, something is off here, and for sure not converged.
      "figs/f2-production-length98-btildes-0.3-apr17-5-34-42-h_iv-sqrerr.jls", # back to small btilde error.  Sigma h-error is large
      "figs/f2-production-length97-btildes-0.3-ball-apr17-5-34-42-h_iv-sqrerr.jls", #  <--- in publication, h-error is small
      "figs/f2-production-length98-btildes-0.3-ball-apr17-5-34-42-h_iv-sqrerr.jls", # Sigma h-error is large
      "figs/f2-production-length995-btildes-0.3-ball-apr17-5-34-42-h_iv-sqrerr.jls", # trying super small length sigmab: looks very good apart from super large sigmas, because the lowest profile is not well fitted
      ][end-2]
@show fl

if true # !isdefined(:sol)
    Revise.track(Base) # base/serialize.jl L971 need to modify it
    sol, sol_prior = open(deserialize, fl, "r")

    gl = sol.theta0.gb.gl
end

if !isdefined(:gl_all_radar)
    glacier = :Unteraar
    gid = BM.ITMIXGlacier(glacier, 2)
    gl_all_radar,gb,pp,pm,pn,pl = BM.init_forward(gid)
end

err, bed, surf, rad = BM.error_h(sol.mh2d, gl_all_radar);
err_max, bed_max = BM.error_h(sol.max_h2d, gl_all_radar);
err_min, bed_min = BM.error_h(sol.min_h2d, gl_all_radar);
_, _,_, _, bed_std = BM.error_h(sol.sh2d, gl_all_radar);

for (ii,sps) in enumerate([[7, 16, 25, 38, 45], [5, 34, 42, 16, 25]]) # fit to 5, 34, 42
    # map view
    p = plot(aspect_ratio=:equal);
    for (ss,s) in enumerate(err.splits)
        plot!(err.x[s], err.y[s], label="", lw=:0.5, lc=:black)
        #annotate!(err.x[s][1], err.y[s][1], Text("$ss"))
    end

    # fitting tracks
    for (ss,s) in enumerate(gl.h.vals.splits)
        plot!(gl.h.vals.x[s], gl.h.vals.y[s], lw=3, label="", lc=:black)
    end

    str = ["A" "B" "C" "D" "E"]
    for (i,ss) in enumerate(sps)
        s = err.splits[ss]
        plot!(gl_all_radar.h.vals.x[s], gl_all_radar.h.vals.y[s], lw=2, label="", lc=:red)
        annotate!(err.x[s][1], err.y[s][1], text("$(str[i])", 15), lc=:green)
    end
    BM.plot_outlines!(gl_all_radar, xyscale=1, lw=0.5, xlabel="", ylabel="", grid=false, framestyle=:none)
    plot!([440e3], [5155e3], xerror=[-1e3,1e3], markercolor=:black, markerstrokecolor=:black, markerstrokewidth=2,
          markersize=10, label="", left_margin=10mm; ((:xlims, (432.2e3,442.7e3)), (:ylims, (5153.8e3, 5162.2e3)))...)
    annotate!(440e3, 5155.4e3, text("2 km", 10))

    #display(p)

    ps = []
    for (i,ss) in enumerate(sps)
        onleft = i==1 || i==4
        onbottom = i>3
        s = err.splits[ss]
        dist = cumsum([0, sqrt.(diff(err.x[s]).^2 + diff(err.y[s]).^2)...])
        push!(ps,
              plot(dist, surf.v[s], lw=2, color=:green,  label= i==length(str) ? "surface" : ""))
        plot!(dist, rad.v[s], lw=2, label= i==length(str) ? "radar bed" : "", markershape=:circle, markersize=1.5,
              color=:red, markercolor=:black)
        plot!(dist, rad.v[s]+27, linecolor=:red, linewidth=0,
              fillrange=rad.v[s]-27, fillcolor=:red, fillalpha=0.3,
              label= i==length(str) ? "   ± 27m" : "")
        plot!(dist, bed.v[s], lw=1.5, linecolor=:blue, label= i==length(str) ? "mean bed" :"",
              xlabel= onbottom ? "distance (m)" : "", ylabel= onleft ? "elevation (m a.s.l.)" : "")
        plot!(dist, bed.v[s]+2bed_std.v[s], linecolor=:blue, linewidth=0,
              fillrange=bed.v[s]-2bed_std.v[s], fillcolor=:blue, fillalpha=0.3,
              label= i==length(str) ? "   2x std" :"")
        plot!(dist, bed_max.v[s], lw=0.5, color=:black, label= i==length(str) ? "    extrema" : "")
        plot!(dist, bed_min.v[s], lw=0.5, color=:black, label="")
        # if i==length(str)
        #     # make labels
        #     x, y = [err.x[1:1],dist[1:1]][2] , [NaN]
        #     # plot!(x,y, lw=2, label="surface")
        #     # plot!(x,y, lw=2, label="radar bed ± 10m")
        #     plot!(x,y, lw=1.5, label= "mean")
        #     plot!(x,y, fillrange=y, linecolor=:blue, linewidth=0,
        #           fillcolor=:blue, fillalpha=0.3,
        #           label="2x std")
        #     plot!(x,y, lw=0.5, label="max")
        #     plot!(x,y, lw=0.5, label="min")
        # end

    end

    l = @layout [grid(1,3)
                 grid(1,3)] # a{0.66w}]
    kwargsj = [(:tickfontsize, 10), (:legendfontsize, 10),
               (:xguidefontsize, 12), (:yguidefontsize, 12), (:titlefontsize, 20)]

    pend = plot(ps..., p, layout=l, title=["A" "B" "C" "D" "E" "F"], title_location=:right, size=(900,600), reuse=false; kwargsj...)
    display(pend)
    if ii==1
        savefig("figs/f4.png")
        savefig("figs/f4.pdf")
    else
        savefig("figs/f4-fitting-profiles.png")
        savefig("figs/f4-fitting-profiles.pdf")
    end

end

#BM.plottheta(sol.thetas, sol.theta0, reuse=false)


text:
"""
    boundaries(W::TriangulationWomble{T}; threshold::T=0.1) where {T <: Number}

Extracts candidate boundaries using calculated rates of change (M) on specified
threshold. Default threshold is 10%, meaning that the top 10% of pixels are
selected as part of the boundaries. This function returns a list of indices
identifying which simplices are part of the boundaries. The NaN values in the
rates of change are not going to be a part of the boundaries.
"""
function boundaries(W::TriangulationWomble{T}; threshold::T=0.1) where {T<:Number}
    limit = floor(Int, length(W.m) * threshold)
    nans = count(isnan, W.m)
    return findall(nans .< denserank(W; rev=true) .<= (limit+nans))
end

"""
    boundaries(W::LatticeWomble{T}; threshold::T=0.1) where {T<:Number}

Extracts candidate boundaries using calculated rates of change (M) on specified
threshold. Default threshold is 10%, meaning that the top 10% of pixels are
selected as part of the boundaries. This function returns an array of Cartesian
indices corresponding to the rate of change matrix, indicating which positions
are part of the boundaries. The NaN values in the rates of change are not going
to be a part of the boundaries.
"""
function boundaries(W::LatticeWomble{T}; threshold::T=0.1) where {T<:Number}
    limit = floor(Int, size(W.m, 2) * size(W.m, 1) * threshold)
    nans = count(isnan, W.m)
    return findall(nans .< denserank(W; rev=true) .<= (limit+nans))
end

ext: jl


text:
# TODO: more tests!

using ChainRulesCore, Test
using LinearAlgebra: Diagonal
using ChainRulesCore: extern, accumulate, accumulate!, store!, @scalar_rule,
    Wirtinger, wirtinger_primal, wirtinger_conjugate, add_wirtinger, mul_wirtinger,
    Zero, add_zero, mul_zero, One, add_one, mul_one, Casted, cast, add_casted, mul_casted,
    DNE, Thunk, Casted, DNERule
using Base.Broadcast: broadcastable

@testset "ChainRulesCore" begin
    include("differentials.jl")
    include("rules.jl")
end


text:
#Test if tf listner works correctly
using RobotOS

# # Initialize node
# init_node("listner", anonymous=true) # delete

path = @__DIR__

# Need to run the listener and broadcaster in parralllel
script = (`julia --project=@. $path/tf2_ros_broadcaster.jl` &
          `julia --project=@. $path/tf2_ros_listener.jl`);
Base.run(script);


text:
struct Boat_sim_pars
    bh2o::Float64
    boat_angles::Boat_angles
    η::Float64
    x1va::Float64
    x1air::Float64
    x′::Float64
    z′::Float64
end

function Q(θ, bsp::Boat_sim_pars)
    x1va, x1air, x′, ba, bh2o = bsp.x1va, bsp.x1air, bsp.x′, bsp.boat_angles, bsp.bh2o
    Q_g(θ, ba.γ_OA, ba.θ_k, ba.θ_t) +
    Q_BUOY(θ, bh2o) +
    Q_VA(θ, x′, bh2o, x1va) +
    Q_AE(θ, x′, bh2o, x1air) +
    Q_ROAE(θ, x′, x1air, ba.θ_t, ba.θ′_t) +
    Q_T(θ, x′, x1va, ba.γ_OA, ba.γ′_OA, bsp.η)
end

function Q_g(θ, γ_OA, θ_k, θ_t)
    sθ, cθ = sincos(θ)
    Xg = m * g * sθ
    Zg = m * g * cθ
    Mg_RG = m_RG * g * r_RG[1] * cθ
    Mg_w = m_w * g * x_w(θ_t) * cθ
    Mg_fs = m_fs * g * x_SB * cθ
    Mg_f = m_f * g * x_SB * cθ
    Mg_k = m_k * g * dot(r_w(θ_t) + r_k_const(θ_t) + 0.5l_k * r_k_∂l_k(θ_k), SA[cθ, -sθ])
    Mg_t =
        m_t * g * ((-0.5l_t * cos(θ_t) + 2x_w(θ_t)) * cos(θ) - 0.5l_t * sin(θ_t) * sθ)
    Mg_s =
        m_s * g * ((0.5l_s * cos(F_0(θ_t)) + 2x_SB) * cos(θ) - 0.5l_t * sin(θ_t) * sθ)
    # Here is also a difference
    Mg_OA = m_OA * g * 0.5(l_OAMX + l_OAMN) * x_OA_∂l_OA(γ_OA) * cθ
    # moa * g * (loamx + loamn) * cos(gamoa) * (cos(θ) + sin(θ))
    @debug "Mg parts are as follows: " Mg_RG Mg_w Mg_fs Mg_f Mg_k Mg_t Mg_s Mg_OA
    Mg = Mg_RG + Mg_w + Mg_fs + 2Mg_f + Mg_k + 2Mg_t + 2Mg_s + 2Mg_OA
    SA[Xg, Zg, Mg]
end

function Q_BUOY(θ, bh2o)
    F_BUOY = -rhoh2o * g * vh2o(θ, bh2o)
    X = F_BUOY * sin(θ)
    Z = F_BUOY * cos(θ)
    M_BUOY = -rhoh2o * g * (szvh2o(θ, bh2o) * sin(θ) + sxvh2o(θ, bh2o) * cos(θ))
    SA[X, Z, M_BUOY]
end

function Q_VA(θ, x1b, bh2o, x1va)
    x1 = x1b - x1va
    sx1 = sign(x1)
    coeff = 0.5rhoh2o * x1^2 * sx1

    R_SHAPE = coeff * cdx * tvmh2o(θ) / cos(θ)
    R_VIS = coeff * cf0 * svbh2o(θ, bh2o) * freh2o
    R_VAVE = coeff * cdw * yth2ov(θ, bh2o)
    R_VA = R_SHAPE + R_VIS + R_VAVE
    X_VA = -R_VA * cos(θ)
    Z_VA = R_VA * sin(θ)
    M_SHAPE = -coeff * cdx * sztvmh2o(θ, bh2o) ./ cos(θ)
    M_VIS = -coeff * cf0 * (szsvbh2o(θ, bh2o) * cos(θ) - sxsvbh2o(θ, bh2o) * sin(θ)) * freh2o
    M_VAVE = -coeff * cdw * (z1yth2ov(θ, bh2o) * cos(θ) - x1yth2ov(θ, bh2o) * sin(θ))
    M_VA = M_SHAPE + M_VIS + M_VAVE
    SA[X_VA, Z_VA, M_VA]
end

function Q_AE(θ, x1b, bh2o, x1air)
    x1 = x1b - x1air

    coeff = 0.5rhoair * x1^2 * sign(x1)
    R_AE = coeff * sbair(θ, bh2o) * cf0 * freair
    X_AE = -R_AE * cos(θ)
    Z_AE = R_AE * sin(θ)
    M_AE = -coeff * (szsbair(θ, bh2o) * cos(θ) - sxsbair(θ, bh2o) * sin(θ)) * cf0 * freair
    SA[X_AE, Z_AE, M_AE]
end

function Q_T(θ, x1b, x1va, γ_OA, γ′_OA, η)
    F_OAR = foa(γ_OA, γ′_OA, x1b, x1va, η)
    T_OAR = 2F_OAR * sin(γ_OA)
    X_T = T_OAR * cos(θ)
    Z_T = -T_OAR * sin(θ)
    SA[X_T, Z_T, 0]
end

function Q_ROAE(θ, x1b, x1air, θ_t, θ′_t)
    x1w = 0
    # or
    # x1w = x′_w(θ_t, θ′_t)
    x1 = x1b - x1air + x1w * cos(θ)

    coeff = 0.5rhoair * x1^2 * sign(x1)
    R_ROAE = coeff * sro * cdro
    X_ROAE = -R_ROAE * cos(θ)
    Z_ROAE = R_ROAE * sin(θ)
    M_ROAE = -coeff * szsro * cdro
    SA[X_ROAE, Z_ROAE, M_ROAE]
end


text:
module nfxptest

    using Base.Test,nfxp

    @testset "simulated data" begin
        N = 50
        T = 150
        mod = nfxp.Param()
        sd = nfxp.simdata(N,T,mod)

        @test length(sd[:id]) == N*T

        @test extrema(sd[:id]) == (1,N)
        @test extrema(sd[:t]) == (1,T)
        @test extrema(sd[:d]) == (0,1)

        # whenever replace, next period mileage must be the same as dx1
        # @test all(sd[:x1][find(sd[:d])] .== sd[:dx1][find(sd[:d])])

        # mileage can never exceed p.n
        @test all(sd[:x] .<= mod.n)

        # P is a proper transition matrix
        @test all(sum(mod.P,2).==1.0)
    end

    @testset "check likelihood gradient" begin
        p = nfxp.Param()
        d = nfxp.simdata(10000,100,p)
        g = zeros(p.n_params)
        f = nfxp.likelihood!(p.theta,g,d,p)
        println(g)
    end

end

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
using PLite
using Base.Test

# todo


text:
@testset "Low-level" begin
    @testset "Definitions" begin
        test(StructDefinition{false}, struct_by_name, :VkPhysicalDeviceProperties, :(
            struct PhysicalDeviceProperties <: ReturnedOnly
                api_version::VersionNumber
                driver_version::VersionNumber
                vendor_id::UInt32
                device_id::UInt32
                device_type::PhysicalDeviceType
                device_name::String
                pipeline_cache_uuid::String
                limits::PhysicalDeviceLimits
                sparse_properties::PhysicalDeviceSparseProperties
            end))

        test(StructDefinition{false}, struct_by_name, :VkApplicationInfo, :(
            struct _ApplicationInfo <: VulkanStruct{true}
                vks::VkApplicationInfo
                deps::Vector{Any}
            end))

        test(StructDefinition{false}, struct_by_name, :VkExtent2D, :(
            struct _Extent2D <: VulkanStruct{false}
                vks::VkExtent2D
            end))

        test(StructDefinition{false}, struct_by_name, :VkExternalBufferProperties, :(
            struct ExternalBufferProperties <: ReturnedOnly
                s_type::StructureType
                next::Ptr{Cvoid}
                external_memory_properties::ExternalMemoryProperties
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkPipelineExecutableInternalRepresentationKHR, :(
            struct PipelineExecutableInternalRepresentationKHR <: ReturnedOnly
                s_type::StructureType
                next::Ptr{Cvoid}
                name::String
                description::String
                is_text::Bool
                data_size::UInt
                data::Ptr{Cvoid}
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkDescriptorSetAllocateInfo, :(
            struct _DescriptorSetAllocateInfo <: VulkanStruct{true}
                vks::VkDescriptorSetAllocateInfo
                deps::Vector{Any}
                descriptor_pool::DescriptorPool
            end
        ))

        test(StructDefinition{false}, struct_by_name, :VkAccelerationStructureBuildGeometryInfoKHR, :(
            struct _AccelerationStructureBuildGeometryInfoKHR <: VulkanStruct{true}
                vks::VkAccelerationStructureBuildGeometryInfoKHR
                deps::Vector{Any}
                src_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
                dst_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
            end
        ))
    end

    @testset "Friendly constructors" begin
        test(Constructor, StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo)), :(
            function _InstanceCreateInfo(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; next=C_NULL, flags=0, application_info=C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                application_info = cconvert(Ptr{VkApplicationInfo}, application_info)
                enabled_layer_names = cconvert(Ptr{Cstring}, enabled_layer_names)
                enabled_extension_names = cconvert(Ptr{Cstring}, enabled_extension_names)
                deps = [
                    next,
                    application_info,
                    enabled_layer_names,
                    enabled_extension_names
                ]
                vks = VkInstanceCreateInfo(VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO,
                                        unsafe_convert(Ptr{Cvoid}, next),
                                        flags,
                                        unsafe_convert(Ptr{VkApplicationInfo}, application_info),
                                        pointer_length(enabled_layer_names),
                                        unsafe_convert(Ptr{Cstring}, enabled_layer_names),
                                        pointer_length(enabled_extension_names),
                                        unsafe_convert(Ptr{Cstring}, enabled_extension_names),
                                        )
                _InstanceCreateInfo(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkSubpassSampleLocationsEXT)), :(
            function _SubpassSampleLocationsEXT(subpass_index::Integer, sample_locations_info::_SampleLocationsInfoEXT)
                _SubpassSampleLocationsEXT(VkSubpassSampleLocationsEXT(subpass_index, sample_locations_info.vks))
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDebugUtilsMessengerCreateInfoEXT)), :(
            function _DebugUtilsMessengerCreateInfoEXT(message_severity::DebugUtilsMessageSeverityFlagEXT, message_type::DebugUtilsMessageTypeFlagEXT, pfn_user_callback::FunctionPtr; next = C_NULL, flags = 0, user_data = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                user_data = cconvert(Ptr{Cvoid}, user_data)
                deps = [next, user_data]
                vks = VkDebugUtilsMessengerCreateInfoEXT(VK_STRUCTURE_TYPE_DEBUG_UTILS_MESSENGER_CREATE_INFO_EXT, unsafe_convert(Ptr{Cvoid}, next), flags, message_severity, message_type, pfn_user_callback, unsafe_convert(Ptr{Cvoid}, user_data))
                _DebugUtilsMessengerCreateInfoEXT(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorSetAllocateInfo)), :(
            function _DescriptorSetAllocateInfo(descriptor_pool, set_layouts::AbstractArray; next = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                set_layouts = cconvert(Ptr{VkDescriptorSetLayout}, set_layouts)
                deps = [next, set_layouts]
                vks = VkDescriptorSetAllocateInfo(VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO, unsafe_convert(Ptr{Cvoid}, next), descriptor_pool, pointer_length(set_layouts), unsafe_convert(Ptr{VkDescriptorSetLayout}, set_layouts))
                _DescriptorSetAllocateInfo(vks, deps, descriptor_pool)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkPipelineShaderStageCreateInfo)), :(
            function _PipelineShaderStageCreateInfo(stage::ShaderStageFlag, _module, name::AbstractString; next = C_NULL, flags = 0, specialization_info = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                name = cconvert(Cstring, name)
                specialization_info = cconvert(Ptr{VkSpecializationInfo}, specialization_info)
                deps = [next, name, specialization_info]
                vks = VkPipelineShaderStageCreateInfo(VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO, unsafe_convert(Ptr{Cvoid}, next), flags, VkShaderStageFlagBits(stage.val), _module, unsafe_convert(Cstring, name), unsafe_convert(Ptr{VkSpecializationInfo}, specialization_info))
                _PipelineShaderStageCreateInfo(vks, deps, _module)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorImageInfo)), :(
            function _DescriptorImageInfo(sampler, image_view, image_layout::ImageLayout)
                _DescriptorImageInfo(VkDescriptorImageInfo(sampler, image_view, image_layout), sampler, image_view)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkDescriptorSetLayoutBinding)), :(
            function _DescriptorSetLayoutBinding(binding::Integer, descriptor_type::DescriptorType, stage_flags::ShaderStageFlag; descriptor_count = 0, immutable_samplers = C_NULL)
                immutable_samplers = cconvert(Ptr{VkSampler}, immutable_samplers)
                deps = [immutable_samplers]
                vks = VkDescriptorSetLayoutBinding(binding, descriptor_type, descriptor_count, stage_flags, unsafe_convert(Ptr{VkSampler}, immutable_samplers))
                _DescriptorSetLayoutBinding(vks, deps)
            end
        ))

        test(Constructor, StructDefinition{false}(struct_by_name(:VkXcbSurfaceCreateInfoKHR)), :(
            function _XcbSurfaceCreateInfoKHR(connection::Ptr{vk.xcb_connection_t}, window::vk.xcb_window_t; next = C_NULL, flags = 0)
                next = cconvert(Ptr{Cvoid}, next)
                connection = cconvert(Ptr{vk.xcb_connection_t}, connection)
                deps = [next, connection]
                vks = VkXcbSurfaceCreateInfoKHR(VK_STRUCTURE_TYPE_XCB_SURFACE_CREATE_INFO_KHR, unsafe_convert(Ptr{Cvoid}, next), flags, unsafe_convert(Ptr{vk.xcb_connection_t}, connection), window)
                _XcbSurfaceCreateInfoKHR(vks, deps)
            end
        ))
    end

    @testset "Manual tweaks" begin
        test(Constructor, StructDefinition{false}(struct_by_name(:VkWriteDescriptorSet)), :(
            function _WriteDescriptorSet(dst_set, dst_binding::Integer, dst_array_element::Integer, descriptor_type::DescriptorType, image_info::AbstractArray, buffer_info::AbstractArray, texel_buffer_view::AbstractArray; next = C_NULL, descriptor_count = max(pointer_length(image_info), pointer_length(buffer_info), pointer_length(texel_buffer_view)))
                next = cconvert(Ptr{Cvoid}, next)
                image_info = cconvert(Ptr{VkDescriptorImageInfo}, image_info)
                buffer_info = cconvert(Ptr{VkDescriptorBufferInfo}, buffer_info)
                texel_buffer_view = cconvert(Ptr{VkBufferView}, texel_buffer_view)
                deps = [next, image_info, buffer_info, texel_buffer_view]
                vks = VkWriteDescriptorSet(VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET, unsafe_convert(Ptr{Cvoid}, next), dst_set, dst_binding, dst_array_element, descriptor_count, descriptor_type, unsafe_convert(Ptr{VkDescriptorImageInfo}, image_info), unsafe_convert(Ptr{VkDescriptorBufferInfo}, buffer_info), unsafe_convert(Ptr{VkBufferView}, texel_buffer_view))
                _WriteDescriptorSet(vks, deps, dst_set)
            end
        ))
    end
end


text:
# XXX Should this permute the dimensions?
dims(H::Hermitian{<:Number,<:Tensor}) = dims(parent(H))

# XXX Should this permute the dimensions?
dim(H::Hermitian{<:Number,<:Tensor}, i::Int) = dim(parent(H), i)

matrix(H::Hermitian{<:Number,<:Tensor}) = Hermitian(matrix(parent(H)))

# XXX Should this permute the indices?
inds(H::Hermitian{<:Number,<:Tensor}) = inds(parent(H))

# XXX Should this permute the indices?
ind(H::Hermitian{<:Number,<:Tensor}, i::Int) = ind(parent(H), i)

# XXX Should this tranpose the block locations?
nnzblocks(H::Hermitian{<:Number,<:Tensor}) = nnzblocks(parent(H))

# XXX Should this tranpose the block locations?
nzblocks(H::Hermitian{<:Number,<:Tensor}) = nzblocks(parent(H))

# XXX Should this tranpose the block locations?
eachnzblock(H::Hermitian{<:Number,<:Tensor}) = eachnzblock(parent(H))

nblocks(H::Hermitian{<:Number,<:Tensor}) = nblocks(parent(H))

blockview(H::Hermitian{<:Number,<:Tensor}, block) = Hermitian(blockview(parent(H), block))


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ aff6e710-27c0-11eb-1782-91c83dd55c13
using Automata

# ╔═╡ c2dd9e00-27c0-11eb-148d-d5fc9de0997f
using Images

# ╔═╡ 8bd97320-27c0-11eb-0a34-535cdfded71b
md"
## Random initial conditions for cellular automata
"

# ╔═╡ d90af470-27c0-11eb-3bc5-abb917b5234d
Gray.(1 .- Float32.(Automata.to_array(elementary(30, 1000))))

# ╔═╡ 774e3e80-27c1-11eb-19f2-25fbb633dc17
mkdir("random_initial")

# ╔═╡ 68a6744e-27c2-11eb-2082-796e0b0a16bb
mkdir("elementary")

# ╔═╡ a59a0d00-27c1-11eb-1db6-8fa51f230fe1
rule = 110

# ╔═╡ ad300970-27c1-11eb-0db8-57b9b7e88b71
steps = 1000

# ╔═╡ 7b389810-27c1-11eb-1c82-492fa5aad6cb
for i = 1:1000
	initial_state = zeros(Bool, Int(ceil(log2(i+1))))
	for j = 1:length(initial_state)
		initial_state[j] = i >> (j - 1) & 1
	end
	render("random_initial/elementary/$i.png", elementary(rule, steps, initial_state))
end

# ╔═╡ 46d0b100-27c3-11eb-2755-75a521586609
mkdir("totalistic")

# ╔═╡ b90af0a0-27c3-11eb-2ea6-1723d1eaf6cd
mkdir("totalistic_diff")

# ╔═╡ 85bee130-27c2-11eb-0209-4be164f25ffb
colors = 4

# ╔═╡ 9a4e8d30-27c2-11eb-06ac-2101484da7d5
code = 2530

# ╔═╡ a3224450-27c3-11eb-3cc3-6b7060b71981
diff_arr(arr) = abs.(diff(arr, dims=2))

# ╔═╡ 780bf3c0-27c2-11eb-062b-c500ef9e2dee
for i = 1:10
	initial_state = zeros(UInt8, Int(ceil(log(colors, i+1))))
	acc = i
	for j = length(initial_state):-1:1
		digit = colors ^ (j - 1)
		initial_state[length(initial_state) - j + 1] = acc ÷ digit
		acc %= digit
	end
	render("random_initial/totalistic_diff/$i.png", totalistic(code, steps, colors, initial_state), [diff_arr])
end

# ╔═╡ Cell order:
# ╟─8bd97320-27c0-11eb-0a34-535cdfded71b
# ╠═aff6e710-27c0-11eb-1782-91c83dd55c13
# ╠═c2dd9e00-27c0-11eb-148d-d5fc9de0997f
# ╠═d90af470-27c0-11eb-3bc5-abb917b5234d
# ╠═774e3e80-27c1-11eb-19f2-25fbb633dc17
# ╠═68a6744e-27c2-11eb-2082-796e0b0a16bb
# ╠═a59a0d00-27c1-11eb-1db6-8fa51f230fe1
# ╠═ad300970-27c1-11eb-0db8-57b9b7e88b71
# ╠═7b389810-27c1-11eb-1c82-492fa5aad6cb
# ╠═46d0b100-27c3-11eb-2755-75a521586609
# ╠═b90af0a0-27c3-11eb-2ea6-1723d1eaf6cd
# ╠═85bee130-27c2-11eb-0209-4be164f25ffb
# ╠═9a4e8d30-27c2-11eb-06ac-2101484da7d5
# ╠═a3224450-27c3-11eb-3cc3-6b7060b71981
# ╠═780bf3c0-27c2-11eb-062b-c500ef9e2dee


text:
using FourierTools

f(x) = sin(x * 2π)
x = [0.0, 0.2, 0.4, 0.6, 0.8]

@show(f.(x))
@show(f.(x .- 0.5))
@show(shift(f.(x), -2.5))

plot(x, f.(x))
#plot!(range(0, 1, length=100), resample(f.(x), (100,)))
plot!(x .- 0.5, f.(x .- 0.5))
plot!(x .- 0.5, shift(f.(x), -2.5))


text:
Base.CoreLogging.disable_logging(Base.CoreLogging.Error)

import PredictMDExtra

PredictMDExtra._print_list_of_package_imports()


text: include("TestModuleBroken.jl")

max_stars_repo_name: tgflynn/IncludeGuards

max_issues_repo_name: tgflynn/IncludeGuards

max_forks_repo_name: tgflynn/IncludeGuards


text:
# functions for neid solar spectra

#function to read a set of NEID solar spectra and filter based on hour angle
#deprecated because data set is too large to load into memory at once - using daily averages instead (see: combine_NEID_daily_obs)
function read_and_filter_all_spectra(data_path::String; verbose::Bool=false, max_solar_hour_angle::Float=0.0)
   if verbose println("# Finding what data files are avaliable.")  end
   if isfile(joinpath(pipeline_output_path_afw5465,"manifest.csv"))
      if verbose println("# Reading in manifest from manifest.csv") end
      df_files  = CSV.read(joinpath(pipeline_output_path_afw5465,"manifest.csv"), DataFrame)
      @assert size(df_files,1) >= 1
      @assert hasproperty(df_files,:Filename)
      @assert hasproperty(df_files,:bjd)
   else
      if verbose println("# Generating manifest file manifest.csv") end
      if inst==NEID
         df_files = make_manifest_NEID(data_path, max_spectra_to_use=max_spectra_to_use)
      else
         df_files = inst.make_manifest(data_path, inst)
      end
      CSV.write(joinpath(pipeline_output_path_afw5465,"manifest.csv"), df_files)
   end
   if @isdefined fits_target_str
      df_files_use = df_files |>
      @filter( _.target == parse(Int64,fits_target_str) ) |>
      @filter( max_solar_hour_angle == 0.0 || abs(_.hour_angle) <= max_solar_hour_angle ) |>
      @take(max_spectra_to_use) |> DataFrame
   else
      df_files_use = df_files |>
      @filter( max_solar_hour_angle == 0.0 || abs(_.hour_angle) <= max_solar_hour_angle ) |>
      @take(max_spectra_to_use) |> DataFrame
   end
   if verbose println("# Reading in ", size(df_files_use,1), " FITS files.")  end
   @time all_spectra = map(NEID.read_data,eachrow(df_files_use))
end

##########################################################################################
#scripts that were aimed towards making a daily template spectrum - keeping this here for example code but should not assume it is working
##########################################################################################

pipeline_plan = PipelinePlan()
verbose = true
reset_all_needs!(pipeline_plan)

if verbose println("# Reading in customized parameters from param.jl.")  end
paths_to_search_for_param = [pwd(),joinpath(pwd(),"inputs"),joinpath(pkgdir(RvSpectMLBase),"..","RvLineList","inputs")]
eval(read_data_paths(paths_to_search=paths_to_search_for_param, filename="param.jl"))

"""
global norm_type
if norm_type != :continuum
   println("Warning: norm_type=:continuum is required by generateEmpiricalMask() but found norm_type=:", String(norm_type), ". Setting global norm_type=:continuum in order to continue.")
   norm_type = :continuum
end
"""

if need_to(pipeline_plan,:read_spectra)
   all_spectra = read_spectra(joinpath(data_path,date_subdir), max_solar_hour_angle=1.0)
   dont_need_to!(pipeline_plan,:read_spectra)
end

#make sure spectra are normalized
#all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_neid_solar_data_20190918.jl"))
#if typeof(get_inst(all_spectra)) <: AnyEXPRES   RvLineList.continuum_normalize_spectra!(all_spectra)   end

order_list_timeseries = extract_orders(all_spectra, pipeline_plan, orders_to_use=orders_to_use, recalc=true);

# First, make line list using default orders, where laser comb calibration is available
#linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",espresso_mask_filename)
#line_list_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=true)
#(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_scale_factor=2.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true)
#line_width_50 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.5)
#line_width_05 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)
#line_list_excalibur_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan,  orders_to_use=43:72, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = line_width_05, recalc=true, verbose=true)
#((ccfs, ccf_vars), v_grid) = ccf_total(order_list_timeseries, line_list_excalibur_df, pipeline_plan,  mask_scale_factor=8.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100, calc_ccf_var=true, recalc=true)
#alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
#rvs_ccf = calc_rvs_from_ccf_total(ccfs, ccf_vars, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)

global f_mean
if need_to(pipeline_plan, :template)  # Compute order CCF's & measure RVs
   if verbose println("# Making template spectra.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
   GC.gc()   # run garbage collector for deallocated memory
   map(i->order_list_timeseries.metadata[i][:rv_est] = 0.0, 1:length(order_list_timeseries) )
   # Smothing is broken with GP interpolation.  Need to fix.  In mean time, here's a Sinc interpolation workaround
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries, smooth_factor=2.0) #, alg=:Sinc)
   #@time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries, alg=:Sinc)
   #if save_data(pipeline_plan, :template)
   using JLD2, FileIO
   save(joinpath(pipeline_output_path_afw5465,date_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   #end
   dont_need_to!(pipeline_plan, :template);
end

using RCall
afs_src = joinpath(pwd(),"AFS.R")
R"source($afs_src)"

function calc_continuum_model(spectrum::AbstractSpectra2D; order_idx::Integer )
    possible_pix = get_pixel_range(get_inst(spectrum),order_idx)
    bad_pix = bad_col_ranges(get_inst(spectrum),order_idx)
    pix_rng = EchelleInstruments.calc_complement_index_ranges(possible_pix,bad_pix)
    pix = mapreduce(p->collect(p),vcat,pix_rng)
    afs_inputs = zeros(Float64,length(pix),2)
    afs_inputs[:,1] .= spectrum.λ[pix,order_idx]
    afs_inputs[:,2] .= spectrum.flux[pix,order_idx]
    @assert !any(isnan.(afs_inputs))
    #=
    wv = mapreduce(p->spec.λ[p,order_idx],vcat,pix_rng)
    @assert !any(isnan.(wv))
    inten = mapreduce(p->convert(Vector{Float64},spec.flux[p,order_idx]),vcat,pix_rng)
    @assert !any(isnan.(inten))
    afs_inputs = hcat(wv,inten)
    =#
    #df = DataFrame("wv"=>wv,"intes"=>inten)
    afs_output_R = R"AFS($afs_inputs,0.95,0.25)"
    afs_output = rcopy(afs_output_R)
    continuum = zeros(eltype(spectrum.flux),size(spectrum.flux[:,order_idx]))
    continuum = fill(NaN,size(spectrum.flux[:,order_idx]))
    continuum[pix] .= afs_output
    return continuum
end

using EchelleInstruments.NEID
function calc_continuum_model(spectrum::AbstractSpectra2D )
    vec_of_orders = pmap(ord->calc_continuum_model(spectrum,order_idx=ord), orders_to_use )
    output = fill(NaN, size(spectrum.flux))
    for (i,ord) in enumerate(min_order(get_inst(spectrum)):max_order(get_inst(spectrum)))
        output[:,ord] .= vec_of_orders[i]
    end
    return output
end

spec = all_spectra[1]
spec.flux = f_mean
continuum = calc_continuum_model(spec)
save(joinpath(pipeline_output_path_afw5465,date_subdir * "_continuum.jld2"), Dict("continuum"=>continuum,"orders_to_use"=>orders_to_use, "AFS_params"=>[0.95,0.25]))

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
using DFTools

filename = ARGS[1]
output = ARGS[2]
curdir = pwd()

if filename[1]!="/"
  filename = curdir*"/"*filename
end

if output[1]!="/"
  output = curdir*"/"*output
end

zmat,yticks = read_KPDOS(filename)
heatmap(zmat_ge,legend=false,color=:BuPu,yticks=(yticks[1][5:5:end]+10,["","","",""]),xticks=([0,101,201],["A","Z","U"])
savefig(output)


text:
using EconometricsLH

function raw_data_files_test(dsName)
	@testset "Raw data files" begin
		ds = make_data_settings(dsName);
        for rf in [CollegeStratData.raw_entry_qual_parental(ds),
            CollegeStratData.raw_transfer_regr(ds)]

            @test isfile(CollegeStratData.data_file(rf));
		end

		# Mapping of moments to files
		missList = missing_file_list(ds);
		@test isempty(missList)

		# copy_raw_data_files(trialRun = true);
    end
end

@testset "DataFiles" begin
	for dsName ∈ CollegeStratData.data_settings_list()
		raw_data_files_test(dsName)
	end
end

# -------------


text:
include("./tensile.jl")
include("./flexure.jl")

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
@test_throws MethodError shuffleobs(kfolds(X1))
@test_throws MethodError splitobs(kfolds(X1))
@test_throws MethodError oversample(kfolds(X1))
@test_throws MethodError undersample(kfolds(X1))
@test_throws MethodError stratifiedobs(kfolds(X1))

@testset "kfolds with Int" begin
    @test typeof(kfolds) <: Function
    @test typeof(kfolds(10)) <: Tuple

    @test_throws ArgumentError kfolds(10,1)
    @test_throws ArgumentError kfolds(10,11)
    @test_throws ArgumentError kfolds(1)
    @test_throws ArgumentError kfolds(2)
    @test_throws ArgumentError kfolds(-1)

    for n in (5,10,15,20,50), k in (2,5)
        for (f1,f2) in (@inferred(kfolds(n,k)), @inferred(kfolds(n)))
            @test typeof(f1) <: Vector{Vector{Int}}
            @test typeof(f2) <: Vector{UnitRange{Int}}
            @test length(unique(vcat(f2...))) == n
            for (i1,i2) in zip(f1,f2)
                @test length(unique(union(i2,i1))) == n
            end
        end
    end
end

@testset "leaveout with Int" begin
    @test typeof(leaveout) <: Function
    @test typeof(leaveout(10)) <: Tuple

    @test_throws ArgumentError leaveout(10,10)
    @test_throws ArgumentError leaveout(10,6)
    @test_throws ArgumentError leaveout(10,-1)
    @test_throws ArgumentError leaveout(1)
    @test_throws ArgumentError leaveout(-1)

    for n in (5,10,15,20,50), s in (1,2)
        for (f1,f2) in (@inferred(leaveout(n,s)), @inferred(leaveout(n)))
            @test typeof(f1) <: Vector{Vector{Int}}
            @test typeof(f2) <: Vector{UnitRange{Int}}
            @test length(unique(vcat(f2...))) == n
            for (i1,i2) in zip(f1,f2)
                @test length(unique(union(i2,i1))) == n
            end
        end
    end
end

println("<HEARTBEAT>")

@testset "FoldsView constructor" begin
    @test FoldsView <: AbstractVector
    @test FoldsView <: LearnBase.DataView
    @test FoldsView{Tuple} <: LearnBase.DataView{Tuple}

    @test_reference "references/FoldsView.txt" @io2str(show(::IO, MIME"text/plain"(), kfolds(rand(10),k=5))) by=matrix_compat_isequal

    @testset "Illegal arguments" begin
        # fold indices out of bounds for the given data
        tx = rand(2,15)
        ty = rand(15)
        @test_throws DimensionMismatch FoldsView(tx, kfolds(16)...)
        @test_throws DimensionMismatch FoldsView((tx,rand(14)), kfolds(15)...)
        # number of folds must match for train and test
        f1,tmp = kfolds(15,3)
        tmp,f2 = kfolds(15,5)
        @test_throws DimensionMismatch FoldsView(tx, f1, f2)
        @test_throws DimensionMismatch FoldsView((tx,ty), f1, f2)
        # Only accept Int arrays as indices
        f1,f2 = kfolds(15)
        @test_throws ArgumentError FoldsView(tx, map(x->Float64.(x),f1), f2)
        @test_throws ArgumentError FoldsView(tx, f1, map(x->Float64.(x),f2))
        # only one fold should fail (or should it not?)
        @test_throws ArgumentError FoldsView(tx, map(x->x[1:1], kfolds(15,2))...)
    end

    # FIXME: in 0.6 change back to (Xs,ys,vars...,tuples...)
    for var in (Xs, ys, X, Xv, yv, XX, tuples[1])
        f1, f2 = kfolds(nobs(var), 5)
        for fv in (@inferred(FoldsView(var, f1, f2)),
                   @inferred(FoldsView(var, f1, f2, ObsDim.Last())),
                   FoldsView(var, f1, f2, obsdim=:last))
            @test sum(length.(fv.val_indices)) == nobs(var)
            @test fv.data === var
            @test length(fv.train_indices) === 5
            @test length(fv.val_indices)  === 5
            @test fv.train_indices === f1
            @test fv.val_indices  === f2
            for ((tr, te), i1, i2) in zip(fv, f1, f2)
                @test tr == datasubset(var, i1)
                @test te == datasubset(var, i2)
            end
        end

        f1, f2 = kfolds(nobs(var), 10)
        fv = @inferred(FoldsView(var, f1, f2))
        @test fv.data === var
        @test length(fv.train_indices) === 10
        @test length(fv.val_indices)  === 10
        @test fv.train_indices === f1
        @test fv.val_indices  === f2
        for ((tr, te), i1, i2) in zip(fv, f1, f2)
            @test tr == datasubset(var, i1)
            @test te == datasubset(var, i2)
        end

        f1, f2 = kfolds(nobs(var), 20)
        fv = @inferred(FoldsView(var, f1, f2))
        cumobs = length(fv.val_indices[1])
        for i = 2:length(fv.val_indices)
            cumobs += length(fv.val_indices[i])
            @test 7 <= length(fv.val_indices[i]) <= length(fv.val_indices[1]) <= 8
            @test 1 <= minimum(fv.val_indices[i-1]) < minimum(fv.val_indices[i]) < nobs(var)
        end
        @test cumobs == nobs(var)

        println("<HEARTBEAT>")
    end

    # check that the stored obsdim is correct
    f1, f2 = kfolds(nobs(X), 5)
    # FIXME: in 0.6 change back to (Xs,ys,vars...)
    for var in (Xs, ys, Xv, yv)
        fv = FoldsView(var, f1, f2)
        @test fv.obsdim == ObsDim.Last()
    end
    # FIXME: in 0.6 change back to (tuples...)
    for var in (tuples[1], tuples[2])
        fv = FoldsView(var, f1, f2)
        @test typeof(fv.obsdim) <: Tuple
        @test all(map(x->typeof(x)<:ObsDim.Last, fv.obsdim))
    end
end

println("<HEARTBEAT>")

@testset "FoldsView getindex, endof, length" begin
    f1, f2 = kfolds(nobs(X), 10)
    for var in vars
        fv = FoldsView(var, f1, f2)
        @test length(fv) == 10
        @test fv[end] == fv[length(fv)]
    end
    for var in (Xs, ys, vars...)
        fv = FoldsView(var, f1, f2)
        @test length(fv) == 10
        @test getobs.(fv[end]) == getobs.(fv[length(fv)])
        println("<HEARTBEAT>")
    end
    f1, f2 = kfolds(nobs(X))
    fv = FoldsView(X, f1, f2)
    @test typeof(@inferred(fv[1])) <: Tuple
    @test typeof(@inferred(fv[2:3])) <: FoldsView
    @test @inferred(getobs(fv,1)) == fv[1]
    @test @inferred(getobs(fv[2:3])) == @inferred(getindex(getobs(fv), 2:3))
    @test @inferred(getobs(fv[2:3])) == [fv[2], fv[3]]
    @test typeof(fv[1][1]) <: SubArray
    @test typeof(fv[1][2]) <: SubArray
    @test size(fv[1][1]) == (4,120)
    @test size(fv[1][2]) == (4,30)
    fv = FoldsView(X', f1, f2, obsdim=1)
    @test typeof(@inferred(fv[1])) <: Tuple
    @test typeof(fv[1][1]) <: SubArray
    @test typeof(fv[1][2]) <: SubArray
    @test size(fv[1][1]) == (120,4)
    @test size(fv[1][2]) == (30,4)
    fv = FoldsView((X',X), f1, f2, obsdim=(1,2))
    @test typeof(@inferred(fv[1])) <: Tuple
    @test typeof(fv[1][1]) <: Tuple
    @test typeof(fv[1][2]) <: Tuple
    @test typeof(fv[1][1][1]) <: SubArray
    @test typeof(fv[1][1][2]) <: SubArray
    @test typeof(fv[1][2][1]) <: SubArray
    @test typeof(fv[1][2][2]) <: SubArray
    @test size(fv[1][1][1]) == (120,4)
    @test size(fv[1][1][2]) == (4,120)
    @test size(fv[1][2][1]) == (30,4)
    @test size(fv[1][2][2]) == (4,30)
end

println("<HEARTBEAT>")

@testset "FoldsView iteration" begin
    # FIXME: in 0.6 change back to (X, Xv, yv, XX, XXX, y)
    for var in (X, Xv, yv, XX, y)
        all_train_indices = Array{Int,1}()
        all_test_indices = Array{Int,1}()
        for (train, test) in FoldsView(var, kfolds(nobs(X),10)...)
            @test nobs(train) == 135
            @test nobs(test)  == 15

            @test typeof(train) <: SubArray
            @test typeof(test) <: SubArray

            @test length(setdiff(train.indices[ndims(train)], test.indices[ndims(test)])) == 135
            @test length(setdiff(test.indices[ndims(test)], train.indices[ndims(train)])) == 15
            append!(all_train_indices, train.indices[ndims(train)])
            append!(all_test_indices, test.indices[ndims(test)])
        end
        @test length(unique(all_train_indices)) == 150
        @test length(unique(all_test_indices)) == 150
        @test length(all_test_indices) == 150
        println("<HEARTBEAT>")
    end
    for var in (Xs, ys)
        all_train_indices = Array{Int,1}()
        all_test_indices = Array{Int,1}()
        for (train, test) in FoldsView(var, kfolds(nobs(X),10)...)
            @test nobs(train) == 135
            @test nobs(test)  == 15

            @test typeof(train) <: DataSubset
            @test typeof(test) <: DataSubset

            @test length(setdiff(train.indices, test.indices)) == 135
            @test length(setdiff(test.indices, train.indices)) == 15

            append!(all_train_indices, train.indices)
            append!(all_test_indices, test.indices)
        end
        @test length(unique(all_train_indices)) == 150
        @test length(unique(all_test_indices)) == 150
        @test length(all_test_indices) == 150
    end
    for var in (Xs, ys)
        for (train, test) in FoldsView(var, kfolds(nobs(X),10)...)
            @test nobs(train) == 135
            @test nobs(test)  == 15

            @test typeof(train) <: DataSubset
            @test typeof(test) <: DataSubset
        end
    end
end

println("<HEARTBEAT>")

@testset "kfolds" begin
    # FIXME: in 0.6 change back to (Xs,ys,vars...,tuples...)
    for var in (Xs, ys, Xv, yv, X, tuples[2])
        for kf in (@inferred(kfolds(var)),
                   @inferred(kfolds(var,ObsDim.Last())),
                   @inferred(kfolds(var,5)),
                   @inferred(kfolds(var,5,ObsDim.Last())),
                   kfolds(var,k=5,obsdim=:last))
            @test typeof(kf) <: FoldsView
            @test length(kf.train_indices) == length(kf.val_indices) == 5
            @test kf.data === var
            @test sum(length.(kf.val_indices)) == nobs(var)
        end
    end
    println("<HEARTBEAT>")
    for v1 in (X, Xv), v2 in (y, yv)
        kf = kfolds((v1, v2), k = 15)
        @test typeof(kf) <: FoldsView
        for ((train_x, train_y), (test_x, test_y)) in kf
            @test typeof(train_x) <: SubArray{Float64, 2}
            @test typeof(test_x) <: SubArray{Float64, 2}
            @test typeof(train_y) <: SubArray{String, 1}
            @test typeof(test_y) <: SubArray{String, 1}
            @test nobs(train_x) == 140
            @test nobs(train_y) == 140
            @test nobs(test_x) == 10
            @test nobs(test_y) == 10
        end
    end
    for var in (Xs, ys, vars...)
        kf = kfolds(var)
        @test typeof(kf) <: FoldsView
    end
end

println("<HEARTBEAT>")

@testset "leaveout" begin
    # FIXME: in 0.6 change back to (Xs,ys,vars...,tuples...)
    for var in (Xs, ys, Xv, yv, tuples[1])
        @test_throws ArgumentError leaveout(var,150)
        for kf in (@inferred(leaveout(var)),
                   @inferred(leaveout(var,ObsDim.Last())),
                   @inferred(leaveout(var,1)),
                   @inferred(leaveout(var,1,ObsDim.Last())),
                   leaveout(var,size=1,obsdim=:last))
            @test typeof(kf) <: FoldsView
            @test length(kf.train_indices) == length(kf.val_indices) == 150
            @test kf.data === var
            @test sum(length.(kf.val_indices)) == nobs(var)
        end
        for kf in (@inferred(leaveout(var,30)),
                   @inferred(leaveout(var,30,ObsDim.Last())),
                   leaveout(var,size=30,obsdim=:last))
            @test typeof(kf) <: FoldsView
            @test length(kf.train_indices) == length(kf.val_indices) == 5
            @test kf.data === var
            @test sum(length.(kf.val_indices)) == nobs(var)
        end
    end
    println("<HEARTBEAT>")
    for v1 in (X, Xv), v2 in (y, yv)
        @test length(leaveout((v1, v2)).val_indices) == 150
        kf = leaveout((v1, v2), size = 10)
        @test typeof(kf) <: FoldsView
        for ((train_x, train_y), (test_x, test_y)) in kf
            @test typeof(train_x) <: SubArray{Float64, 2}
            @test typeof(test_x) <: SubArray{Float64, 2}
            @test typeof(train_y) <: SubArray{String, 1}
            @test typeof(test_y) <: SubArray{String, 1}
            @test nobs(train_x) == 140
            @test nobs(train_y) == 140
            @test nobs(test_x) == 10
            @test nobs(test_y) == 10
        end
    end
    for var in (Xs, ys, vars...)
        kf = @inferred leaveout(var)
        @test typeof(kf) <: FoldsView
    end
end

println("<HEARTBEAT>")

@testset "nest DataView" begin
    # FIXME: in 0.6 change back to (Xs,ys,vars...,tuples...)
    for var in (Xs, ys, Xv, yv, tuples[2])
        A = ObsView(var)
        kf = @inferred kfolds(A,10)
        @test kf == kfolds(var,10)
        @test typeof(kf) == typeof(kfolds(var,10))

        A = BatchView(var,15)
        kf = @inferred kfolds(A,10)
        @test kf == kfolds(var,10)
        @test typeof(kf) == typeof(kfolds(var,10))

        A = BatchView(obsview(var),15)
        kf = @inferred kfolds(A,10)
        @test kf == kfolds(var,10)
        @test typeof(kf) == typeof(kfolds(var,10))
    end
end


text:
@recipe function f(sim::ConvergenceSimulation)
  if ndims(collect(values(sim.errors))[1])>1 #Monte Carlo
    vals = [mean(x,1)' for x in values(sim.errors)]
  else #Deterministic
    vals = [x for x in values(sim.errors)]
  end
  seriestype --> :path
  label  --> reshape([string(key) for key in keys(sim.errors)],1,length(keys(sim.errors)))
  xguide  --> "Convergence Axis"
  yguide  --> "Error"
  xscale --> :log10
  yscale --> :log10
  sim.convergence_axis, vals
end

@recipe function f(shoot::Shootout)
  seriestype --> :bar
  legend := false
  xguide --> "Algorithms"
  yguide --> "Efficiency"
  shoot.names,shoot.effs
end

@recipe function f(wp::WorkPrecision)
  seriestype --> :path
  label -->  wp.name
  lw --> 3
  yguide --> "Time (s)"
  xguide --> "Error"
  xscale --> :log10
  yscale --> :log10
  wp.errors,wp.times
end

@recipe function f(wp_set::WorkPrecisionSet)
  seriestype --> :path
  label -->  reshape(wp_set.names,1,length(wp_set))
  lw --> 3
  yguide --> "Time (s)"
  xguide --> "Error"
  xscale --> :log10
  yscale --> :log10
  errors = Vector{Any}(0)
  times = Vector{Any}(0)
  for i in 1:length(wp_set)
    push!(errors,wp_set[i].errors)
    push!(times,wp_set[i].times)
  end
  errors,times
end

@recipe function f(tab::ODERKTableau;dx=1/100,dy=1/100,xlim=[-6,1],ylim=[-5,5])
  x = xlim[1]:dx:xlim[2]
  y = ylim[1]:dy:ylim[2]
  f = (u,v)-> abs(stability_region(u+v*im,tab))<1
  seriestype --> :contour
  fill --> true
  cbar --> false
  x,y,f
end

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
# # 2D forced-dissipative turbulence
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb).
#
# A simulation of forced-dissipative two-dimensional turbulence. We solve the
# two-dimensional vorticity equation with stochastic excitation and dissipation in
# the form of linear drag and hyperviscosity.

using FourierFlows, Printf, Plots

using FourierFlows: parsevalsum
using Random: seed!
using FFTW: irfft

import GeophysicalFlows.TwoDNavierStokes
import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy

# ## Choosing a device: CPU or GPU

dev = CPU()    # Device (CPU/GPU)
nothing # hide

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

 n, L  = 256, 2π             # grid resolution and domain length
 ν, nν = 2e-7, 2             # hyperviscosity coefficient and hyperviscosity order
 μ, nμ = 1e-1, 0             # linear drag coefficient
    dt = 0.005               # timestep
nsteps = 4000                # total number of steps
 nsubs = 20                  # number of steps between each plot
nothing # hide

# ## Forcing

# We force the vorticity equation with stochastic excitation that is delta-correlated
# in time and while spatially homogeneously and isotropically correlated. The forcing
# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and
# width $\delta k_f$, and it injects energy per unit area and per unit time equal
# to $\varepsilon$.

forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
forcing_bandwidth  = 1.5     # the width of the forcing spectrum
ε = 0.1                      # energy input rate by the forcing

grid = TwoDGrid(dev, n, L)

K = @. sqrt(grid.Krsq)
forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
forcing_spectrum[K .< ( 2 * 2π/L)] .= 0 # no power at low wavenumbers
forcing_spectrum[K .> (20 * 2π/L)] .= 0 # no power at high wavenumbers
ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
@. forcing_spectrum *= ε/ε0             # normalize forcing to inject energy at rate ε

seed!(1234)
nothing # hide

# Next we construct function `calcF!` that computes a forcing realization every timestep
function calcF!(Fh, sol, t, clock, vars, params, grid)
  ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
  ξ[1, 1] = 0
  @. Fh = ξ * sqrt(forcing_spectrum)
  nothing
end
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. The
# `stepper` keyword defines the time-stepper to be used.
prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4",
                                calcF=calcF!, stochastic=true)
nothing # hide

# Define some shortcuts for convenience.
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

x, y = grid.x, grid.y
nothing # hide

# First let's see how a forcing realization looks like. Function `calcF!()` computes
# the forcing in Fourier space and saves it into variable `vars.Fh`, so we first need to
# go back to physical space.
calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)

heatmap(x, y, irfft(vars.Fh, grid.nx),
     aspectratio = 1,
               c = :balance,
            clim = (-200, 200),
           xlims = (-L/2, L/2),
           ylims = (-L/2, L/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "a forcing realization",
      framestyle = :box)

# ## Setting initial conditions

# Our initial condition is a fluid at rest.
TwoDNavierStokes.set_zeta!(prob, zeros(grid.nx, grid.ny))

# ## Diagnostics

# Create Diagnostics; the diagnostics are aimed to probe the energy budget.
E  = Diagnostic(energy,                prob, nsteps=nsteps) # energy
Z  = Diagnostic(enstrophy,             prob, nsteps=nsteps) # enstrophy
diags = [E, Z] # a list of Diagnostics passed to `stepforward!` will  be updated every timestep.
nothing # hide

# ## Visualizing the simulation

# We initialize a plot with the vorticity field and the time-series of
# energy and enstrophy diagnostics. To plot energy and enstrophy on the same
# axes we scale enstrophy with $k_f^2$.

p1 = heatmap(x, y, vars.zeta,
         aspectratio = 1,
                   c = :balance,
                clim = (-40, 40),
               xlims = (-L/2, L/2),
               ylims = (-L/2, L/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "vorticity, t="*@sprintf("%.2f", clock.t),
          framestyle = :box)

p2 = plot(2, # this means "a plot with two series"
               label = ["energy E(t)" "enstrophy Z(t) / k_f²"],
              legend = :right,
           linewidth = 2,
               alpha = 0.7,
              xlabel = "μ t",
               xlims = (0, 1.1 * μ * nsteps * dt),
               ylims = (0, 0.55))

l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout = l, size = (900, 420))

# ## Time-stepping the `Problem` forward

# Finally, we time-step the `Problem` forward in time.

startwalltime = time()
anim = @animate for j = 0:Int(nsteps/nsubs)

  log = @sprintf("step: %04d, t: %d, E: %.4f, Z: %.4f, walltime: %.2f min",
      clock.step, clock.t, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

  if j%(1000/nsubs)==0; println(log) end

  p[1][1][:z] = vars.zeta
  p[1][:title] = "vorticity, μ t = "*@sprintf("%.2f", μ * clock.t)
  push!(p[2][1], μ * E.t[E.i], E.data[E.i])
  push!(p[2][2], μ * Z.t[Z.i], Z.data[Z.i] / forcing_wavenumber^2)

  stepforward!(prob, diags, nsubs)
  TwoDNavierStokes.updatevars!(prob)

end

mp4(anim, "twodturb_forced.mp4", fps=18)


text:
using CUDA, SparseArrays

P   = Symbol[Symbol("P$(i)") for i in 1:num_places]
T   = Symbol[Symbol("T$(i)") for i in 1:num_transitions]
C   = sprand(Bool, num_places, num_transitions, 0.3)
ΔF  = sprand(Int, num_places, num_transitions, 0.3)
eft = zeros(Float64, num_transitions)
lft = ones(Float64, num_transitions)
m₀  = zeros(Int, num_places)
h₀  = zeros(Float64, num_transitions)


text:
using RayTracer
using Test

@testset "tuple_creation" begin
    @testset "w_1_is_point" begin
        a = RayTracer.tuple(4.3,-4.2,3.1,1)
        @test a.x ≈ 4.3
        @test a.y ≈ -4.2
        @test a.z ≈ 3.1
        @test a.w == 1
        @test RayTracer.is_point(a) == true
        @test RayTracer.is_vector(a) == false
    end

    @testset "w_0_is_vector" begin
        a = RayTracer.tuple(4.3,-4.2,3.1,0)
        @test a.x ≈ 4.3
        @test a.y ≈ -4.2
        @test a.z ≈ 3.1
        @test a.w == 0
        @test RayTracer.is_point(a) == false
        @test RayTracer.is_vector(a) == true
    end

    @testset "point_creates_tuple" begin
        p = RayTracer.point(4,-4,3)
        @test p == RayTracer.tuple(4,-4,3,1)
    end

    @testset "vector_creates_tuple" begin
        v = RayTracer.vector(4,-4,3)
        @test v == RayTracer.tuple(4,-4,3,0)
    end
end

@testset "tuple_operations" begin
    @testset "adding_tuples" begin
        a1 = RayTracer.tuple(3,-2,5,1)
        a2 = RayTracer.tuple(-2,3,1,0)
        @test a1 + a2 == RayTracer.tuple(1,1,6,1)
    end

    @testset "subtracting_two_points" begin
        p1 = RayTracer.point(3,2,1)
        p2 = RayTracer.point(5,6,7)
        @test p1 - p2 == RayTracer.vector(-2,-4,-6)
    end

    @testset "subtracting_vector_from_point" begin
        p = RayTracer.point(3,2,1)
        v = RayTracer.vector(5,6,7)
        @test p - v == RayTracer.point(-2,-4,-6)
    end

    @testset "subtracting_two_vectors" begin
        v1 = RayTracer.vector(3,2,1)
        v2 = RayTracer.vector(5,6,7)
        @test v1 - v2 == RayTracer.vector(-2,-4,-6)
    end

    @testset "subtracting_vector_from_zero_vector" begin
        zero = RayTracer.vector(0,0,0)
        v = RayTracer.vector(1,-2,3)
        @test zero - v == RayTracer.vector(-1,2,-3)
    end

    @testset "negating_tuple" begin
        a = RayTracer.tuple(1,-2,3,-4)
        @test -a == RayTracer.tuple(-1,2,-3,4)
    end

    @testset "multiplying_tuple_by_scalar" begin
        a = RayTracer.tuple(1,-2,3,-4)
        @test a * 3.5 ≈ RayTracer.tuple(3.5,-7,10.5,-4)
        @test a * 0.5 ≈ RayTracer.tuple(0.5,-1,1.5,-4)
    end

    @testset "dividing_tuple_by_scalar" begin
        a = RayTracer.tuple(1,-2,3,-4)
        @test a / 2 ≈ RayTracer.tuple(0.5,-1,1.5,-4)
    end

    @testset "modulus_unit_vectors" begin
        v1 = RayTracer.vector(1,0,0)
        v2 = RayTracer.vector(0,1,0)
        v3 = RayTracer.vector(0,0,1)
        v4 = RayTracer.vector(1,2,3)
        v5 = RayTracer.vector(-1,-2,-3)
        @test RayTracer.modulus(v1) ≈ 1
        @test RayTracer.modulus(v2) ≈ 1
        @test RayTracer.modulus(v3) ≈ 1
        @test RayTracer.modulus(v4) ≈ √14
        @test RayTracer.modulus(v5) ≈ √14
    end

    @testset "normalization_vectors" begin
        v1 = RayTracer.vector(4,0,0)
        v2 = RayTracer.vector(1,2,3)
        @test RayTracer.normalize(v1) ≈ RayTracer.vector(1,0,0)
        @test RayTracer.normalize(v2) ≈ RayTracer.vector(1/√14,2/√14,3/√14)
    end

    @testset "modulus_normalized_vector" begin
        v = RayTracer.vector(1,2,3)
        norm = RayTracer.normalize(v)
        @test RayTracer.modulus(norm) ≈ 1
    end

    @testset "dot_product" begin
        a = RayTracer.vector(1,2,3)
        b = RayTracer.vector(2,3,4)
        @test RayTracer.dot(a,b) ≈ 20
    end

    @testset "cross_product" begin
        a = RayTracer.vector(1,2,3)
        b = RayTracer.vector(2,3,4)
        @test RayTracer.cross(a,b) ≈ RayTracer.vector(-1,2,-1)
        @test RayTracer.cross(b,a) ≈ RayTracer.vector(1,-2,1)
    end
end

@testset "colors" begin
    @testset "clors_are_tuples" begin
        c = RayTracer.color(-0.5,0.4,1.7)
        @test c.red ≈ -0.5
        @test c.green ≈ 0.4
        @test c.blue ≈ 1.7
    end
    @testset "adding colors" begin
        c1 = RayTracer.color(0.9,0.6,0.75)
        c2 = RayTracer.color(0.7,0.1,0.25)
        @test c1 + c2 ≈ RayTracer.color(1.6,0.7,1.0)
    end
    @testset "subtracting_colors" begin
        c1 = RayTracer.color(0.9,0.6,0.75)
        c2 = RayTracer.color(0.7,0.1,0.25)
        @test c1 - c2 ≈ RayTracer.color(0.2,0.5,0.5)
    end
    @testset "multiplying_color_by_scalar" begin
        c = RayTracer.color(0.2,0.3,0.4)
        @test c * 2 ≈ RayTracer.color(0.4,0.6,0.8)
        @test 2 * c ≈ RayTracer.color(0.4,0.6,0.8)
    end
    @testset "multiplying_colors" begin
        c1 = RayTracer.color(1,0.2,0.4)
        c2 = RayTracer.color(0.9,1,0.1)
        @test c1 * c2 ≈ RayTracer.color(0.9,0.2,0.04)
    end
end

@testset "reflection_at_45_degrees" begin
    v = RayTracer.vector(1,-1,0)
    n = RayTracer.vector(0,1,0)
    r = RayTracer.reflect(v,n)
    @test r ≈ RayTracer.vector(1,1,0)
end

@testset "reflection_slanted_surface" begin
    v = RayTracer.vector(0,-1,0)
    n = RayTracer.vector(√2/2,√2/2,0)
    r = RayTracer.reflect(v,n)
    @test ≈(r,RayTracer.vector(1,0,0),1e-10)
end


text:
@testset "SizedMemory" begin
    for str in ["Hello, world!", "", "αβγ", "    "]
        mem = SizedMemory(str)
        @test pointer(mem) == pointer(str)
        @test length(mem) == ncodeunits(str)

        arr = collect(codeunits(str))
        mem = SizedMemory(arr)
        @test pointer(mem) == pointer(arr)
        @test length(mem) == length(arr)
    end

    for str in ["Hello, world!", "αβγδ", "          "]
        for start in (1, 3), len in (0, 3, 5)
            substr = SubString(str, start:start + len - 1)
            mem = SizedMemory(substr)
            @test pointer(mem) == pointer(substr)
            @test length(mem) == ncodeunits(substr)
        end
    end
end

@testset "memchr" begin
    for T in (nothing, ScanByte.v128, ScanByte.v256)
        for byte in [0x00, 0xa0, 0xda, 0xff]
            @test ScanByte._memchr(T, SizedMemory(UInt8[]), byte) === nothing
            not_bytes = [i for i in 0x00:0xff if i != byte]
            for (len, ind) in [(2, 2), (20, 14), (500, 431)]
                v = rand(not_bytes, len)
                @test ScanByte._memchr(T, SizedMemory(v), byte) === nothing
                v[ind] = byte
                @test ScanByte._memchr(T, SizedMemory(v), byte) == ind
            end
        end
    end

    # Also test memchr directly (instead of _memchr)
    @test memchr(SizedMemory([]), 0x01) === nothing
    @test memchr(SizedMemory([1,2,3]), 0x02) === 9
    @test memchr(SizedMemory([1,2,3]), 0x04) === nothing

    # Test the memchr method with pointer
    bytes = [1,2,3]
    @test memchr(pointer(bytes), UInt(11), 0x02) === 9
    @test memchr(pointer(bytes), UInt(11), 0x03) === nothing
end

@testset "Scanning" begin
    # We test the inverted bytesets, because the codegen is built to look for
    # the bytes its NOT scanning for. Just to make the test easier
    for T in (ScanByte.v128, ScanByte.v256, nothing)
        for inv_byteset in [
            ByteSet(), # empty
            ByteSet(0x00:0xff), # full
            ByteSet(0x41), # one member
            ~ByteSet(0x6a), # all but one member
            ByteSet([0xa2, 0xb4, 0xc8, 0xf9, 0xf6]), # unique lower nibble
            ~ByteSet([0x38, 0x40, 0x51, 0x79, 0x94]), # inv uniq nibble
            ByteSet(0x49:0x71), # Contiguous
            ~ByteSet(0x02:0x31), # Inverted contiguous
            ByteSet(rand(0x80:0xff, 20)), # 128 - 255
            ByteSet(rand(0x00:0x7f, 20)), # 0 - 127
            ByteSet(rand(0x5a:0x94, 20)), # within 127 of each other
            ~ByteSet(rand(0x5a:0x94, 20)), # inv within 127
            ByteSet([0x38, 0x40, 0x90, 0xba, 0xc5, 0xc8, 0xe7]), # at most 8 elements
            ByteSet([i for i in 0x00:0xff if rand(Bool)]) # fallback
        ]
            byteset = ~inv_byteset

            # Empty
            @test ScanByte._memchr(T, SizedMemory(UInt8[]), Val(byteset)) === nothing

            # 1000 bytes not in the set
            if length(byteset) != 256
                bytes = rand(collect(~byteset), 1000)
                @test ScanByte._memchr(T, SizedMemory(bytes), Val(byteset)) === nothing

                if !isempty(byteset)
                    bytes[500] = first(iterate(byteset))
                    @test ScanByte._memchr(T, SizedMemory(bytes), Val(byteset)) == 500

                    bytes[25] = first(iterate(byteset))
                    @test ScanByte._memchr(T, SizedMemory(bytes), Val(byteset)) == 25
                end

                resize!(bytes, 20)
                @test ScanByte._memchr(T, SizedMemory(bytes), Val(byteset)) === nothing
            end

            # Short vector of bytes in the set
            if !isempty(byteset)
                bytes = rand(collect(byteset), 128)
                @test ScanByte._memchr(T, SizedMemory(bytes), Val(byteset)) == 1
            end
        end
    end

    # Also test gen_scan_function directly (instead of _gen_scan_function)
    byteset = ByteSet([0x02, 0x09, 0x11])
    @test memchr(SizedMemory([1,2,3]), Val(byteset)) == 9
    @test memchr(SizedMemory([1,3,4]), Val(byteset)) === nothing
    @test memchr(SizedMemory([1,5,6,7,4,17,13]), Val(byteset)) == 41

    # Test the memchr method with pointer
    bytes = [1,5,6,7,4,17,13]
    @test memchr(pointer(bytes), UInt(45), Val(byteset)) === 41
    @test memchr(pointer(bytes), UInt(33), Val(byteset)) === nothing
end


text:
#
###
# Lagrange polynomial matrices
###

"""
    Compute the Lagrange interpolation matrix from xi to xo.
    lagrange_poly_interp_mat(xₒ,xᵢ)
"""
function lagrange_interp_mat(xo,xi)

    no = length(xo)
    ni = length(xi)

    a = ones(1,ni)
    for i=1:ni
        for j=1:(i-1)  a[i]=a[i]*(xi[i]-xi[j]); end
        for j=(i+1):ni a[i]=a[i]*(xi[i]-xi[j]); end
    end
    a = 1 ./ a # Barycentric weights

    J = zeros(no,ni)
    s = ones(1,ni)
    t = ones(1,ni)
    for i=1:no
        x = xo[i]
        for j=2:ni
            s[j]      = s[j-1]    * (x-xi[j-1]   )
            t[ni+1-j] = t[ni+2-j] * (x-xi[ni+2-j])
        end
        J[i,:] = a .* s .* t
    end

    return J
end

"""
 Compute derivative matrix for lagrange
 interpolants on points [x]
"""
function lagrange_deriv_mat(x)

    n = length(x)

    a = ones(1,n)
    for i=1:n
        for j=1:(i-1) a[i]=a[i]*(x[i]-x[j]) end
        for j=(i+1):n a[i]=a[i]*(x[i]-x[j]) end
    end
    a = 1 ./ a # Barycentric weights

    # diagonal elements
    D = x .- x'
    for i=1:n D[i,i] = 1. end
    D = 1 ./ D
    for i=1:n
        D[i,i] = 0.
        D[i,i] = sum(D[i,:])
    end

    # off-diagonal elements
    for j=1:n for i=1:n
        if(i!=j) D[i,j] = a[j] / (a[i]*(x[i]-x[j])) end
    end end

    return D
end


text:
using BenchmarkTools, IntervalArithmetic

# don't really want to measure interval performance, so we're using the
# fastest interval rounding mode
setrounding(Interval, :none)
using McCormick

println("Begin Allocation Check")

# Allocations are expected for the following operations:
# - MC(Float64, Interval, 1): Nonstatic definition of static type.
# - exp2, exp10, tanh, atanh, acosh: Underlying interval operation may allocate.
for T in (NS, Diff, MV)
    !iszero(allocs(@benchmark MC{$5,$T}(2.1) samples = 1)) && println("MC(Float64) Allocates")
    !iszero(allocs(@benchmark MC{$5,$T}(0.1,2.0) samples = 1)) && println("MC(Float64, Float64) Allocates")
    !iszero(allocs(@benchmark MC{$5,$T}((Interval{Float64}(0.1,2.0))) samples = 1)) && println("MC(Interval) Allocates")
    a = MC{5,T}(1.0,(Interval{Float64}(0.1,2.0)),2)
    for op in (+, -, exp, expm1, log, log2, log10, log1p,
               abs, one, zero, real, max, min, inv, cosh, sqrt,
               isone, isnan)
        !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("$op Allocates")
    end
           # domain violations asin, asind, acos, acosd, atan, atand, tan, tand, sec, secd,
           for op in (sin, sind, cos, cosd, step, sign, sinh)
               a = MC{5,T}(0.4,(Interval{Float64}(0.1,0.9)),2)
               !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("Positive $op Allocates")
               a = MC{5,T}(0.1,(Interval{Float64}(-0.4,0.3)),2)
               !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("0 in X $op Allocates")
               a = MC{5,T}(-0.5,(Interval{Float64}(-0.9,-0.1)),2)
               !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("Negative $op Allocates")
              end
           for op in (min, max, *, -, +, /)
               x = MC{5,T}(0.4,(Interval{Float64}(0.1,0.9)),2)
               y = MC{5,T}(0.5,(Interval{Float64}(0.3,0.9)),2)
               q = 1.34534
               !iszero(allocs(@benchmark ($op)($x, $q) samples = 1)) && println("$op(x,q) Allocates")
               !iszero(allocs(@benchmark ($op)($q, $x) samples = 1)) && println("$op(q,x) Allocates")
               !iszero(allocs(@benchmark ($op)($x, $y) samples = 1)) && println("$op(x,y) Allocates")
           end
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^2 samples = 1)) && println("a^2 Pos Allocates")
           a = MC{5,T}(0.1,Interval{Float64}(-0.4,0.3),2);   !iszero(allocs(@benchmark $a^2 samples = 1)) && println("a^2 Zero Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^2 samples = 1)) && println("a^2 Neg Allocates")
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^3 samples = 1)) && println("a^3 Pos Allocates")
           a = MC{5,T}(0.1,Interval{Float64}(-0.4,0.3),2);   !iszero(allocs(@benchmark $a^3 samples = 1)) && println("a^3 Zero Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^3 samples = 1)) && println("a^3 Neg Allocates")
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^4 samples = 1)) && println("a^4 Pos Allocates")
           a = MC{5,T}(0.1,Interval{Float64}(-0.4,0.3),2);   !iszero(allocs(@benchmark $a^4 samples = 1)) && println("a^4 Zero Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^4 samples = 1)) && println("a^4 Neg Allocates")
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^(-1) samples = 1)) && println("a^-1 Pos Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^(-1) samples = 1)) && println("a^-1 Neg Allocates")
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^(-2) samples = 1)) && println("a^-2 Pos Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^(-2) samples = 1)) && println("a^-2 Neg Allocates")
           a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^(-3) samples = 1)) && println("a^-3 Pos Allocates")
           a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^(-3) samples = 1)) && println("a^-3 Neg Allocates")
   end
end

println("Allocation Check Complete!")


text:
# Try to approximate exponential-sine-squared kernel:
using PyPlot
using Optim

function expsin2_approx()
nexp = 5
ngamma = 100
coeff_grid = zeros(ngamma,nexp)
#gamma_grid = logspace(-2,0.5,ngamma)
gamma_grid = linspace(0.01,3.0,ngamma)
ntime = 200
time = collect(linspace(0,2pi,ntime))

gamma = 1.0
gaussian=exp(-gamma.*sin(time).^2)

function compute_chi(x)
model = zeros(ntime)
for k=1:nexp
  model += x[k].*cos((k-1)*2.*time)
end
return sum((gaussian-model).^2)
end

npar = nexp

for i=1:ngamma
gamma = gamma_grid[i]
gaussian=exp(-gamma.*sin(time).^2)

abest=zeros(npar)
abest = [0.6,0.2,0.01,0.05,.025,.025]

atrial = zeros(npar)
for k=1:npar
  atrial[k] = abest[k]
end
model = zeros(ntime)
for k=1:nexp
  model += atrial[k].*cos((k-1)*2.*time)
end
clf()
plot(time,gaussian)
plot(time,model)
#chibest=sum((gaussian-model).^2)
chibest=compute_chi(atrial)

println("Initial chi-square: ",chibest)
#read(STDIN,Char)

result = optimize(compute_chi, atrial, BFGS(), OptimizationOptions(autodiff = true))

println(Optim.minimizer(result),Optim.minimum(result))
abest= Optim.minimizer(result)
clf()
plot(time,gaussian)
model = zeros(ntime)
for k=1:nexp
  model += abest[k].*cos((k-1)*2.*time)
end
plot(time,model)
plot(time,(gaussian-model))
chi=sum((gaussian-model).^2)
println("New minimum: ",chi," ",abest," ",std(gaussian-model)," ",maximum(abs(gaussian-model)))
#read(STDIN,Char)
coeff_grid[i,:] = abest
end

#coeff1 = log10(vec(coeff_grid[:,1]))

#function coeff1_fit(x)
## Fit the first coefficient:
#model= x[1]+x[2].*gamma_grid.^x[3]
#return sum((coeff1-model).^2)
#end

#x1 = [0.0,-1.75,0.2]
#result = optimize(coeff1_fit, x1, BFGS(), OptimizationOptions(autodiff = true))
#clf()
#loglog(gamma_grid,coeff_grid[:,1])
#x1 = Optim.minimizer(result)
##loglog(gamma_grid,10.^(-x1[2].*exp(log10(gamma_grid)*x1[1])),"r-")
#loglog(gamma_grid,10.^(x1[1]+x1[2].*gamma_grid.^x1[3]),"r-")
#println(Optim.minimizer(result),Optim.minimum(result))
#println("Fit to first coefficient")
#read(STDIN,Char)

# Now fit the grid of coefficients:
coeff2 = vec(coeff_grid[:,1])

ord = 1.0
function coeff2_fit(x)
# Fit the second-fourth coefficient:
model= x[1]+x[2].*exp(-gamma_grid.^x[3].*x[4])
return sum((coeff2-model).^2)
end

x1 = [0.45,0.4,1.0,1.0]
result = optimize(coeff2_fit, x1, BFGS(), OptimizationOptions(autodiff = true))
x1 = Optim.minimizer(result)
clf()
plot(gamma_grid,coeff2)
plot(gamma_grid,x1[1]+x1[2].*exp(-gamma_grid.^x1[3].*x1[4]))
println(Optim.minimizer(result),Optim.minimum(result))
println("Fit to first coefficient")
read(STDIN,Char)

coeff2 = vec(coeff_grid[:,2])
ord = 1.0
x2 = [0.45,0.4,1.0,1.0]
result = optimize(coeff2_fit, x2, BFGS(), OptimizationOptions(autodiff = true))
x2 = Optim.minimizer(result)
clf()
plot(gamma_grid,coeff2)
plot(gamma_grid,x2[1]+x2[2].*exp(-gamma_grid.^x2[3].*x2[4]))
println(Optim.minimizer(result),Optim.minimum(result))
println("Fit to second coefficient")
read(STDIN,Char)

coeff2 = vec(coeff_grid[:,3])
ord = 2.0
x3 = [0.45,0.4,1.0,1.0]
result = optimize(coeff2_fit, x3, BFGS(), OptimizationOptions(autodiff = true))
x3 = Optim.minimizer(result)
clf()
plot(gamma_grid,coeff2)
plot(gamma_grid,x3[1]+x3[2].*exp(-gamma_grid.^x3[3].*x3[4]))
println(Optim.minimizer(result),Optim.minimum(result))
println("Fit to third coefficient")
read(STDIN,Char)

coeff2 = vec(coeff_grid[:,4])
ord = 3.0
x4 = [0.45,0.4,1.0,1.0]
result = optimize(coeff2_fit, x4, BFGS(), OptimizationOptions(autodiff = true))
x4 = Optim.minimizer(result)
clf()
plot(gamma_grid,coeff2)
plot(gamma_grid,x4[1]+x4[2].*exp(-gamma_grid.^x4[3].*x4[4]))
println(Optim.minimizer(result),Optim.minimum(result))
println("Fit to fourth coefficient")
read(STDIN,Char)

coeff2 = vec(coeff_grid[:,5])
ord = 4.0
x5 = [0.45,0.4,1.0,1.0]
result = optimize(coeff2_fit, x5, BFGS(), OptimizationOptions(autodiff = true))
x5 = Optim.minimizer(result)
clf()
plot(gamma_grid,coeff2)
plot(gamma_grid,x5[1]+x5[2].*exp(-gamma_grid.^x5[3].*x5[4]))
println(Optim.minimizer(result),Optim.minimum(result))
println("Fit to fifth coefficient")
read(STDIN,Char)

# Now test with random numbers:
ntest = 10000
dev = zeros(ntest)
for i=1:ntest
  time = rand()*2pi
  gamma = rand()*3.0
  model = 0.0
  model += x1[1]+x1[2].*exp(-gamma.^x1[3].*x1[4])
  model += (x2[1]+x2[2].*exp(-gamma.^x2[3].*x2[4]))*cos(2*time)
  model += (x3[1]+x3[2].*exp(-gamma.^x3[3].*x3[4]))*cos(4*time)
  model += (x4[1]+x4[2].*exp(-gamma.^x4[3].*x4[4]))*cos(6*time)
  model += (x5[1]+x5[2].*exp(-gamma.^x5[3].*x5[4]))*cos(8*time)
  value = exp(-gamma*sin(time)^2)
  dev[i]=value-model
  println(i," ",gamma," ",time," ",value," ",model," ",value-model)
end
println("Maximum deviation:  ",maximum(abs(dev)))
println("Standard deviation: ",std(dev))
read(STDIN,Char)

return gamma_grid,x1,x2,x3,x4,x5
end

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
@testset "model.query" begin

@testset "equivalence" begin
    @test QueryFormula(:v, :a) == QueryFormula(:v, :a)
    @test QueryFormula(:v, :b) != QueryFormula(:v, :c)
    @test QueryFormula(:x, 1) != QueryFormula(:v, 1)

    and1 = QueryFormula(:and, Set([
        QueryFormula(:x, "hoge"),
        QueryFormula(:y, "fuga"),
        QueryFormula(:z, "piyo")
    ]))
    and2 = QueryFormula(:and, Set([
        QueryFormula(:z, "piyo"),
        QueryFormula(:y, "fuga")
    ]))
    @test and1 == and1
    @test and1 != and2
    @test issubset(and1, and2)  # seems not intuitive, but and1 result records should be less than and2
    @test !issubset(and2, and1)
    @test issubset(and1, QueryFormula(:x, "hoge"))
    @test !issubset(QueryFormula(:x, "hoge"), and1)

    or1 = QueryFormula(:or, Set([
        QueryFormula(:x, "hoge"),
        QueryFormula(:x, "fuga"),
        QueryFormula(:x, "piyo")
    ]))
    or2 = QueryFormula(:or, Set([
        QueryFormula(:x, "piyo"),
        QueryFormula(:x, "fuga")
    ]))
    @test or1 == or1
    @test or1 != or2
    @test issubset(or2, or1)  # opposit to the `and` case
    @test !issubset(or1, or2)

    nested1 = QueryFormula(:and, Set([
        QueryFormula(:x, 1),
        QueryFormula(:or, Set([
            QueryFormula(:v, :a),
            QueryFormula(:v, :b)
        ])),
        QueryFormula(:or, Set([
            QueryFormula(:n, true),
            QueryFormula(:m, true)
        ]))
    ]))
    nested2 = QueryFormula(:and, Set([
        QueryFormula(:x, 1),
        QueryFormula(:or, Set([
            QueryFormula(:v, :a),
            QueryFormula(:v, :b),
            QueryFormula(:v, :c)
        ])),
        QueryFormula(:or, Set([
            QueryFormula(:n, true),
            QueryFormula(:m, true),
            QueryFormula(:o, true)
        ]))
    ]))
    nested3 = QueryFormula(:and, Set([
        QueryFormula(:x, 1),
        QueryFormula(:v, :a),
        QueryFormula(:or, Set([
            QueryFormula(:n, true),
            QueryFormula(:m, true)
        ]))
    ]))
    @test nested1 == nested1
    @test nested1 != nested2
    @test issubset(nested1, nested2)
    @test !issubset(nested2, nested1)
    @test issubset(nested3, nested1)
    @test !issubset(nested1, nested3)

    nota = QueryFormula(:not, QueryFormula(:v, :a))
    notb = QueryFormula(:not, QueryFormula(:v, :b))
    @test nota == nota
    @test nota != notb

    or3 = QueryFormula(:or, Set([
        QueryFormula(:x, 1),
        QueryFormula(:x, 2)
    ]))
    or4 = QueryFormula(:or, Set([
        QueryFormula(:v, 1),
        QueryFormula(:x, 1),
        QueryFormula(:x, 2)
    ]))
    not3 = QueryFormula(:not, QueryFormula(:x, 3))
    @test or3 != not3
    @test issubset(or3, not3)
    @test !issubset(not3, or3)
    @test !issubset(or4, not3)

    rec1 = QueryFormula(:recursive, "[NH2]C")
    rec2 = QueryFormula(:recursive, "[NH2]CC")
    rec3 = QueryFormula(:recursive, "[NH2+0]C")
    rec4 = QueryFormula(:recursive, "[N,O;H2]C")
    and3 = QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :N)
        QueryFormula(:or, Set([
            QueryFormula(:hydrogenconnected, 2),
            QueryFormula(:hydrogenconnected, 3)
        ]))
    ]))
    @test issubset(rec2, rec1)
    @test !issubset(rec1, rec2)
    @test issubset(rec3, rec1)
    @test !issubset(rec1, rec3)
    @test !issubset(rec4, rec1)
    @test issubset(rec1, rec4)
    @test issubset(rec1, and3)
    @test !issubset(rec4, and3)
end

@testset "tidyformula" begin
    and1 = QueryFormula(:and, Set([
        QueryFormula(:x, :a),
        QueryFormula(:and, Set([
            QueryFormula(:y, :b),
            QueryFormula(:z, :c)
        ]))
    ]))
    @test tidyformula(and1) == QueryFormula(:and, Set([
        QueryFormula(:x, :a),
        QueryFormula(:y, :b),
        QueryFormula(:z, :c)
    ]))

    or1 = QueryFormula(:or, Set([
        QueryFormula(:or, Set([
            QueryFormula(:x, :a),
            QueryFormula(:y, :a)
        ])),
        QueryFormula(:z, :a)
    ]))
    @test tidyformula(or1) == QueryFormula(:or, Set([
        QueryFormula(:x, :a),
        QueryFormula(:y, :a),
        QueryFormula(:z, :a)
    ]))

    and2 = QueryFormula(:and, Set([
        QueryFormula(:and, Set([
            QueryFormula(:x, 1),
            QueryFormula(:y, 1)
        ])),
        QueryFormula(:v, :a),
        QueryFormula(:or, Set([
            QueryFormula(:z, 0),
            QueryFormula(:z, nothing)
        ])),
        QueryFormula(:not, QueryFormula(:a, "hoge"))
    ]))
    @test tidyformula(and2) == QueryFormula(:and, Set([
        QueryFormula(:v, :a),
        QueryFormula(:x, 1),
        QueryFormula(:y, 1),
        QueryFormula(:or, Set([
            QueryFormula(:z, 0),
            QueryFormula(:z, nothing)
        ])),
        QueryFormula(:not, QueryFormula(:a, "hoge"))
    ]))

    or2 = QueryFormula(:or, Set([
        QueryFormula(:and, Set([
            QueryFormula(:x, :a),
            QueryFormula(:y, :b)
        ])),
        QueryFormula(:and, Set([
            QueryFormula(:x, :c),
            QueryFormula(:y, :d)
        ])),
        QueryFormula(:or, Set([
            QueryFormula(:x, :e),
            QueryFormula(:y, :f)
        ]))
    ]))
    @test tidyformula(or2) == QueryFormula(:or, Set([
        QueryFormula(:and, Set([
            QueryFormula(:x, :a),
            QueryFormula(:y, :b)
        ])),
        QueryFormula(:and, Set([
            QueryFormula(:x, :c),
            QueryFormula(:y, :d)
        ])),
        QueryFormula(:x, :e),
        QueryFormula(:y, :f)
    ]))

    not1 = QueryFormula(:not, QueryFormula(:v, :a))
    @test tidyformula(not1) == not1

    and3 = QueryFormula(:and, Set([
        QueryFormula(:and, Set([
            QueryFormula(:and, Set([
                QueryFormula(:a, 1),
                QueryFormula(:and, Set([
                    QueryFormula(:b, 1),
                    QueryFormula(:c, 1)
                ]))
            ])),
            QueryFormula(:d, 1)
        ])),
        QueryFormula(:e, 1)
    ]))
    @test tidyformula(and3) == QueryFormula(:and, Set([
        QueryFormula(:a, 1),
        QueryFormula(:b, 1),
        QueryFormula(:c, 1),
        QueryFormula(:d, 1),
        QueryFormula(:e, 1)
    ]))

    abs1 = QueryFormula(:and, Set([
        QueryFormula(:x, "important"),
        QueryFormula(:or, Set([
            QueryFormula(:x, "important"),
            QueryFormula(:z, "optional")
        ]))
    ]))
    @test tidyformula(abs1) == QueryFormula(:x, "important")

    abs2 = QueryFormula(:and, Set([
        QueryFormula(:or, Set([
            QueryFormula(:a, 1),
            QueryFormula(:a, 2),
            QueryFormula(:a, 3)
        ]))
        QueryFormula(:or, Set([
            QueryFormula(:a, 1),
            QueryFormula(:a, 2)
        ]))
    ]))
    @test tidyformula(abs2) == QueryFormula(:or, Set([
        QueryFormula(:a, 1),
        QueryFormula(:a, 2)
    ]))

    dist1 = QueryFormula(:or, Set([
        QueryFormula(:and, Set([
            QueryFormula(:a, 1),
            QueryFormula(:b, 2),
            QueryFormula(:c, 3),
        ])),
        QueryFormula(:and, Set([
            QueryFormula(:a, 1),
            QueryFormula(:b, 2),
            QueryFormula(:c, 4)
        ]))
    ]))
    @test tidyformula(dist1) == QueryFormula(:and, Set([
        QueryFormula(:a, 1),
        QueryFormula(:b, 2),
        QueryFormula(:or, Set([
            QueryFormula(:c, 3)
            QueryFormula(:c, 4)
        ]))
    ]))

    dist2 = QueryFormula(:and, Set([
        QueryFormula(:or, Set([
            QueryFormula(:x, false),
            QueryFormula(:y, nothing)
        ])),
        QueryFormula(:or, Set([
            QueryFormula(:x, false),
            QueryFormula(:z, 3),
            QueryFormula(:z, 4)
        ]))
    ]))
    @test tidyformula(dist2) == QueryFormula(:or, Set([
        QueryFormula(:x, false),
        QueryFormula(:and, Set([
            QueryFormula(:y, nothing)
            QueryFormula(:or, Set([
                QueryFormula(:z, 3),
                QueryFormula(:z, 4)
            ]))
        ]))
    ]))

    dist3 = QueryFormula(:and, Set([
        QueryFormula(:hydrogenconnected, 1),
        QueryFormula(:or, Set([
            QueryFormula(:and, Set([
                QueryFormula(:atomsymbol, :O),
                QueryFormula(:isaromatic, false)
            ])),
            QueryFormula(:and, Set([
                QueryFormula(:atomsymbol, :N),
                QueryFormula(:isaromatic, false)
            ]))
        ]))
    ]))
    @test tidyformula(dist3) == QueryFormula(:and, Set([
        QueryFormula(:hydrogenconnected, 1),
        QueryFormula(:or, Set([
            QueryFormula(:atomsymbol, :O),
            QueryFormula(:atomsymbol, :N)
        ])),
        QueryFormula(:isaromatic, false)
    ]))

    abs3 = QueryFormula(:and, Set([
        QueryFormula(:any, true),
        QueryFormula(:a, 1)
    ]))
    @test tidyformula(abs3) == QueryFormula(:a, 1)

    abs4 = QueryFormula(:or, Set([
        QueryFormula(:any, true),
        QueryFormula(:a, 1)
    ]))
    @test tidyformula(abs4) == QueryFormula(:any, true)

    abs5 = QueryFormula(:and, Set([
        QueryFormula(:any, false),
        QueryFormula(:a, 1)
    ]))
    @test tidyformula(abs5) == QueryFormula(:any, false)

    abs6 = QueryFormula(:or, Set([
        QueryFormula(:any, false),
        QueryFormula(:a, 1)
    ]))
    @test tidyformula(abs6) == QueryFormula(:a, 1)

    abs7 = QueryFormula(:and, Set([
        QueryFormula(:or, Set([
            QueryFormula(:or, Set([
                QueryFormula(:any, true),
                QueryFormula(:c, 3)
            ])),
            QueryFormula(:b, 2)
        ])),
        QueryFormula(:a, 1)
    ]))
    @test tidyformula(abs7) == QueryFormula(:a, 1)
end

@testset "convertnotquery" begin
    not1 = QueryFormula(:not, QueryFormula(:atomsymbol, :H))
    @test convertnotquery(not1) == QueryFormula(:any, true)

    not2 = QueryFormula(:not, QueryFormula(:isringbond, true))
    @test convertnotquery(not2) == QueryFormula(:isringbond, false)

    not3 = QueryFormula(:not, QueryFormula(:sssrcount, 0))
    @test convertnotquery(not3) == QueryFormula(:or, Set([
        QueryFormula(:sssrcount, 1),
        QueryFormula(:sssrcount, 2),
        QueryFormula(:sssrcount, 3)
    ]))

    not4 = QueryFormula(:and, Set([
        QueryFormula(:not, QueryFormula(:hydrogenconnected, 0)),
        QueryFormula(:not, QueryFormula(:hydrogenconnected, 1)),
        QueryFormula(:not, QueryFormula(:hydrogenconnected, 2))
    ]))
    @test tidyformula(convertnotquery(not4)) == QueryFormula(:or, Set([
        QueryFormula(:hydrogenconnected, 3),
        QueryFormula(:hydrogenconnected, 4)
    ]))
end

@testset "querymol" begin
    q = querymol(SmartsAtom, SmartsBond)
    addatom!(q, SmartsAtom(QueryFormula(:any, true)))
    addatom!(q, SmartsAtom(QueryFormula(:any, true)))
    addbond!(q, 1, 2, SmartsBond(QueryFormula(:any, true)))
    @test atomcount(q) == 2
    @test bondcount(q) == 1
    @test hasbond(q, 1, 2)
end

end # model.query


text:
using JuMP
using CPLEX
include("relaxed.jl")
include("base.jl")
include("DW_mp_dual.jl")
include("heuristic.jl")

function test(str)
	Relaxed("di-yuan/di-yuan_1/")
	t=time()
	Relaxed(str)
	t = round(time()-t, digits=2)
	println(t)

	return 0
end

function func()
	dir_tab=["pdh","di-yuan","abilene","atlanta","dfn-bwin","dfn-gwin","newyork","nobel-germany","nobel-us","polska"]

	for dir in dir_tab

		files = readdir("Instances/"*dir)
		for file in files

			f = open("perf.txt","w")

			str=dir*"/"*file*"/"

			println(str)

			write(f,file*" & ")

			t = time()
			val1 = Relaxed(str)[1]
			t = round(time()-t, digits=2)
			write(file,string(t)*" & ")

			t = time()
			val2 = Static(str)[1]
			t = round(time()-t, digits=2)
			write(file,string(t)*" & ")

			if heuristic(str)

				t = time()
				heuristic(str)
				t = round(time()-t, digits=2)
				write(file,string(t)*" & ")

				ov, x_i, x_fi, x_ikf, e = heuristic(str,res=true)
				X_i=Array{Float64,2}(zeros(1,size(x_i,1)))
				X_i[1,:]=x_i
				X_fi=Array{Float64,3}(zeros(1,size(x_fi,1),size(x_fi,2)))
				X_fi[1,:,:]=x_fi
				X_ikf=Array{Float64,4}(zeros(1,size(x_ikf,1),size(x_ikf,2),size(x_ikf,3)))
				X_ikf[1,:,:,:]=x_ikf
				E=Array{Float64,5}(zeros(1,size(e,1),size(e,2),size(e,3),size(e,4)))
				E[1,:,:,:,:]=e

				t = time()
				DW_md(str,X_i,X_fi,X_ikf,E)
				t = round(time()-t, digits=2)
				write(file,string(t)*" & ")
			else
				write(file," \emptyset & \empty set & ")
			end
			write(file,string(val2)*" & "*string(val1)*'\n')

			close(f)

		end
	end

end

func()

max_stars_repo_name: boumpteryx/RORT

max_issues_repo_name: boumpteryx/RORT

max_forks_repo_name: boumpteryx/RORT


text:
## Features.jl
## -----------
## Feature Types and basic functionality
## related directly to IO of Feature Types
using DataFrames
using IntervalTrees

abstract Feature
abstract SequenceFeature <: Feature

type SequencePos
     start::Int32
     stop::Int32

     function SequencePos( start::String,stop::String)
        start_int = parseint(start)
        stop_int  = parseint(stop)
        start     = int32(start_int)
        stop      = int32(stop_int)
        new(start, stop)
    end
    function SequencePos(start::Int,stop::Int)
       new(int32(start),int32(stop))
    end
end

#DEPRECATE
type GenomePos
     chr::Int16
     start::Int32
     stop::Int32

     function GenomePos( chr::String,start::String,stop::String, organism="none")
        start_int = parseint(start)
        stop_int  = parseint(stop)
        start     = int32(start_int)
        stop      = int32(stop_int)
        if organism == "none"
            chr       = chr_convert_int16(chr)
        else
            chr  = chr_convert_int16(chr, organism)
        end
        new(chr, start, stop)
    end
    function GenomePos(chr::Int,start::Int,stop::Int)
       new(int16(chr),int32(start),int32(stop))
    end
    function GenomePos(chr::String, start::Int,stop::Int)
       new( chr_convert_int16(chr),int32(start),int32(stop))
    end
end

type CytosineCount <: SequenceFeature
    pos::SequencePos
    strand::Bool
    C_count::Int16
    T_count::Int16
end

type CytosineMethylation <: SequenceFeature
     pos::SequencePos
     strand::Bool
     mc::Float32
     hmc::Float32
     function CytosineMethylation(pos::SequencePos,strand::Bool, mc,hmc)
       new(pos,strand,float32(mc),float32(hmc))
    end
end

type CytosineCount2
    pos::GenomePos
    strand::Bool
    C_count::Int16
    T_count::Int16
end

type CytosineMethylation2
    pos::GenomePos
    strand::Bool
    mc::Float32
    hmc::Float32
    #constructor to ensure float32s
    function CytosineMethylation(pos::GenomePos,strand::Bool, mc,hmc)
       new(pos,strand,float32(mc),float32(hmc))
    end
end

# This holds features in a dataframe
type FeatureTable
  description::Dict{ASCIIString,Any}
  dataframe::DataFrame
end

# This holds Features organised by sequence and ordered by co-ordinate
# Idea being this might be efficient for holding many features
type SequenceFeatures
    description::Dict{ASCIIString,Any}
    features::Dict{Any,Any}
end

function add_sequence_feature(sequence_features::SequenceFeatures,  sequence_name::String, feature::SequenceFeature )

     if haskey( sequence_features.features, sequence_name )
         sequence_feature_vector = sequence_features.features[sequence_name]
         push!( sequence_feature_vector, feature)
     else
         vec = Array(SequenceFeature,0)
         push!(vec, feature)
         sequence_features.features[sequence_name] = vec
     end
end

# Similar to SequenceFeatures
# Based on an interval tree
type FeatureRegions
  description::Dict{ASCIIString,Any}
  regions::Dict{Any,Any}

  function FeatureRegions(FT::FeatureTable)
     sequence_regions = Dict()
     DF = FT.dataframe
     nrows = size(DF,1)
     for row = 1:nrows
         gene = DF[row, :]
         #println(gene)
         gene_name = gene[:name][1]
         start     = gene[:start][1]
         stop      = gene[:stop][1]
         chr       = gene[:chr][1]
         insertintervals!(sequence_regions, chr, start, stop, gene )
     end
     new( ["organism"=>""], sequence_regions )
  end
end

# I think the reason we can't just stick anything into value
# and we have multiple functions for with different value types is because
# IntervalTree didn't let us do that. Check that.
function insertintervals!( sequence_regions_dict, chr, start, stop, value::DataFrame)
    if haskey(sequence_regions_dict, chr)
        interval_tree = sequence_regions_dict[chr]
        interval_tree[(start,stop)] = value
    else
        try
           interval_tree = IntervalTree{Int32,DataFrame}()
           interval_tree[(start,stop)] = value
           sequence_regions_dict[chr]  = interval_tree
        catch e
           println("Exception chr: $chr, start: $start, stop: $stop, value $value ")
           println( join( [  "Exception, ", typeof(chr), " ", typeof(start), " ", typeof(stop) ], "\t" ) )
           error(e)
        end
    end
end

function insertintervals!( sequence_regions_dict, chr, start, stop, value::ASCIIString)
    if haskey(sequence_regions_dict, chr)
        interval_tree = sequence_regions_dict[chr]
        interval_tree[(start,stop)] = value
    else
        interval_tree = IntervalTree{Int32,ASCIIString}()
        interval_tree[(start,stop)] = value
        sequence_regions_dict[chr] = interval_tree
    end
end

function getoverlaps(feature_regions,chr,start,stop)
     # if (haskey(genome_intervals, chr) == true )
          interval_tree = feature_regions.regions[chr]
          iter = intersect(interval_tree, (start,stop))
          return iter
      #else
      #   println("WARN: No key found for $key")
      #   return []
      #end
end

function chr_convert_int16( chr)
    rgx1 = match(r"chr(.)",chr)
    if rgx1 != nothing
      # if it's just numbers
        chr_name = rgx1.captures[1]
        if chr_name == "X"
            return int16(200)
         end
        if chr_name == "Y"
            return int16(201)
        end
        if chr_name == "M"
            return int16(202)
        end
        rgx2 = match(r"^[0-9]+$", string( chr_name ))
        if rgx2 != nothing
           #chr_num = rgx2.captures[1]
            return convert(Int16, parseint(chr_name))
        else
           error("no valid chr number for $chr")
        end
    else
         error("Chromosome field not in valid format e.g. chr1")
    end
end

function feature_table_to_bed( features::FeatureTable, path, delim="\t"; overwrite=false)
        if isfile(path) && !overwrite
             println("feature_table_to_bed: $path exists")
             return 0;
         end

         df = features.dataframe
         nrows = size(df,1)
         iostream = open(path,"w")
         count = 0
         for row =  1:nrows
         	  row    = df[row,:]
            name   = row[:name][1]
            start  = row[:start][1]
            stop   = row[:stop][1]
            strand = row[:strand][1]
            strand = strand ? "+" : "-"
         	  #chr    = CURRENT_ANNOTATION_CHR_NAMES["mm10"][ row[:chr][1] ]
            chr    = row[:chr][1]
            line   = join([chr,start,stop,name,0,strand,"\n"],delim)
            print(iostream, line)
            count += 1
         end
         close(iostream)
         return count
end


text:
struct NotImplementedError <: Exception
    msg::String
end

function Base.showerror(io::IO, err::NotImplementedError)
    print(io, "NotImplementedError: ")
    print(io, err.msg)
end

"""
    maxdiff(a, b)

Get the maximum difference between vectors. Equivalent to the infinity norm of the
difference, but avoids forming the temporary difference vector.
"""
function maxdiff(a, b)
    err = -Inf
    for i = min(length(a), length(b))
        err = max(err, abs(a[i] - b[i]))
    end
    return err
end

# Add non-allocating lu! method that allows you to pass in the pivot vector
using LinearAlgebra: BlasFloat, BlasInt
using LinearAlgebra.BLAS: @blasfunc
@static if VERSION < v"1.7"
    const libblas = LinearAlgebra.LAPACK.liblapack
else
    const libblas = LinearAlgebra.LAPACK.libblastrampoline
end
for (getrf,elty) in ((:dgetrf_, :Float64), (:sgetrf, :Float32))
    @eval begin
        function LinearAlgebra.LAPACK.getrf!(A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt})
            LinearAlgebra.require_one_based_indexing(A)
            LinearAlgebra.chkstride1(A)
            m, n = size(A)
            @assert length(ipiv) == min(m, n)
            lda  = max(1,stride(A, 2))
            info = Ref{BlasInt}()
            ccall((@blasfunc($getrf), libblas), Cvoid,
                    (Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                    Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                    m, n, A, lda, ipiv, info)
            LAPACK.chkargsok(info[])
            A, ipiv, info[] #Error code is stored in LU factorization type
        end
    end
end
function LinearAlgebra.lu!(A::StridedMatrix{T}, ipiv::AbstractVector{BlasInt}; check::Bool=true) where T
    lpt = LAPACK.getrf!(A, ipiv)
    check && LinearAlgebra.checknonsingular(lpt[3])
    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end