text:
using Vulkan
using VulkanExamples
using VulkanShaders
using ColorTypes
using FileIO
using Meshes

function main(output_png, points, texture; width = 1000, height = 1000, uv_coords = nothing)
    device, queue, messenger = init(enabled_features = PhysicalDeviceFeatures(:sampler_anisotropy))
    props = get_physical_device_properties(device.physical_device)

    uv_coords = something(uv_coords, [Point2f((1 .+ coordinates(p)) / 2) for p ∈ points])
    VertexType = PosUV{eltype(points),eltype(uv_coords)}
    format = FORMAT_R16G16B16A16_SFLOAT

    # define render pass
    target_attachment = AttachmentDescription(
        format,
        SAMPLE_COUNT_1_BIT,
        ATTACHMENT_LOAD_OP_CLEAR,
        ATTACHMENT_STORE_OP_STORE,
        ATTACHMENT_LOAD_OP_DONT_CARE,
        ATTACHMENT_STORE_OP_DONT_CARE,
        IMAGE_LAYOUT_UNDEFINED,
        IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL,
    )
    render_pass = RenderPass(
        device,
        [target_attachment],
        [
            SubpassDescription(
                PIPELINE_BIND_POINT_GRAPHICS,
                [],
                [AttachmentReference(0, IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL)],
                [],
            ),
        ],
        [
            SubpassDependency(
                SUBPASS_EXTERNAL,
                0;
                src_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_access_mask = ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
            ),
        ],
    )

    # create image and framebuffer
    fb_image = Image(
        device,
        IMAGE_TYPE_2D,
        format,
        Extent3D(width, height, 1),
        1,
        1,
        SAMPLE_COUNT_1_BIT,
        IMAGE_TILING_OPTIMAL,
        IMAGE_USAGE_COLOR_ATTACHMENT_BIT | IMAGE_USAGE_TRANSFER_SRC_BIT,
        SHARING_MODE_EXCLUSIVE,
        [0],
        IMAGE_LAYOUT_UNDEFINED,
    )
    fb_image_memory = DeviceMemory(fb_image, MEMORY_PROPERTY_DEVICE_LOCAL_BIT)
    fb_image_view = ImageView(
        fb_image.device,
        fb_image,
        IMAGE_VIEW_TYPE_2D,
        format,
        ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
        ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
    )
    framebuffer = Framebuffer(render_pass.device, render_pass, [fb_image_view], width, height, 1)

    # prepare vertex and index data
    p = PolyArea(Meshes.CircularVector(points))
    mesh = discretize(p, FIST())
    vdata = VertexType.(points, uv_coords)

    # prepare vertex buffer
    vbuffer = Buffer(device, buffer_size(vdata), BUFFER_USAGE_VERTEX_BUFFER_BIT, SHARING_MODE_EXCLUSIVE, [0])
    vmemory = DeviceMemory(vbuffer, vdata)

    # prepare shaders
    vert_shader = Shader(device, ShaderFile(joinpath(@__DIR__, "texture_2d.vert"), FormatGLSL()), [])
    frag_shader = Shader(
        device,
        ShaderFile(joinpath(@__DIR__, "texture_2d.frag"), FormatGLSL()),
        [DescriptorBinding(DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, 0, 1)],
    )
    shaders = [vert_shader, frag_shader]
    dset_layouts = create_descriptor_set_layouts(shaders)

    # build graphics pipeline
    shader_stage_cis = PipelineShaderStageCreateInfo.(shaders)
    vertex_input_state = PipelineVertexInputStateCreateInfo(
        [VertexInputBindingDescription(eltype(vdata), 0)],
        VertexInputAttributeDescription(eltype(vdata), 0),
    )
    input_assembly_state = PipelineInputAssemblyStateCreateInfo(PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP, false)
    viewport_state = PipelineViewportStateCreateInfo(
        viewports = [Viewport(0, 0, width, height, 0, 1)],
        scissors = [Rect2D(Offset2D(0, 0), Extent2D(width, height))],
    )
    rasterizer = PipelineRasterizationStateCreateInfo(
        false,
        false,
        POLYGON_MODE_FILL,
        FRONT_FACE_CLOCKWISE,
        false,
        0.0,
        0.0,
        0.0,
        1.0,
        cull_mode = CULL_MODE_BACK_BIT,
    )
    multisample_state = PipelineMultisampleStateCreateInfo(SAMPLE_COUNT_1_BIT, false, 1.0, false, false)
    color_blend_attachment = PipelineColorBlendAttachmentState(
        true,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD;
        color_write_mask = COLOR_COMPONENT_R_BIT | COLOR_COMPONENT_G_BIT | COLOR_COMPONENT_B_BIT,
    )
    color_blend_state = PipelineColorBlendStateCreateInfo(
        false,
        LOGIC_OP_CLEAR,
        [color_blend_attachment],
        Float32.((0.0, 0.0, 0.0, 0.0)),
    )
    pipeline_layout = PipelineLayout(device, dset_layouts, [])
    (pipeline, _...), _ = unwrap(
        create_graphics_pipelines(
            device,
            [
                GraphicsPipelineCreateInfo(
                    shader_stage_cis,
                    rasterizer,
                    pipeline_layout,
                    render_pass,
                    0,
                    0;
                    vertex_input_state,
                    multisample_state,
                    color_blend_state,
                    input_assembly_state,
                    viewport_state,
                ),
            ],
        ),
    )

    # load texture into buffer
    tdata = RGBA{Float16}.(load(texture))
    tbuffer = Buffer(
        device,
        buffer_size(tdata),
        BUFFER_USAGE_TRANSFER_DST_BIT | BUFFER_USAGE_TRANSFER_SRC_BIT,
        SHARING_MODE_EXCLUSIVE,
        [0],
    )
    tmemory = DeviceMemory(tbuffer, tdata)

    # create texture image
    timage = Image(
        device,
        IMAGE_TYPE_2D,
        format,
        Extent3D(size(tdata)..., 1),
        1,
        1,
        SAMPLE_COUNT_1_BIT,
        IMAGE_TILING_OPTIMAL,
        IMAGE_USAGE_TRANSFER_DST_BIT | IMAGE_USAGE_SAMPLED_BIT,
        SHARING_MODE_EXCLUSIVE,
        [0],
        IMAGE_LAYOUT_UNDEFINED,
    )
    timage_memory = DeviceMemory(timage, MEMORY_PROPERTY_DEVICE_LOCAL_BIT)

    command_pool = CommandPool(device, 0)

    # upload texture to image
    cbuffer, _... = unwrap(
        allocate_command_buffers(device, CommandBufferAllocateInfo(command_pool, COMMAND_BUFFER_LEVEL_PRIMARY, 1)),
    )
    begin_command_buffer(cbuffer, CommandBufferBeginInfo())
    cmd_pipeline_barrier(
        cbuffer,
        PIPELINE_STAGE_TOP_OF_PIPE_BIT,
        PIPELINE_STAGE_TRANSFER_BIT,
        [],
        [],
        [
            ImageMemoryBarrier(
                AccessFlag(0),
                ACCESS_TRANSFER_WRITE_BIT,
                IMAGE_LAYOUT_UNDEFINED,
                IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                QUEUE_FAMILY_IGNORED,
                QUEUE_FAMILY_IGNORED,
                timage,
                ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
            ),
        ],
    )
    cmd_copy_buffer_to_image(
        cbuffer,
        tbuffer,
        timage,
        IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
        [
            BufferImageCopy(
                0,
                size(tdata)...,
                ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1),
                Offset3D(0, 0, 0),
                Extent3D(size(tdata)..., 1),
            ),
        ],
    )
    cmd_pipeline_barrier(
        cbuffer,
        PIPELINE_STAGE_TRANSFER_BIT,
        PIPELINE_STAGE_FRAGMENT_SHADER_BIT,
        [],
        [],
        [
            ImageMemoryBarrier(
                ACCESS_TRANSFER_WRITE_BIT,
                AccessFlag(0),
                IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL,
                QUEUE_FAMILY_IGNORED,
                QUEUE_FAMILY_IGNORED,
                timage,
                ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
            ),
        ],
    )
    end_command_buffer(cbuffer)

    tuploaded = Fence(device)
    unwrap(queue_submit(queue, [SubmitInfo([], [], [cbuffer], [])]; fence = tuploaded))

    timage_view = ImageView(
        timage.device,
        timage,
        IMAGE_VIEW_TYPE_2D,
        format,
        ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
        ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
    )

    tsampler = Sampler(
        device,
        FILTER_LINEAR,
        FILTER_LINEAR,
        SAMPLER_MIPMAP_MODE_LINEAR,
        SAMPLER_ADDRESS_MODE_CLAMP_TO_BORDER,
        SAMPLER_ADDRESS_MODE_CLAMP_TO_BORDER,
        SAMPLER_ADDRESS_MODE_CLAMP_TO_BORDER,
        0,
        true,
        props.limits.max_sampler_anisotropy,
        false,
        COMPARE_OP_ALWAYS,
        0,
        0,
        BORDER_COLOR_FLOAT_OPAQUE_BLACK,
        false,
    )

    # create local image for transfer
    local_image = Image(
        device,
        IMAGE_TYPE_2D,
        format,
        Extent3D(width, height, 1),
        1,
        1,
        SAMPLE_COUNT_1_BIT,
        IMAGE_TILING_LINEAR,
        IMAGE_USAGE_TRANSFER_DST_BIT,
        SHARING_MODE_EXCLUSIVE,
        [0],
        IMAGE_LAYOUT_UNDEFINED,
    )

    local_image_memory = DeviceMemory(local_image, MEMORY_PROPERTY_HOST_COHERENT_BIT | MEMORY_PROPERTY_HOST_VISIBLE_BIT)

    command_buffer, _... = unwrap(
        allocate_command_buffers(device, CommandBufferAllocateInfo(command_pool, COMMAND_BUFFER_LEVEL_PRIMARY, 1)),
    )

    descriptor_pool = DescriptorPool(device, 1, [DescriptorPoolSize(DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, 1)])
    dsets = unwrap(allocate_descriptor_sets(device, DescriptorSetAllocateInfo(descriptor_pool, dset_layouts)))
    update_descriptor_sets(
        device,
        [
            WriteDescriptorSet(
                first(dsets),
                1,
                0,
                DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER,
                [DescriptorImageInfo(tsampler, timage_view, IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL)],
                [],
                [];
            ),
        ],
        [],
    )

    begin_command_buffer(command_buffer, CommandBufferBeginInfo())
    cmd_bind_vertex_buffers(command_buffer, [vbuffer], [0])
    cmd_bind_descriptor_sets(command_buffer, PIPELINE_BIND_POINT_GRAPHICS, pipeline_layout, 0, dsets, [])
    cmd_bind_pipeline(command_buffer, PIPELINE_BIND_POINT_GRAPHICS, pipeline)
    cmd_pipeline_barrier(
        command_buffer,
        PIPELINE_STAGE_TOP_OF_PIPE_BIT,
        PIPELINE_STAGE_TRANSFER_BIT,
        [],
        [],
        [
            ImageMemoryBarrier(
                AccessFlag(0),
                ACCESS_MEMORY_READ_BIT,
                IMAGE_LAYOUT_UNDEFINED,
                IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                QUEUE_FAMILY_IGNORED,
                QUEUE_FAMILY_IGNORED,
                local_image,
                ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
            ),
        ],
    )
    cmd_begin_render_pass(
        command_buffer,
        RenderPassBeginInfo(
            render_pass,
            framebuffer,
            Rect2D(Offset2D(0, 0), Extent2D(width, height)),
            [ClearValue(ClearColorValue((0.1f0, 0.1f0, 0.15f0, 1.0f0)))],
        ),
        SUBPASS_CONTENTS_INLINE,
    )
    cmd_draw(command_buffer, 4, 1, 0, 0)
    cmd_end_render_pass(command_buffer)
    cmd_pipeline_barrier(
        command_buffer,
        PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
        PIPELINE_STAGE_TRANSFER_BIT,
        [],
        [],
        [
            ImageMemoryBarrier(
                ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
                ACCESS_MEMORY_READ_BIT,
                IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL,
                IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL,
                QUEUE_FAMILY_IGNORED,
                QUEUE_FAMILY_IGNORED,
                fb_image,
                ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
            ),
        ],
    )
    cmd_copy_image(
        command_buffer,
        fb_image,
        IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL,
        local_image,
        IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
        [
            ImageCopy(
                ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1),
                Offset3D(0, 0, 0),
                ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1),
                Offset3D(0, 0, 0),
                Extent3D(width, height, 1),
            ),
        ],
    )
    end_command_buffer(command_buffer)

    GC.@preserve tbuffer tmemory timage timage_memory cbuffer unwrap(wait_for_fences(device, [tuploaded], true, 0))
    unwrap(queue_submit(queue, [SubmitInfo([], [], [command_buffer], [])]))
    GC.@preserve vbuffer vmemory timage timage_view timage_memory fb_image fb_image_view fb_image_memory local_image descriptor_pool command_pool command_buffer dsets unwrap(
        queue_wait_idle(queue),
    )

    # map image into a Julia array
    image = download_data(Vector{RGBA{Float16}}, local_image_memory, (width, height))

    save(output_png, image)

    image
end

main(
    joinpath(@__DIR__, "render.png"),
    Point2f[(-1.0, -1.0), (1.0, -1.0), (-1.0, 1.0), (1.0, 1.0)],
    joinpath(@__DIR__, "texture_2d.png"),
    width = 2048,
    height = 2048,
)


text:
using NIFTy

dom = RGDomain(500,500, harmonic = true)
pdom = PowerDomain(dom)

f1 = randn(dom)
f2 = randn(dom)
kernel1= Field(pdom, map(x -> 1/(1+x)^1.4, pdom._kvec))
kernel2 = Field(pdom, map(x -> 1/(1+x)^3, pdom._kvec))
full_field = MultiField(MultiField(f1, kernel1), MultiField(f2, kernel2))

distributor = DistributionOperator((dom, pdom), dom, pdom._pindices)
op1 = exp(hartley(distributor))
op2 = scale(target(op1), 80)(op1)
full_op = op1 + op2

using Plots
heatmap(full_op(full_field))

ext: jl


text:
function make_intake_info_compatible!(model::MetNet)
    for (intake, info) in model.intake_info
        if !(intake in model.rxns)
            delete!(model.intake_info, intake)
        end
    end
    return model
end


text:
using PyPlot, JUDI,LinearAlgebra, FFTW, DSP, SpecialFunctions, JOLI
using IterativeSolvers

using Shampoo

n = (201, 201)   # (x,y,z) or (x,z)
d = (5f0, 5f0)
o = (0f0, 0f0)

extentx = (n[1]-1)*d[1]
extentz = (n[2]-1)*d[2]

v0 = 2f0*ones(Float32,n)
v = deepcopy(v0)
v[101,101] = 2.1f0  # dirac perturbation

# Slowness squared [s^2/km^2]
m = convert(Array{Float32,2},(1f0 ./ v).^2)	# true model
m0 = convert(Array{Float32,2},(1f0 ./ v0).^2)
dm = vec(m-m0)

model = Model(n, d, o, m; nb = 200)
model0 = Model(n, d, o, m0; nb = 200)

# Model structure

dtS = 2f0
timeS = 2000f0
nt = Int64(timeS/dtS)+1

fmin = 10f0
fmax = 60f0
cfreqs = (fmin, fmin+5f0, fmax-5f0, fmax) # corner frequencies1
wavelet = reshape(cfreq_wavelet(500, nt, dtS/1f3, cfreqs; edge=hanning),:,1)
#wavelet = ricker_wavelet(timeS, dtS, 0.03f0)

################################### Source/receiver geometry ######################################

idx_wb = 30
Tm = judiTopmute(model0.n, idx_wb, 10)  # Mute water column
S = judiDepthScaling(model0)
Mr = S*Tm

# Set up receiver geometry

nrec = n[1]

xrec = range(0f0, stop=Float32((n[1]-1)*d[1]), length=nrec)
yrec = 0f0
zrec = range(14f0, stop=14f0, length=nrec)

# receiver sampling and recording time
timeR = timeS   # receiver recording time [ms]
dtR   = dtS
# Set up receiver structure
nsrc = 8

xsrc = convertToCell(range(0f0, stop=Float32((n[1]-1)*d[1]), length=nsrc))
# xsrc = convertToCell(range((n[1]-1)*d[1]/2f0, stop=(n[1]-1)*d[1]/2f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(9f0, stop=9f0, length=nsrc))

recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

opt = Options(free_surface=false,isic=false)
F0 = judiModeling(info, model0, srcGeometry, recGeometry; options=opt)
J = judiJacobian(F0, q)

PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure(); imshow(reshape(sqrt.(1f0./m), n[1], n[2])',extent=(0,extentx,extentz,0));
PyPlot.scatter(xrec,zrec,marker=".",color="red",label="receivers",s=1);
PyPlot.scatter(xsrc,zsrc,marker=".",color="white",label="sources",s=1);
legend(loc=3,fontsize=8);
xlabel("Lateral position [m]", fontsize=8);ylabel("Horizontal position [m]", fontsize=8);
fig = PyPlot.gcf();
title("Experimental setup")

# left preconditioners
P = FractionalIntegrationOp(nt,nsrc,nrec,0.5)

d_obs = J*dm

dm1 = J'*d_obs
dm2 = J'*P'*P*d_obs

figure();
subplot(2,1,1);
PyPlot.plot(wavelet);title("zoomed-in wavelet in time domain");
subplot(2,1,2);
PyPlot.plot(abs.(fftshift(fft(wavelet))));title("wavelet in frequency domain");

figure();
subplot(2,2,1)
imshow(dm1.data');title("J'*J*dm")
subplot(2,2,2)
imshow(dm2.data');title("J'*P'*P*J*dm")
subplot(2,2,3)
imshow(abs.(fftshift(fft(dm1.data')))/norm(abs.(fftshift(fft(dm1.data'))),Inf),cmap="jet",vmin=0,vmax=1)
subplot(2,2,4)
imshow(abs.(fftshift(fft(dm2.data')))/norm(abs.(fftshift(fft(dm2.data'))),Inf),cmap="jet",vmin=0,vmax=1)

maxit = 5

x1 = 0f0 .* dm
x1,his1 = lsqr!(x1,P*J*Mr,P*d_obs,atol=0f0,btol=0f0,conlim=0f0,maxiter=maxit,log=true,verbose=true)

x2 = 0f0 .* dm
x2,his2 = lsqr!(x2,J*Mr,d_obs,atol=0f0,btol=0f0,conlim=0f0,maxiter=maxit,log=true,verbose=true)

res1 = his1[:resnorm]
res2 = his2[:resnorm]

u1 = 20*log.(res1./norm(d_obs))
u2 = 20*log.(res2./norm(P*d_obs))

figure();
xlabel("iterations")
ylabel("ratio")
PyPlot.plot(res1/norm(d_obs),label="LSQR")
PyPlot.plot(res2/norm(P*d_obs),label="P-LSQR")
title("normalized least squares residual")
legend()

figure();
xlabel("iterations")
ylabel("dB")
PyPlot.plot(u1,label="LSQR")
PyPlot.plot(u2,label="P-LSQR")
title("normalized log-based least squares residual")
legend()

figure();
subplot(2,2,1)
imshow(reshape(Mr*x1,n)');title("LSQR")
subplot(2,2,2)
imshow(reshape(Mr*x2,n)');title("P-LSQR")
subplot(2,2,3)
imshow(abs.(fftshift(fft(reshape(Mr*x1,n)')))/norm(abs.(fftshift(fft(reshape(Mr*x1,n)'))),Inf),cmap="jet",vmin=0,vmax=1)
subplot(2,2,4)
imshow(abs.(fftshift(fft(reshape(Mr*x2,n)')))/norm(abs.(fftshift(fft(reshape(Mr*x2,n)'))),Inf),cmap="jet",vmin=0,vmax=1)


text:
"""
Constraint handler solving the semi-subproblem (fixed y, binary z).
"""
mutable struct SemiSubHdlr <: SCIP.AbstractConstraintHandler
    scip::SCIP.Optimizer
    solver
    finaltime::Float64 # seconds
    counter::Int
    primal_bound::Float64
    best_solution::Union{CandSol, Nothing}

    function SemiSubHdlr(scip, solver, finaltime)
        return new(scip, solver, finaltime, 0, Inf, nothing)
    end
end

"""
Constraint data, referencing the topology and binary edges variables.
"""
mutable struct SemiSubCons <: SCIP.AbstractConstraint{SemiSubHdlr}
    inst::Instance
    topo::Topology
    y::Array{JuMP.VariableRef}
    q::Array{JuMP.VariableRef}
end

function solve_semi_sub(ch::SemiSubHdlr, cons::SemiSubCons,
                        sol::Ptr{SCIP.SCIP_SOL}=C_NULL)
    ch.counter += 1

    # preparation
    inst = cons.inst
    topo = cons.topo
    ndiams = length(inst.diameters)
    narcs = length(topo.arcs)

    # fetch solution values
    ysol = SCIP.sol_values(ch.scip, cons.y, sol)
    qsol = SCIP.sol_values(ch.scip, cons.q, sol)
    @debug "  cand. sol:$(findall(!iszero, qsol))"

    # complete solution candidate
    zcand = fill(false, narcs, ndiams)
    zcand[ysol .> 0.5, 1] .= true  # pretend to use smallest diameter
    cand = CandSol(zcand, qsol, qsol.^2)

    # solve subproblem (from scratch)
    submodel, candarcs, z = make_semisub(inst, topo, cand, ch.solver)
    # TODO: solver.writemodels && writeLP(submodel, "sub_$(counter).lp", genericnames=false)
    settimelimit!(submodel, ch.solver, ch.finaltime - time())
    JuMP.optimize!(submodel)
    substatus = JuMP.termination_status(submodel)
    if substatus == MOI.OPTIMAL
        # have found improving solution?
        newobj = JuMP.objective_value(submodel)
        if newobj < ch.primal_bound
            ch.primal_bound = newobj
            znew = fill(false, narcs, ndiams)
            znew[candarcs, :] = JuMP.value.(z) .> 0.5
            ch.best_solution = CandSol(znew, qsol, qsol.^2)
            @debug "  new sol: $(ch.primal_bound)"
            return true # found improvement
        end
    elseif substatus != MOI.INFEASIBLE
        error("Unexpected status: $(:substatus)")
    end

    # TODO: what can we do about dual/primal stopping criterion?
    # TODO: add cut on objective var?

    return false # did not find improvement
end

function SCIP.check(ch::SemiSubHdlr,
                    constraints::Array{Ptr{SCIP.SCIP_CONS}},
                    sol::Ptr{SCIP.SCIP_SOL},
                    checkintegrality::Bool,
                    checklprows::Bool,
                    printreason::Bool,
                    completely::Bool)
    @assert length(constraints) == 1
    cons::SemiSubCons = SCIP.user_constraint(constraints[1])

    # Always solve semi-subproblem (to find improving solutions)
    # TODO: cache candidates to avoid re-solving?
    improvement = solve_semi_sub(ch, cons, sol)

    # Never accept the solution candidate.
    return SCIP.SCIP_INFEASIBLE
end

function enforce_semi_sub(ch::SemiSubHdlr, cons::SemiSubCons)
    # Always solve semi-subproblem (to find improving solutions)
    # TODO: cache candidates to avoid re-solving?
    improvement = solve_semi_sub(ch, cons)

    # Always add no-good to cut off this candidate.
    ysol = SCIP.sol_values(ch.scip, cons.y)
    active = ysol .> 0.5
    coefs = 2.0 * active .- 1.0
    ci = @cb_constraint(ch.scip, coefs ⋅ cons.y ≤ sum(active) - 1)

    return SCIP.SCIP_CONSADDED
end

function SCIP.enforce_lp_sol(ch::SemiSubHdlr,
                             constraints::Array{Ptr{SCIP.SCIP_CONS}},
                             nusefulconss::Cint,
                             solinfeasible::Bool)
    @assert length(constraints) == 1
    cons::SemiSubCons = SCIP.user_constraint(constraints[1])
    return enforce_semi_sub(ch, cons)
end

function SCIP.enforce_pseudo_sol(ch::SemiSubHdlr,
                                 constraints::Array{Ptr{SCIP.SCIP_CONS}},
                                 nusefulconss::Cint,
                                 solinfeasible::Bool,
                                 objinfeasible::Bool)
    @assert length(constraints) == 1
    cons::SemiSubCons = SCIP.user_constraint(constraints[1])
    return enforce_semi_sub(ch, cons)
end

function SCIP.lock(ch::SemiSubHdlr,
                   constraint::Ptr{SCIP.SCIP_CONS},
                   locktype::SCIP.SCIP_LOCKTYPE,
                   nlockspos::Cint,
                   nlocksneg::Cint)
    cons::SemiSubCons = SCIP.user_constraint(constraint)
    for y in cons.y
        # TODO: understand why lock is called during SCIPfree, after the
        # constraint should have been deleted already. Does it mean we should
        # implement CONSTRANS?
        var_::Ptr{SCIP.SCIP_VAR} = SCIP.var(ch.scip, JuMP.index(y))
        var_ != C_NULL || continue  # avoid segfault!

        SCIP.@SC SCIP.SCIPaddVarLocksType(
            ch.scip, var_, locktype, nlockspos + nlocksneg, nlockspos + nlocksneg)
    end
end


text:
module Qubits

using ..Clifford
using ..Waveforms
using InstrumentControl

export Qubit
export cosInit
export gaussInit
export generalInit

# From clifford
export benchmark1Qubit
export benchmark2Qubit
export Pulse
export gateNames

#========================== The Qubit DataType ================================#
# For the purposes of this code, the only info needed about a qubit is its
# resonant frequency relative to the local oscillator, which boards/lines
# are able to communicate with it, and the on-board data for gates.

# Depending on if lines permit, we can control our qubits with only XY control
# or with full XYZ.  If the Z line is absent, the underlying type is a QubitNoZ,
# and if present it is a QubitWithZ.  The user only creates "Qubit" objects,
# and depending on whether a Z line is specified, the program generates the
# proper concrete type.

abstract Qubit

type QubitWithZ <: Qubit
  IFreq::Float64
  lineXYI::Tuple{Instrument,Int} # (board, channel)
  lineXYQ::Tuple{Instrument,Int} #  "
  lineZ::Tuple{Instrument,Int}   #  "
  pulseConvert::Matrix        # 10x3 matrix, which maps the 10 basic pulses
                              # to how the DAC's memory refers to them

  waveforms::Dict{Pulse, ExactWaveform} # See below.
end

type QubitNoZ <: Qubit
  IFreq::Float64
  lineXYI::Tuple{Instrument,Int} # (board, channel)
  lineXYQ::Tuple{Instrument,Int} #  "
  pulseConvert::Matrix        # 7x2 mapping basic pulses to XYI and XYQ labels
                              # in the DAC's memory (for AWG, this is strings,
                              # for reasons that only frustrate me..)

  waveforms::Dict{Pulse, ExactWaveform} # See below.
end

# Accessing a qubit as a dictionary will alter the contained dictionary of
# pulses.  Though the dictionary itself only contains ExactWaveform values,
# the methods below provide ways to convert FloatWaveform and Vector{Float64}
# values to ExactWaveform.  An explicit convert method cannot be written, since
# the ExactWaveform produced depends on the Pulse and on the IFreq of the qubit.
import Base.setindex!
function setindex!(q::Qubit, w::Vector{Float64}, p::Pulse)
  if p[1] == 7
    q.waveforms[p] = ExactWaveform(UInt16[], UInt16[], UInt16[], true)
  elseif length(w) != floatIdleLength
      error("FloatWaveform pulses must contain exactly "*
                                 string(floatIdleLength)*" points")
  elseif p[1] < 7
    (xyi, xyq) = IQgen(q.IFreq, p, w)
    q.waveforms[p] = ExactWaveform(xyi, xyq, UInt16[], true)
  elseif p[1] < 11 && isa(q, QubitWithZ)
    q.waveforms[p] = ExactWaveform(UInt16[], UInt16[],
                    map(x -> UInt16(offsetValue + round(x)), w), true)
  else # Trying to set Z gates on a QubitNoZ object
    println("Warning: no gates set by this operation.")
  end
end

# A helper method for the above that uses the IFreq and the intended phase to
# determine the I and Q pulses.
function IQgen(IFreq::Float64, pulse::Pulse, window::Vector{Float64})
  if pulse[1] > 6
    error("Not a pulse in need of IQ mixing")
  else
    phase = im ^ pulse[1]
    helix = phase * [exp(im * 2π * IFreq * p / sampleRate) * window[p]
                                  for p in 1:length(window)]
    waveformI = map(z->UInt16(offsetValue + round(real(z))), helix)
    waveformQ = map(z->UInt16(offsetValue + round(imag(z))), helix)
    waveformI, waveformQ
  end
end

# For FloatWaveforms, use the contained wavedata field.
function setindex!(q::Qubit, w::FloatWaveform, p::Pulse)
  setindex!(q, w.wavedata, p)
end

# For ExactWaveforms, use them as they are.
function setindex!(q::Qubit, w::ExactWaveform, p::Pulse)
  if p[1] < 7 || isa(QubitWithZ)
    w.dirty = true
    q.waveforms[p] = w
  else # Trying to set Z gates on a QubitNoZ object
    println("Warning: no gates set by this operation.")
  end
end

# Always display the ExactWaveform for the user to see.
import Base.getindex
function getindex(q::Qubit, p::Pulse)
  return q.waveforms[p]
end

# ================= Initializing a Qubit's information ========================#

# A user shouldn't have to specify the matrix or the dictionary to define a
# qubit.  All that they should have to input is IFreq and the line info.
# A user can specify "False" for the Z line to generate a QubitNoZ object,
# or omit it entirely.
function Qubit(IFreq::Float64, lineXYI::Tuple{Int,Int}, lineXYQ::Tuple{Int,Int},
               lineZ=false)
  if isa(lineZ, Bool) && !lineZ
    ret = QubitNoZ(IFreq, lineXYI, lineXYQ, fill(-1, (7,2)), Dict())
  elseif isa(lineZ, Tuple)
    ret = QubitWithZ(IFreq, lineXYI, lineXYQ, lineZ, fill(-1, (10,3)), Dict())
  end
  println("To initialize pulse shapes for this Qubit, please run one of the "*
           "init routines:\n\tgaussInit\n\tcosInit\n\tgeneralInit")
  ret
end

# Initializing the dictionary shouldn't require manually inputting up to 10
# different pulse shapes, so we provide quick ways to make gaussian-envelope and
# cos-envelope pulses.  These are all FloatWaveform objects with all the
# assumptions mentioned above.  We also provide a generalInit function that
# takes a shape for the window, though it must be floatIdleLength points exactly
# or an error will be thrown down the line.

# These functions assume that half-amplitude pulses provide half the rotation.
# This is a good starting point but is not accurate for high-fidelity gates.

# The XY and Z arguments specify whether the pulses created are for the XY
# control or the Z control.  By default, both are on, but it is likely the
# ideal amplitude for each is different.  For QubitNoZ objects, setting Z to
# true or false makes no difference.
function gaussInit(q::Qubit, amplitude, sigma, XY::Bool = true, Z::Bool = true)
  pulseShape = [amplitude*exp(-(x - (floatIdleLength + 1)/2)^2 / (2*sigma^2))
                                for x in 1:floatIdleLength]
  pulseShape -= pulseShape[1] # To eliminate the tails
  generalInit(q, pulseShape, XY, Z)
end

function cosInit(q::Qubit, amplitude, XY::Bool = true, Z::Bool = true)
  pulseShape = [amplitude*(1-cos(x/floatIdleLength)/2)
                                for x in 1:floatIdleLength]
  generalInit(q, pulseShape, XY, Z)
end

function generalInit(q::Qubit, pulseShape, XY::Bool = true, Z::Bool = true)
  if XY
    q[Xpi] = q[Ypi] = pulseShape
    q[Xpi2] = q[Ypi2] = q[X3pi2] = q[Y3pi2] = pulseShape/2
  end
  if Z && isa(q, QubitWithZ)
    q[Zpi] = pulseShape
    q[Zpi2] = pulseShape/2
    q[Z3pi2] = -pulseShape/2
  end
  if q.pulseConvert[7] == -1
    q[I] = [0.0]
  end
end

end # End module Qubit


text:
"""
    Double64s is a package

exported types: [`Double64`](@ref), [`ComplexD64`](@ref)

exported field accessors: [`Hi`](@ref), [`Lo`](@ref), [`HiLo`](@ref)
"""
module Double64s

export Double64, ComplexD64,
  Hi, Lo, HiLo,
  NaND64, InfD64,
  signs, signbits,
  fastabs, fastabs2

using ErrorfreeArithmetic, Quadmath

# online help text for the types

"""
   fastabs(x)

abs(Hi(x))
""" fastabs

"""
   fastabs2(x)

abs2(Hi(x))
""" fastabs2

include("types/double64.jl")
include("types/constructors.jl")
include("support/fast_minmax.jl")
include("types/comparison.jl")
include("types/promotion.jl")
include("types/predicates.jl")

include("math/prearith.jl")
include("math/arith.jl")
include("math/rounding.jl")
include("math/fma.jl")
include("math/roots.jl")

for F in (:abs, :acos, :acosh, :asin, :asinh, :atan, :atanh, :ceil, :cos, :cosh,
          :exp, :exp2, :exp10, :expm1, :log, :log2, :log10, :log1p,
          :sin, :sinh, :sqrt)
    @eval Base.$F(x::Double64) = Double64($F(Float128(x)))
end

end  # Double64s


text:
# Calculator.jl
module Calculator

export interest, rate

"""
 interest(amount, rate)

Calculate interest from an `amount` and interest rate of `rate`.
"""
function interest(amount, rate)
 return amount * (1 + rate)
end

"""
 rate(amount, interest)

Calculate interest rate based on an `amount` and `interest`.
"""
function rate(amount, interest)
 return interest / amount
end

end # module


text:
# test FluxModelsPose2Pose2

using Test
using LinearAlgebra
using Flux
using RoME

@testset "Basic funtionality tests on FluxModelsPose2Pose2" begin

##

mdls = [RoME.buildPose2OdoNN_01_FromElements( randn(4,8),
                                              randn(8),
                                              randn(8,48),
                                              randn(8),
                                              randn(2,8),
                                              randn(2)) for i in 1:10];

##

# start with a basic factor graph

mvnNaive = MvNormal(zeros(3), diagm([1.0;1.0;0.01]))

jvd = zeros(25,4)
pp = FluxModelsPose2Pose2(mdls, jvd, mvnNaive, 0.5)

##

fg = generateCanonicalFG_ZeroPose2()
addVariable!(fg, :x1, Pose2)

addFactor!(fg, [:x0;:x1], pp)

##

pts = approxConv(fg, :x0x1f1, :x1)

#

end


text:
export StaticString, ODEString, AccelString, StaticNEB, ODENEB, AccelNEB

neb_string_shared_docs =  """
### Shared Parameters
* `tol` : residual tolerance
* `maxnit` : maximum number of iterations
* `precon_scheme` : preconditioner scheme (localPrecon(), globalPrecon())
* `path_traverse` : how to travserse path to compute energy (serial(), palindrome())
* `fixed_ends` : true/false for fixing the ends of the path, default is false
* `verbose` : how much information to print (0: none, 1:end of iteration, 2:each iteration, 3:file log)
"""

@def neb_string_params begin
   tol::Float64 = 1e-5
   maxnit::Int = 1000
   maxtol::Float64 = 1e3
   precon_scheme = localPrecon()
   path_traverse = serial()
   fixed_ends = false
   verbose::Int = 2
end

neb_only_docs = """
### NEB only Parameters:
* `k` : spring constant
* `interp` : degree of interpolant (1: central differences, >1: splines of degree interp)
"""
@def neb_params begin
   k::Float64 = 0.1
   interp = 1
end

"""
`StaticString`: the most basic string variant, integrating potential gradient
flow with a fixed step-size.

### Parameters:
* `alpha` : Euler method integration step length

$(neb_string_shared_docs)
"""
@with_kw type StaticString
   alpha::Float64
   # ------ shared parameters ------
   @neb_string_params
end

"""
`ODEString`: string method with adaptive ODE step size selection.

### Parameters:
* `reltol` : ode solver relative tolerance
* `threshold` : threshold for error estimate
* `a0` : initial step, if a0 is not passed then use a default

$(neb_string_shared_docs)
"""
@with_kw type ODEString
   reltol::Float64  # please think about reducing it to one tol parameter `odetol`
   threshold::Float64 = 1.0   # threshold for error estimate; we want to get rid of this
   a0 = nothing      # if a0 is not passed then use a default
   # ------ shared parameters ------
   @neb_string_params
end

"""
`AccelString`: string method using momentum descent to accelerate energy minimisation.

### Parameters:
* `a0` : initial step, if a0 is not passed then use a default
* `b` : momentum term damping coefficient
* `finite_diff_scheme` : choice of discretisation of dumped wave equation

$(neb_string_shared_docs)
"""
@with_kw type AccelString
   h = nothing      # if h is not passed then use a default
   a0 = nothing      # if a0 is not passed then use a default
   b = nothing      # if b is not passed then optimal value is used
   reltol::Float64  # please think about reducing it to one tol parameter `odetol`
   finite_diff_scheme = :central
   redistrib = :canonical
   # ------ shared parameters ------
   @neb_string_params
end

"""
`StaticNEB`: the most basic NEB variant, integrating potential gradient
flow with a fixed step-size.

### Parameters:
* `alpha` : Euler method integration step length

$(neb_only_docs)

$(neb_string_shared_docs)
"""
@with_kw type StaticNEB
   alpha::Float64
   @neb_params
   # ------ shared parameters ------
   @neb_string_params
end

"""
`ODEString`: NEB method with adaptive ODE step size selection.

### Parameters:
* `reltol` : ode solver relative tolerance
* `threshold` : threshold for error estimate
* `a0` : initial step, if a0 is not passed then use a default

$(neb_only_docs)

$(neb_string_shared_docs)
"""
@with_kw type ODENEB
   reltol::Float64  # please think about reducing it to one tol parameter `odetol`
   threshold::Float64 = 1.0   # threshold for error estimate; we want to get rid of this
   a0 = nothing      # if a0 is not passed then use a default
   @neb_params
   # ------ shared parameters ------
   @neb_string_params
end

"""
`AccelNEB`: NEB method, using momentum descent to accelerate energy minimisation.

### Parameters:
* `a0` : initial step, if a0 is not passed then use a default
* `b` : momentum term damping coefficient
* `finite_diff_scheme` : choice of discretisation of dumped wave equation

$(neb_only_docs)

$(neb_string_shared_docs)
"""
@with_kw type AccelNEB
    h = nothing      # if h is not passed then use a default
    a0 = nothing      # if a0 is not passed then use a default
    b = nothing      # if b is not passed then optimal value is used
    reltol::Float64  # please think about reducing it to one tol parameter `odetol`
    finite_diff_scheme = :central
    redistrib = :canonical
    @neb_params
    # ------ shared parameters ------
    @neb_string_params
end

function String(step, args...; kwargs...)
   if step == :static
      return StaticString(args...; kwargs...)
   elseif step == :ode
      return ODEString(args...; kwargs...)
   elseif step == :accel
       return AccelString(args...; kwargs...)
   else
      error("`String`: unknown step selection mechanism $(step)")
   end
end

function NEB(step, args...; kwargs...)
   if step == :static
      return StaticNEB(args...; kwargs...)
   elseif step == :ode
      return ODENEB(args...; kwargs...)
   elseif step == :accel
       return AccelNEB(args...; kwargs...)
   else
      error("`NEB`: unknown step selection mechanism $(step)")
   end
end

solver(method::StaticString) = Euler(h=method.alpha)
solver(method::ODEString) = ODE12r(rtol=method.reltol, threshold=method.threshold, h=method.a0)
# solver(method::AccelString) = momentum_descent(h=method.a0, b=method.b, finite_diff=method.finite_diff_scheme)
solver(method::AccelString) = momentum_descent(h=method.h, h0=method.a0, rtol=method.reltol, b=method.b,
fd_scheme=method.finite_diff_scheme, redistrib=method.redistrib)
solver(method::StaticNEB) = Euler(h=method.alpha)
solver(method::ODENEB) = ODE12r(rtol=method.reltol, threshold=method.threshold, h=method.a0)
# solver(method::AccelNEB) = momentum_descent(h=method.a0, b=method.b, finite_diff=method.finite_diff_scheme)
solver(method::AccelNEB) = momentum_descent(h=method.a0, b=method.b, rtol=method.reltol,
fd_scheme=method.finite_diff_scheme, redistrib=method.redistrib)


text:
# state names ad abbreviations

states = ["Alabama","AL","Alaska","AK","Arizona","AZ","Arkansas","AR","California","CA","Colorado","CO","Connecticut","CT","Delaware","DE","District of Columbia","DC","Florida","FL","Georgia","GA","Hawaii","HI","Idaho","ID","Illinois","IL","Indiana","IN","Iowa","IA","Kansas","KS","Kentucky","KY","Louisiana","LA","Maine","ME","Maryland","MD","Massachusetts","MA","Michigan","MI","Minnesota","MN","Mississippi","MS","Missouri","MO","Montana","MT","Nebraska","NE","Nevada","NV","New Hampshire","NH","New Jersey","NJ","New Mexico","NM","New York","NY","North Carolina","NC","North Dakota","ND","Ohio","OH","Oklahoma","OK","Oregon","OR","Pennsylvania","PA","Rhode Island","RI","South Carolina","SC","South Dakota","SD","Tennessee","TN","Texas","TX","Utah","UT","Vermont","VT","Virginia","VA","Washington","WA","West Virginia","WV","Wisconsin","WI","Wyoming","WY"]

println(states)


text:
# #############################################################################
# #############################################################################
#Cavity analysis functions and global variables

#This file contains the main driving functions for evaluating a cavity, as well as the state variables and setting/getting functions for the simulation and cavity parameters. Typically, a control file will run on a multiprocess julia instance

#   shell-prompt> julia -p <num of parallel procs> <ctrlfile>

#The file starts with the declaration

#   @everywhere using CavChaos

#to initiate the module's methods into the namespace on every process. The following functions can then be used in a control file to set up a simulation run:

#>> Cavity properties

#   set_cavity_bnd!(bnd::Boundary)
#   Sets the desired cavity boundary to <bnd>.

#   bnd = get_cavity_bnd()
#   Returns the current cavity boundary.

#   set_cavity_idx!(idx::RefractiveIndex)
#   Sets the desired cavity index distribution to <idx>.

#   idx = get_cavity_idx()
#   Returns the current cavity index distribution.

#>> Solver parameters

#   set_solver_params!(;kwargs...)
#   Sets the solver parameter represented by keywords to the numerical value associated to it in <kwargs>. Appropriate keywords are: maxpathlength, maxbounces, relativetolerance, absolutetolerance.

#   value = get_solver_param(key::Union(AbstractString,Symbol))
#   Returns the current value of the solver parameter represented by the symbol or string <key>.

#>> Initial conditions

#   set_init_params!(;kwargs...)
#   Sets the initial conditions autogeneration parameter represented by keywords to the numerical value associated to it in <kwargs>. Appropriate keywords are: symmetry, thetares, chires.

#   value = get_init_param(key::Union(AbstractString,Symbol))
#   Returns the current value of the initial conditions autogeneration parameter represented by the symbol or string <key>.

#   set_pssinitarray!(A::Array{Float64,2})
#   Sets the custom initial conditions as an array <A> of initial (theta,chi) values as its rows. If nonempty, the autogeneration parameters will be ignored, and the initial conditions fed to the solver will be computed from <A> directly. This is useful for choosing points on the Poincare Surface of Section for performing simulation runs.

#   pssinitarray = get_pssinitarray()
#   Returns the current custom/autogenerated (theta,chi)-based initial conditions array <pssinitarray>.

#   set_initarray!(A::Array{Float64,2})
#   Sets the custom initial conditions as an array <A> of initial (r,theta,phi) values as its rows. If nonempty, the autogeneration parameters and the (theta,chi)-based custom initial conditions array will be ignored, and the initial conditions fed to the solver will be taken from <A> directly. This is useful for fine control over the initial ray position and direction when performing simulation runs.

#   initarray = get_initarray()
#   Returns the current custom/autogenerated (r,theta,phi)-based initial conditions array <initarray>.

#>> Results directories

#   set_resultsroot!(dir::AbstractString)
#   Sets the results root directory (in which a "results" subdirectory stores all output data) to <dir>.

#   get_resultsroot()
#   Returns the results root directory.

#   get_resultsdir()
#   Returns the configuration-specific results directory

#>> Results to record

#   set_results_params!(;kwargs...)
#   Sets the result represented by keywords to a boolean value associated to it in <kwargs>, indicating whether that result will be recorded. Appropriate keywords are: pssinitarray, initarray, initialconditions, raypath, bounceindices, bouncepoints, cavityimage, rayimage, pathlengths, actions, modeproperties, farfield.

#   set_results_params!(args...)
#   unset_results_params!(args...)
#   For a convenient syntax, if the above results keywords are passed as regular symbol or string arguments to set_results!, then they will be recorded. If passed to unset_results!, they will not.

#   boolean = get_results_param(key::Union(AbstractString,Symbol))
#   Returns whether the result represented by the keyword represented by the string or symbol <key> will be recorded.

#Finally, the actual run is started by a call to the function run_cavity().

# #############################################################################
# #############################################################################
#Initiate

#Dependencies
#require("util.jl")
#require("boundary.jl")
#require("refractive_index.jl")
#require("solver.jl")
#require("io.jl")
#require("ray_plot.jl")
#require("ray_analysis.jl")
#require("ray_combined.jl")

# #############################################################################
# #############################################################################
#Settings as global state

#Cavity objects
cavity_bnd = FourierBnd()
cavity_idx = UniformIdx()

#Solver parameters
solver = Dict{Symbol,Number}()
solver[:maxpathlength] = 200.0
solver[:maxbounces] = 500
solver[:relativetolerance] = 1.0e-12
solver[:absolutetolerance] = 1.0e-12

#Initial conditions array
init = Dict{Symbol,Number}()
init[:symmetry] = 1
init[:thetares] = 25
init[:chires]= 25
pssinitarray = Array(Float64,0,2)
initarray = Array(Float64,0,3)

#Root directory to store results
resultsroot = "."

#Results to record
results = Dict{Symbol,Bool}()
results[:pssinitarray] = false
results[:initarray] = true
results[:initialconditions] = false
results[:raypath] = false
results[:bounceindices] = false
results[:bouncepoints] = true
results[:cavityimage] = true
results[:rayimage] = true
results[:pathlengths] = true
results[:actions] = true
results[:modeproperties] = false
results[:farfield] = false

# #############################################################################
# #############################################################################
#Setters and unsetters and getters

#Cavity parameters
function set_cavity_bnd!(bnd::Boundary)
    global cavity_bnd = bnd
    return nothing
end

function get_cavity_bnd()
    return cavity_bnd::Boundary
end

function set_cavity_idx!(idx::RefractiveIndex)
    global cavity_idx = idx
    return nothing
end

function get_cavity_idx()
    return cavity_idx::RefractiveIndex
end

#Solver parameters
function set_solver_params!(;kwargs...)
    for i=1:length(kwargs)
        if haskey(solver,kwargs[i][1])
            solver[kwargs[i][1]] = kwargs[i][2]
        else
            error("No solver parameter $(kwargs[i][1]).")
        end
    end
    return nothing
end

get_solver_param(key::Symbol) = solver[key]::Number
get_solver_param(key::AbstractString) = solver[symbol(key)]::Number

#Initial conditions parameters
function set_init_params!(;kwargs...)
    for i=1:length(kwargs)
        if haskey(init,kwargs[i][1])
            init[kwargs[i][1]] = kwargs[i][2]
        else
            error("No initial conditions parameter $(kwargs[i][1]).")
        end
    end
    return nothing
end

get_init_param(key::Symbol) = init[key]::Number
get_init_param(key::AbstractString) = init[symbol(key)]::Number

function set_pssinitarray!(A::Array{Float64,2})
    if size(A,2) != 2
        error("PSS initial conditions array must have (theta,chi) pairs as rows.")
    end
    global pssinitarray = A
    return nothing
end

function get_pssinitarray()
    return pssinitarray::Array{Float64,2}
end

function set_initarray!(A::Array{Float64,2})
    if size(A,2) != 3
        error("Initial conditions array must have (r,theta,phi) values as rows.")
    end
    global initarray = A
    return nothing
end

function get_initarray()
    return initarray::Array{Float64,2}
end

#Root directory to store results
function set_resultsroot!(directory::AbstractString)
    global resultsroot = directory
    return nothing
end

function get_resultsroot()
    return resultsroot::AbstractString
end

#Results directory for the current configuration
function get_resultsdir()
    #Autogenerate initarray
    initarray = get_initarray()
    if isempty(initarray)
        pssinitarray = get_pssinitarray()
        if isempty(pssinitarray)
            pssinitarray = gen_pssinitarray(get_init_param(:thetares),
                                            get_init_param(:chires),
                                            get_init_param(:symmetry))
        end
        initarray = pss2initarray(get_cavity_bnd(),pssinitarray)
    end

    #Get run parameter hash
    run_params_hash::UInt64 =
        get_runparamshash(  initarray,
                            get_solver_param(:maxpathlength),
                            get_solver_param(:maxbounces),
                            get_solver_param(:relativetolerance),
                            get_solver_param(:absolutetolerance),
                            get_results_param(:pssinitarray),get_results_param(:initarray),
                            get_results_param(:initialconditions),get_results_param(:raypath),
                            get_results_param(:bounceindices),get_results_param(:bouncepoints),
                            get_results_param(:cavityimage),get_results_param(:rayimage),
                            get_results_param(:pathlengths),get_results_param(:actions),
                            get_results_param(:modeproperties),get_results_param(:farfield))

    return getresultsdir(   run_params_hash,get_cavity_bnd(),get_cavity_idx(),
                            get_resultsroot(),makedir=false)
end

#Results parameters
function set_results_params!(;kwargs...)
    for i=1:length(kwargs)
        if haskey(results,kwargs[i][1])
            results[kwargs[i][1]] = kwargs[i][2]
        else
            error("No result named $(kwargs[i][1]).")
        end
    end
    return nothing
end

function set_results_params!(args...)
    set_results!(;[(symbol(arg),true) for arg in args]...)
end

function unset_results_params!(args...)
    set_results!(;[(symbol(arg),false) for arg in args]...)
end

get_results_param(key::Symbol) = results[key]::Bool
get_results_param(key::AbstractString) = results[symbol(key)]::Bool

# #############################################################################
# #############################################################################
#Main driver

#Profile a single cavity by running multiple rays in it
function run_cavity()

    run_rays(   #Cavity parameters
                get_cavity_bnd(),get_cavity_idx();

                #Initial conditions
                pssinitarray = get_pssinitarray(),
                initarray = get_initarray(),

                #Initial conditions generator parameters
                symmetry = get_init_param(:symmetry),
                thetares = get_init_param(:thetares),
                chires= get_init_param(:chires),

                #Solver parameters
                maxpathlength = get_solver_param(:maxpathlength),
                maxbounces = get_solver_param(:maxbounces),
                relativetolerance = get_solver_param(:relativetolerance),
                absolutetolerance = get_solver_param(:absolutetolerance),

                #Results root directory
                resultsroot = get_resultsroot(),

                #Results to record
                record_pssinitarray = get_results_param(:pssinitarray),
                record_initarray = get_results_param(:initarray),
                record_initialconditions = get_results_param(:initialconditions),
                record_raypath = get_results_param(:raypath),
                record_bounceindices = get_results_param(:bounceindices),
                record_bouncepoints = get_results_param(:bouncepoints),
                record_cavityimage = get_results_param(:cavityimage),
                record_rayimage = get_results_param(:rayimage),
                record_pathlengths = get_results_param(:pathlengths),
                record_actions = get_results_param(:actions),
                record_modeproperties = get_results_param(:modeproperties),
                record_farfield = get_results_param(:farfield),
            )

end


text:
using PlotGraphviz
using Test

@testset "PlotGraphviz.jl" begin
    # Write your tests here.
end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
# Include packages
using Polynomials: polyfit, polyval
using DelimitedFiles: readdlm

####################
## main
####################

# define the filepath & filename
fdir = "./data"
fname = "windspeed.dat"

# load ascii file
dataorg = readdlm(joinpath(fdir,fname))
years = dataorg[:,1]
wind = dataorg[:,2]
stdval = dataorg[:,3]
# regression analysis
lin_p = polyfit(years, wind, 1)
quad_p = polyfit(years, wind, 2)
years_ext = years[end]:2100.0

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

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()
#plotlyjs()

# original data
plot(years, wind, linetype=:scatter, ms=10., tickfont=12, lab="Raw data",
     legendfont=font(14,"sans-serif"), legend=:bottomleft,
     title="Regression analysis", titlefont=font("sans-serif",16),
     xlabel="Year", ylabel="Wind speed (m/s)",guidefont=font("sans-serif",16),
     size=(800,600),
     )
# linear
plot!(years, polyval(lin_p,years), line=(:solid,2), linecolor=:red, lab="Linear")
plot!(years_ext, polyval(lin_p,years_ext), line=(:dash,2), linecolor=:red, lab="Linear to 2100")
# quad
plot!(years, polyval(quad_p,years), line=(:solid,2), linecolor=:green, lab="Quad")
plot!(years_ext, polyval(quad_p,years_ext), line=(:dash,2), linecolor=:green, lab="Quad to 2100")
xlims!(1945, 2105)

# save figure
savefig(joinpath(figdir,"regression.png"))
####################


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir)

df = DataFrame!(CSV.File(scm_path("..", "data", "WaffleDivorce.csv"), delim=';'));
df = DataFrame(
  :s => df[:, :South],
  :a => df[:, :MedianAgeMarriage],
  :m => df[:, :Marriage],
  :w => df[:, :WaffleHouses],
  :d => df[:, :Divorce]
);

fname = scm_path("..", "examples", "SR", "SR6.4.3", "sr6.4.3.dot")
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

d = OrderedDict(
  :s => [:a, :m, :w],
  :a => [:m, :d],
  :m => [:d],
  :w => [:d]
);

dag = DAG("sr6.4.3", d, df=df);

@testset "sr6_4_3b" begin

  adjustmentsets = adjustment_sets(dag, :w, :d)
  @test adjustmentsets == Array{Symbol,1}[]

  c = [:s]
  ag = ancestral_graph(dag; c=c)
  @test names(ag, 1) == [:a, :m, :w, :d]
  @test Array(ag) == [
    0  10  10  10;
   10   0  10  10;
   10  10   0  10;
   10  10  10   0
  ]

end


text:
#######################################################################################
# Simplest Loop method
# Makes the simplest possible method for teaching and performance testing
#######################################################################################
struct LoopRK4 <: DiffEqBase.AbstractODEAlgorithm end
export LoopRK4

# Out-of-place
# No caching, good for static arrays, bad for arrays
@muladd function DiffEqBase.__solve(prob::ODEProblem{uType,tType,false},
                          alg::LoopRK4;
                          dt = error("dt is required for this algorithm"),
                          save_everystep = true,
                          save_start = true,
                          adaptive = false,
                          dense = false,
                          save_end = true,
                          kwargs...) where {uType,tType}
  @assert !adaptive
  @assert !dense
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  t = tspan[1]
  tf = prob.tspan[2]
  ts = tspan[1]:dt:tspan[2]

  if save_everystep && save_start
      us = Vector{typeof(u0)}(undef,length(ts))
      us[1] = u0
  elseif save_everystep
      us = Vector{typeof(u0)}(undef,length(ts)-1)
  elseif save_start
      us = Vector{typeof(u0)}(undef,2)
      us[1] = u0
  else
      us = Vector{typeof(u0)}(undef,1) # for interface compatibility
  end

  u = u0
  half = convert(eltype(u0),1//2)
  sixth = convert(eltype(u0),1//6)

  for i in 2:length(ts)
      uprev = u; t = ts[i]
      k1 = f(u,p,t)
      tmp = uprev+dt*half*k1
      k2 = f(tmp, p, t+half*dt)
      tmp = uprev+dt*half*k2
      k3 = f(tmp, p, t+half*dt)
      tmp = uprev+dt*k3
      k4 = f(tmp, p, t+dt)
      u = uprev+dt*sixth*(k1+2k2+2k3+k4)
      save_everystep && (us[i] = u)
  end

  !save_everystep && save_end && (us[end] = u)

  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  k = nothing, destats = nothing,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end

# In-place
# Good for mutable objects like arrays
# Use DiffEqBase.@.. for simd ivdep
@muladd function DiffEqBase.solve(prob::ODEProblem{uType,tType,true},
                          alg::LoopRK4;
                          dt = error("dt is required for this algorithm"),
                          save_everystep = true,
                          save_start = true,
                          adaptive = false,
                          dense = false,
                          save_end = true,
                          kwargs...) where {uType,tType}
  @assert !adaptive
  @assert !dense
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  t = tspan[1]
  tf = prob.tspan[2]
  ts = tspan[1]:dt:tspan[2]

  if save_everystep && save_start
      us = Vector{typeof(u0)}(undef,length(ts))
      us[1] = u0
  elseif save_everystep
      us = Vector{typeof(u0)}(undef,length(ts)-1)
  elseif save_start
      us = Vector{typeof(u0)}(undef,2)
      us[1] = u0
  else
      us = Vector{typeof(u0)}(undef,1) # for interface compatibility
  end

  u = copy(u0); uprev = copy(u0)
  k1 = zero(u0); k2 = zero(u0); k3 = zero(u0); k4 = zero(u0)
  half = convert(eltype(u0),1//2)
  sixth = convert(eltype(u0),1//6)

  for i in 2:length(ts)
      uprev .= u; t = ts[i]
      f(k1,u,p,t)
      DiffEqBase.@.. u = uprev+dt*half*k1
      f(k2, u, p, t+half*dt)
      DiffEqBase.@.. u = uprev+dt*half*k2
      f(k3, u, p, t+half*dt)
      DiffEqBase.@.. u = uprev+dt*k3
      f(k4, u, p, t+dt)
      DiffEqBase.@.. u = uprev+dt*sixth*(k1+2k2+2k3+k4)
      save_everystep && (us[i] = copy(u))
  end

  !save_everystep && save_end && (us[end] = u)

  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  k = nothing, destats = nothing,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end


text:
# Containers for SCIP structures
import Base: pointer

{% for type_name in parser.structs|sort %}type {{ type_name }}_t
    array_ptr_{{ type_name|lower }}::{{ 'Array{Ptr{%s}}'|format(type_name) }}
end
{% endfor %}
# Pointer and array access functions
{% for type_name in parser.structs|sort %}array({{ type_name|lower }}::{{ type_name }}_t) = {{ type_name|lower }}.array_ptr_{{ type_name|lower }}
pointer({{ type_name|lower }}::{{ type_name }}_t) = array({{ type_name|lower }})[1]
{% endfor %}


text:
# https://www.codewars.com/kata/54ff3102c1bad923760001f3

strtochars = Base.Fix2(split, "")
isvowel = Base.Fix2(occursin, "aeiou")
onlyvowels = Base.Fix1(filter, isvowel)

getcount(str::String) =
  str |>
  strtochars |>
  onlyvowels |>
  length

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
module DashTextareaAutocomplete
using Dash, Pkg.Artifacts

resources_path() = artifact"dash_textarea_autocomplete_resources"
const version = "1.3.2"

include("jl/dashtextareaautocomplete.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_textarea_autocomplete",
            resources_path(),
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_textarea_autocomplete.min.js",
    external_url = "https://unpkg.com/dash-textarea-autocomplete@1.3.2/dash_textarea_autocomplete/dash_textarea_autocomplete.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_textarea_autocomplete.min.js.map",
    external_url = "https://unpkg.com/dash-textarea-autocomplete@1.3.2/dash_textarea_autocomplete/dash_textarea_autocomplete.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end


text:
module EqualitySaturation

using Yolk
import Yolk: allow, YolkOptimizer
using InteractiveUtils

using Metatheory
using Metatheory.Library
using Metatheory.EGraphs

using BenchmarkTools

@metatheory_init()

#####
##### Exceedingly simple.
#####

f(x) = begin
    z = (x - x) + (10 * 15)
    y = z + 10
    return y
end

#####
##### Define theory.
#####

# Automatic.
fold = @theory begin
    Base.:(+)(a::Number, b::Number) |> a + b
end
display(fold)

# Manual.
th = []
left = Metatheory.PatTerm(:call, [PatLiteral(Base.:(-)), PatVar(:a), PatVar(:a)])
right = Metatheory.PatLiteral(0)
new = Metatheory.RewriteRule(left, right)
push!(th, new)
left = Metatheory.PatTerm(:call, [PatLiteral(Base.:(+)), PatTypeAssertion(PatVar(:a), Number), PatTypeAssertion(PatVar(:b), Number)])
right =  fold[1].right
new = Metatheory.DynamicRule(left, right)
push!(th, new)
left = Metatheory.PatTerm(:call, [PatLiteral(Base.:(*)), PatTypeAssertion(PatVar(:a), Number), PatTypeAssertion(PatVar(:b), Number)])
new = Metatheory.DynamicRule(left, Expr(:call, Base.:(*), :a, :b))
push!(th, new)
display(th)

# Optimize.
println("Pre-opt:")
src = opt(f, Tuple{Int}; ctx = YolkOptimizer(th), opt = false)
display(src)
println("Julia native opt:")
src = opt(f, Tuple{Int}; ctx = YolkOptimizer(th, false), opt = true)
display(src)
println("(Metatheory) + Julia native opt:")
src = opt(f, Tuple{Int}; ctx = YolkOptimizer(th), opt = true)
display(src)
end # module


text:
using ApproxFun

#The following solves the Airy ODE with dirichlet boundary conditions

x=Fun(identity,[-1000.,15.])   # Fun corresponding to multiplication by x, on [-100,15]
d=domain(x)
D=Derivative()             # The derivative operator
B=dirichlet()              # Dirichlet boundary conditions, [u(-100),u(15)]

#Construct operator

A=[B;D^2-x]                # This is dirichlet conditions and u'' - x u
b=[airyai(first(d)),0.]         # We want it to equal airyai(-100) at -100, and 0 at
                           # 10, with 0 rhs
#Solve ODE

u=A\b                      # u satisfies A*u = b, or in other words,
                           # B*u = [airyai(-100.),0.] and (D2 - x)*u = 0.

# Check the accuracy
norm(u - Fun(airyai,d))

## We now solve with Neumann conditions

B=neumann()
A=[B;D^2-x]
b=[airyaiprime(first(d)),0.,0.]

u=A\b

# Check the accuracy
norm(u - Fun(airyai,d))

max_stars_repo_path: examples/Airy equation.jl

max_issues_repo_path: examples/Airy equation.jl

max_forks_repo_path: examples/Airy equation.jl


text:
#=

    DeterministicModel.jl

    Perform inference using the deterministic model. I.e., take the 'typical' approach
    and neglect parameter heterogeneity.

=#

using Inference
using Model

using CSV, DataFrames, DataFramesMeta
using Distributions
using ForwardDiff
using JLD2
using LinearAlgebra
using Optim
using StatsBase
using StatsPlots

##############################################################
## Data
##############################################################

# Load all data
data = CSV.read("Data/CSV/DualMarker.csv",DataFrame)

# Load summarised data
data_summary = @subset(CSV.read("Data/CSV/DualMarker_Summary.csv",DataFrame), :Temp .== 37.0, :Time .> 0.0)

# Observation times
T = sort(unique(data_summary.Time))

# Get pooled standard deviation
σ = pooled_std(data_summary,[:Quenched,:Time],[:Signal_Cy5])[1]

# Get quenching efficiency
η = 1 - (@df @subset(data,:Time .== 120.0,:Temp .== 4.0,:Quenched) geomean(:Signal_Cy5)) /
        (@df @subset(data,:Time .== 120.0,:Temp .== 4.0,(!).(:Quenched)) geomean(:Signal_Cy5))

# Get autofluorescence
e = @df @subset(data,:Time .== 0.0, :Quenched .== false, :Temp .== 37.0) geomean(:Signal_Cy5)

##############################################################
## Fit ODE model to data
##############################################################

# Likelihood
function loglike(θ)
    loglike = 0.0
    for t in T
        # Solve ODE model
        A,I = solve_ode_model(t,θ[1:3])
        # Predicted unquenched signal
        m_unquenched = A + e / θ[4]
        # Predicted quenched signal
        m_quenched = I + (1 - η) * (A - I) + e / θ[4]
        # Data ~ N(model,σ²)
        loglike += loglikelihood(Normal(θ[4]*m_unquenched,σ),@subset(data_summary,:Time .== t,(!).(:Quenched)).Signal_Cy5)
        loglike += loglikelihood(Normal(θ[4]*m_quenched,σ),@subset(data_summary,:Time .== t,:Quenched).Signal_Cy5)
    end
    return loglike
end

# Obtain maximum likelihood estimate
θdet = optimize(p -> -loglike(p),ones(4)).minimizer

# Convert to named tuple for easy manipulation
θdet = NamedTuple{(:λ,:β,:p,:α)}(θdet)

##############################################################
## Confidence intervals
##############################################################

# Use observed FIM to approximate confidence intervals
H = ForwardDiff.hessian(loglike,collect(θdet))
se = sqrt.(diag(inv(-H)))
ci = [θdet[i] .+ [-1.0,1.0] * se[i] * quantile(Normal(),0.975) for i = 1:4]

# Tabulate estimates
tab1 = DataFrame(
    Parameter = [:λ,:β,:p,:α],
    Estimate  = round.(collect(θdet),sigdigits=3),
    CI95      = [round.(ci[i],sigdigits=3) for i = 1:length(ci)]
)

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
export NetworkInput,
       NetworkSFunc

struct NetworkInput{T} end

tupler(x::Array) = length(x) == 1 ? x[1] : tuple(x...)
tupler(x::Tuple) = length(x) == 1 ? x[1] : x
tupler(x) = x

get_place_type(p::NetworkInput{T}) where T = T

"""
    struct NetworkSFunc{I,O} <: SFunc{I,O}

An sfunc that combines multiple sfuncs in a network structure.

# Arguments
    input_placeholders A vector of placeholders indicating the types of network inputs. The type parameter `I` is computed from these.
    sfuncs The sfuncs to combine.
    parents A `Dict` that maps sfuncs to their parent sfuncs. Note that this parallels networks,
         except that we are mapping sfuncs to lists of sfuncs directly rather than variables to vectors of variables.
    output A vector of output sfuncs, determining the `O` type parameter.

# Additional supported operators
- `sample`
- `sample_logcpdf`
- `logcpdf`
"""
struct NetworkSFunc{I,O} <: SFunc{I,O}
    input_placeholders::NTuple{N,NetworkInput} where N
    sfuncs::NTuple{N,SFunc} where N
    parents::Dict{SFunc,Vector} where M
    outputs::Tuple
    """
        function NetworkSFunc(input_placeholders, sfuncs, parents, outputs)

    TODO
    """
    function NetworkSFunc(input_placeholders,
                          sfuncs,
                          parents,
                          outputs)
        in_types = [get_place_type(placeholder) for placeholder in input_placeholders]
        # I = length(in_types) == 1 ? in_types[1] : Tuple{in_types...}
        I = Tuple{in_types...}

        out_types = [output_type(sf) for sf in outputs]
        O = length(out_types) == 1 ? out_types[1] : Tuple{out_types...}

        # TODO (MRH): Validate types of parents
        # TODO (MRH): Propagate param types of constituent SFuncs to type parameter P
        return new{I,O}(input_placeholders,
                                sfuncs,
                                parents,
                                outputs)
    end
end

@impl begin
    struct NetworkSFuncSample end
    function sample(sf::NetworkSFunc{I,O}, input::I)::O where {I,O}
        network = sf
        sample_cache = Dict{Union{SFunc,NetworkInput},Any}(s_inp => walk_inp for (s_inp, walk_inp) in zip(network.input_placeholders, input))

        # Assume network.funcs is topologically sorted. TODO?
        for sfunc in network.sfuncs
            sample_cache[sfunc] = sample(sfunc, Tuple([sample_cache[sinp] for sinp in network.parents[sfunc]]))
        end

        s = [sample_cache[o] for o in network.outputs]
        return tupler(s)
    end
end

@impl begin
    struct NetworkSFuncSampleLogcpdf end
    function sample_logcpdf(sf::NetworkSFunc{I,O}, input::I)::Tuple{O, AbstractFloat} where {I,O}
        network = sf
        sample_cache = Dict{Union{SFunc,NetworkInput},Any}(s_inp => (walk_inp, 0.0) for (s_inp, walk_inp) in zip(network.input_placeholders, input))

        # Assume network.funcs is topologically sorted. TODO?
        for sfunc in network.sfuncs
            par_samples = [sample_cache[sinp][1] for sinp in network.parents[sfunc]]
            this_sample, logcpdf = sample_logcpdf(sfunc, Tuple(par_samples))
            logcpdf = logcpdf + sum([sample_cache[sinp][2] for sinp in network.parents[sfunc]])
            sample_cache[sfunc] = (this_sample, logcpdf)
        end
        joint_sample = tuple([sample_cache[o][1] for o in network.outputs]...)
        joint_logcpdf = sum([sample_cache[o][2] for o in network.outputs])

        return (tupler(joint_sample), joint_logcpdf)
    end
end

@impl begin
    struct NetworkSFuncLogcpdf end
    function logcpdf(sf::NetworkSFunc{I,O}, input::I, output::O)::AbstractFloat where {I,O}
        # I think this returns a sample x of a distribution s.t. log(Expectation(exp(x))) gives the logcpdf
        # In the special case where all NetworkSFunc nodes are in outputs then the calculation is deterministic

        sample_cache = merge(Dict{Union{SFunc,NetworkInput},Any}(s_inp => (walk_inp, 0.0) for (s_inp, walk_inp) in zip(sf.input_placeholders, input)),
                        Dict{Union{SFunc,NetworkInput},Any}(s_out => (walk_out, nothing) for (s_out, walk_out) in zip(sf.outputs, output)))

        # Assume network.funcs is topologically sorted. TODO?

        # Is this even right? Seems neat
        for sfunc in sf.sfuncs
            if haskey(sample_cache, sfunc)
                par_samples = [sample_cache[sinp][1] for sinp in sf.parents[sfunc]]
                cumlogcpdf = logcpdf(sfunc, Tuple(par_samples), sample_cache[sfunc][1]) +
                        sum([sample_cache[sinp][2] for sinp in sf.parents[sfunc]])
                sample_cache[sfunc] = (sample_cache[sfunc][1], cumlogcpdf)
            else
                par_samples = [sample_cache[sinp][1] for sinp in sf.parents[sfunc]]
                this_sample = sample(sfunc, Tuple(par_samples))
                cumlogcpdf = sum([sample_cache[sinp][2] for sinp in sf.parents[sfunc]])
                sample_cache[sfunc] = (this_sample, cumlogcpdf)
            end
        end

        return sum([sample_cache[o][2] for o in sf.outputs])
    end

end


text:
using Plots, Statistics, StatPlots
pyplot(show=true, reuse=false)
# map R^3 to param space
function θ2ϕ(θ)
    ϕ = similar(θ)
    ϕ[1] = exp(θ[1])
    ϕ[2] = 1.0 / (1.0 + exp(-θ[2]))
    ϕ[3] = exp(θ[3])
    return ϕ
end
# map param space to R^3
function ϕ2θ(ϕ)
    θ = similar(ϕ)
    θ[1] = log(ϕ[1])
    θ[2] = -log(1.0/ϕ[2] - 1.0)
    θ[3] = log(ϕ[3])
    return θ
end
# simple discrete time SV model
function SVmodel(n, θ, randdraws = 0.0)
    ϕ = θ2ϕ(θ)
    α = ϕ[1] # mean of log variance: positive
    ρ = ϕ[2] # AR of log variance: between 0 and 1
    σᵤ = ϕ[3] # std. dev. of log variance: positive
    burnin = 1000
    ys = zeros(n)
    hlag = 0.0
    # generate shocks if not provided
    if randdraws == 0.0
        randdraws = randn(burnin+n,2)
    end
    @inbounds for t = 1:burnin+n
        h = α + ρ*(hlag-α) + σᵤ*randdraws[t,1] # log variance follows AR(1)
        y = exp(h/2.0)*randdraws[t,2]
        if t > burnin
            ys[t-burnin] = y
        end
        hlag = h
    end
    σ = exp(hlag/2.0) # the latent volatility in last period
    return ys, σ
end

# auxiliary model: HAR-RV(p)
# Corsi, Fulvio. "A simple approximate long-memory model
# of realized volatility." Journal of Financial Econometrics 7,
# no. 2 (2009): 174-196.
function HAR(y,p)
    σhat = std(y)
    y ./= σhat
    RV = abs.(y)
    RVlags = lags(RV,p)
    X = [ones(size(y,1)) RVlags]
    # drop missings
    RV = RV[p+1:end]
    X = X[p+1:end,:]
    βhat = X\RV
    return vcat(βhat,σhat)
end

function II_moments(θ, ϕhat, n, randdraws)
    S = size(randdraws,1)
    n_auxparams = size(ϕhat,1)
    p = n_auxparams - 2 # number of lags in HAR
    ϕhatS = zeros(S, n_auxparams)
    for s = 1:S
        yₛ, junk = SVmodel(n, θ, randdraws[s,:,:])
        ϕhatS[s,:] = HAR(yₛ,p)
    end
    ms = ϕhat' .- ϕhatS
    return  ms # the moments, in a SxG matrix
end

function main()
# generate the sample
ϕ₀ = [0.7, 0.95, 0.2] # true param values, on param space
θₒ  = ϕ2θ(ϕ₀) # true parameter values, on R^2
n = 1000 # sample size
y, σ = SVmodel(n, θₒ) # generate the sample
plot(y)
density(y)

# Estimation by indirect inference
S = 20 # number of simulation reps
randdraws = randn(S,n+1000,2) # fix the shocks to control "chatter" (includes the burnin period)
p = 4
ϕhat = HAR(y,p) # statistics using real sample
moments = θ -> II_moments(θ, ϕhat, n, randdraws)
thetahat, junk, junk, junk, junk = gmm(moments, θₒ, 1.0)
gmmresults(moments, thetahat, "", "This minimizes m(θ)'inv[Ω(θ)]m(θ), but std errs, etc. are not correct")
θ2ϕ(thetahat)
end
main()


text:
using IJulia
using Reactive, Interact

notebook();

min_val = slider(-300:1:300,label="Min value:")
max_val = slider(-300:1:300,label="Max value:")
points_number = slider(1:1:100,label="Points number:")

include("convex_hull.jl")

points = generate_points(signal(points_number), signal(min_val), signal(max_val))
main(points)

max_stars_repo_name: MrShemek/ConvexHullTDD

max_issues_repo_name: MrShemek/ConvexHullTDD

max_forks_repo_name: MrShemek/ConvexHullTDD


text:
using RoundedNumbersForHumans
using Test
import RoundedNumbersForHumans: Scientific, Common

@testset "string" begin
    function test_common(n, str_common)
        @test string(RoundedNumber(n)) == str_common
        @test string(RoundedNumber(n, Common())) == str_common
    end
    function test_common_sci(n, str_common, str_sci = str_common)
        test_common(n, str_common)
        @test string(RoundedNumber(n, Scientific())) == str_sci
    end

    test_common_sci(300, "300")
    test_common_sci(3_000, "3k")
    test_common_sci(3_000_000, "3M")
    test_common_sci(3_000_000_000, "3B", "3G")
    if Int === Int64
        test_common_sci(3_000_000_000_000, "3T")
        test_common_sci(3_000_000_000_000_000, "3Qd", "3P")
        test_common_sci(3_000_000_000_000_000_000, "3Qn", "3E")
    end

    @testset "sigdigits" begin
    	@testset "default" begin
            test_common(3_100, "3.1k")
            test_common(3_150, "3.15k")
            test_common(3_153, "3.15k")
            test_common(31_532, "31.5k")
    	end
    	@testset "specified" begin
		    @test string(RoundedNumber(31_532),sigdigits=1) == "30k"
		    @test string(RoundedNumber(31_532),sigdigits=2) == "32k"
		    @test string(RoundedNumber(31_532),sigdigits=3) == "31.5k"
		    @test string(RoundedNumber(31_532),sigdigits=4) == "31.53k"
		    @test string(RoundedNumber(31_532),sigdigits=5) == "31.532k"
    	end
    end

    @testset "float" begin
        test_common_sci(300.3, "300")
        test_common_sci(300.6, "301")
        test_common_sci(3_000.6, "3k")
        @test string(RoundedNumber(3_000.6),sigdigits=4) == "3.001k"

        @testset "negative exponent" begin
            test_common_sci(10^-1, "100m")
            test_common_sci(10^-12, "1p")
            test_common_sci(10^-24, "1.0e-6a", "1y")
        end
    end

    @testset "zero" begin
        test_common_sci(0, "0")
    end

    @testset "negative" begin
        test_common_sci(-30, "-30")
        test_common_sci(-3000, "-3k")
        test_common_sci(-30_000, "-30k")
    end

end

@testset "convert" begin
    for n in [100, 100_000, 100_000_000]
        r = RoundedNumber(n)
        @test Int(r) === n
        @test BigInt(r) == big(n)
        @test float(r) === float(n)
    end
end


text:
using MAT, Meganet
# BLAS.set_num_threads(1)
n = 512;
Y_train,C_train,Y_test,C_test = getCIFAR10(n,Pkg.dir("Meganet")*"/data/CIFAR10/");

# using PyPlot
# y = Y_train[:,50]; y = y - minimum(y); y = y./maximum(y);
# y = reshape(y,32,32,3);
# y[:,:,1] = y[:,:,1]';y[:,:,2] = y[:,:,2]';y[:,:,3] = y[:,:,3]';
# figure(); imshow(y)

miniBatchSize = 64;
nImg = [32; 32]
cin  = 3
nc   = [16;32;64;64]
nt   = 2*[1;1;1]
h    = [1.;1.;1.]

TYPE = Float32;

getConvKernel = (nImg,sK) -> getConvGEMMKernel(TYPE,nImg,sK);
# getConvKernel = (nImg,sK) -> getConvFFTKernel(TYPE,nImg,sK);
#getConvKernel = (nImg,sK) -> getSparseConvKernel2D(TYPE,nImg,sK);

# opening layer
K1 = getConvKernel(nImg,[3,3,cin,nc[1]]);

nL = getBatchNormLayer(TYPE,[prod(nImg);nc[1]],isTrainable=true);
blocks = [getSingleLayer(TYPE,K1,nL)]

for k=1:length(nt)
    # ResNN layers
    K2 = getConvKernel(nImg,[3,3,nc[k],nc[k]])
    nL = getBatchNormLayer(TYPE,[prod(nImg);nc[k]],isTrainable=true)
    L2 = getDoubleSymLayer(TYPE,K2,nL)
    RN  = getResNN(TYPE,L2,nt[k],h[k])

    if k<length(nt)
        RN.outTimes *=0
    end
    blocks = [blocks;RN]
    # change channels
    Kc = getConvKernel(nImg,[1,1,nc[k],nc[k+1]]);

    nL = getBatchNormLayer(TYPE,[prod(nImg);nc[k+1]],isTrainable=true)
    blocks = [blocks; getSingleLayer(TYPE,Kc,nL)]

    if k<length(nt)
        Kp = getAverageMatrix(TYPE,nImg,nc[k+1])
        blocks = [blocks; getConnector(TYPE,Kp)]
        # nImg ./=2
		nImg = div.(nImg,2)
    end
end

# Connector block
B      = kron(speye(TYPE,nc[end]),ones(TYPE, prod(nImg)))/prod(nImg);
blocks = [blocks; getConnector(TYPE,B')];
blocks[end].outTimes=1

net = getNN(blocks)
theta = initTheta(net);

display(net)

# @time Zj = apply(net,theta,Y_train[:,1:2],true)
# @time Zj = apply(net,theta,Y_train[:,1:miniBatchSize],true)

# regularizers
pRegTh = getTikhonovReg(TYPE;alpha=4e-4)
pRegW = getTikhonovReg(TYPE;alpha=4e-4)
pLoss = getSoftMaxLoss(TYPE);
objFun = dnnObjFctn(net,pLoss,pRegTh,pRegW)
opt = getSGDsolver(TYPE,learningRate=1e-2,maxEpochs=1,miniBatch=miniBatchSize,out=true, nesterov=true)

W      = 0.1*vec(randn(TYPE,10,nFeatOut(net)+1));
W = min.(W,.2);
W = max.(W,-.2);
W = convert(Array{TYPE},W);

solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)
@time solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)

# Profile.clear()
# Profile.clear_malloc_data()
# Profile.init(n = 10^7, delay = 0.01)
# @profile solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y_train,C_train,Y_test,C_test)

# open("/tmp/EREsNN_CIFAR10.txt", "w") do s
#     Profile.print(IOContext(s, :displaysize => (24, 500)))
# end


text:
"""
    runexamples()

Run all *_examples.jl files in the working folder.
"""
function runexamples()
    origcwd = pwd()

    pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))

    println("\nExamples in folder $(pwd())\n")

    for ex1 in readdir(".")
        if occursin(r".*_examples.jl", ex1)
            println("\nExample $ex1 in $(pwd())\n")
            include(ex1);
            n, ext = splitext(ex1)
            Meta.eval(Meta.parse("Main." * n * ".allrun()"))
        end
    end

    cd(origcwd)
end


text:
function calc_tderivative!(integrator, cache, dtd1, repeat_step)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack du2,fsalfirst,dT,tf,linsolve_tmp = cache

    # Time derivative
    if !repeat_step # skip calculation if step is repeated
      if DiffEqBase.has_tgrad(f)
        f.tgrad(dT, uprev, p, t)
      else
        tf.uprev = uprev
        tf.p = p
        derivative!(dT, tf, t, du2, integrator, cache.grad_config)
      end
    end

    f(fsalfirst, uprev, p, t)
    @. linsolve_tmp = fsalfirst + dtd1*dT
  end
end

function calc_tderivative(integrator, cache)
  @unpack t,dt,uprev,u,f,p = integrator

  # Time derivative
  if DiffEqBase.has_tgrad(f)
    dT = f.tgrad(uprev, p, t)
  else
    tf = cache.tf
    tf.u = uprev
    tf.p = p
    dT = derivative(tf, t, integrator)
  end
  dT
end

"""
    calc_J!(integrator,cache,is_compos)

Interface for calculating the jacobian.

For constant caches, a new jacobian object is returned whereas for mutable
caches `cache.J` is updated. In both cases, if `integrator.f` has a custom
jacobian update function, then it will be called for the update. Otherwise,
either ForwardDiff or finite difference will be used depending on the
`jac_config` of the cache.
"""
function calc_J(integrator, cache::StochasticDiffEqConstantCache, is_compos)
  @unpack t,dt,uprev,u,f,p = integrator
  if DiffEqBase.has_jac(f)
    J = f.jac(uprev, p, t)
  else
    J = jacobian(cache.uf,uprev,integrator)
  end
  is_compos && (integrator.eigen_est = opnorm(J, Inf))
  return J
end

function calc_J!(integrator, cache::StochasticDiffEqMutableCache, is_compos)
  @unpack t,dt,uprev,u,f,p = integrator
  J = cache.J
  if has_jac(f)
    f.jac(J, uprev, p, t)
  else
    @unpack du1,uf,jac_config = cache
    uf.t = t
    uf.p = p
    jacobian!(J, uf, uprev, du1, integrator, jac_config)
    if is_compos
      integrator.eigen_est = norm(J, Inf)
    end
  end
end

"""
    WOperator(mass_matrix,gamma,J)

A linear operator that represents the W matrix of an ODEProblem, defined as

```math
W = MM - \\gamma J
```

where `MM` is the mass matrix (a regular `AbstractMatrix` or a `UniformScaling`),
`γ` is a real number proportional to the time step, and `J` is the Jacobian
operator (must be a `AbstractDiffEqLinearOperator`). A `WOperator` can also be
constructed using a `*DEFunction` directly as

    WOperator(f,gamma)

`f` needs to have a jacobian and `jac_prototype`, but the prototype does not need
to be a diffeq operator --- it will automatically be converted to one.

`WOperator` supports lazy `*` and `mul!` operations, the latter utilizing an
internal cache (can be specified in the constructor; default to regular `Vector`).
It supports all of `AbstractDiffEqLinearOperator`'s interface.
"""
mutable struct WOperator{T,
  MType <: Union{UniformScaling,AbstractMatrix},
  GType <: Real,
  JType <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  } <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  mass_matrix::MType
  gamma::GType
  J::JType
  _func_cache           # cache used in `mul!`
  _concrete_form        # non-lazy form (matrix/number) of the operator
  WOperator(mass_matrix, gamma, J) = new{eltype(J),typeof(mass_matrix),
    typeof(gamma),typeof(J)}(mass_matrix,gamma,J,nothing,nothing)
end
function WOperator(f::DiffEqBase.AbstractSDEFunction, gamma)
  @assert DiffEqBase.has_jac(f) "f needs to have an associated jacobian"
  if isa(f, SplitSDEFunction)
    error("WOperator does not support SplitSDEFunction yet")
  end
  # Convert mass matrix, if needed
  mass_matrix = f.mass_matrix
  if !isa(mass_matrix, Union{AbstractMatrix,UniformScaling})
    mass_matrix = convert(AbstractMatrix, mass_matrix)
  end
  # Convert jacobian, if needed
  J = deepcopy(f.jac_prototype)
  if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
    J = DiffEqArrayOperator(J; update_func=f.jac)
  end
  return WOperator(mass_matrix, gamma, J)
end

set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W)
DiffEqBase.update_coefficients!(W::WOperator,u,p,t) = (update_coefficients!(W.J,u,p,t); W)
function Base.convert(::Type{AbstractMatrix}, W::WOperator)
  if W._concrete_form === nothing
    # Allocating
    W._concrete_form = W.mass_matrix - W.gamma * convert(AbstractMatrix,W.J)
  else
    # Non-allocating
    copyto!(W._concrete_form, W.mass_matrix)
    axpy!(-W.gamma, convert(AbstractMatrix,W.J), W._concrete_form)
  end
  W._concrete_form
end
Base.convert(::Type{Number}, W::WOperator) = W.mass_matrix - W.gamma * convert(Number,W.J)
Base.size(W::WOperator, args...) = size(W.J, args...)
Base.getindex(W::WOperator, i::Int) = W.mass_matrix[i] - W.gamma * W.J[i]
Base.getindex(W::WOperator, I::Vararg{Int,N}) where {N} =
  W.mass_matrix[I...] - W.gamma * W.J[I...]
Base.:*(W::WOperator, x::Union{AbstractVecOrMat,Number}) = W.mass_matrix*x - W.gamma * (W.J*x)
function Base.:\(W::WOperator, x::Union{AbstractVecOrMat,Number})
  if size(W) == () # scalar operator
    convert(Number,W) \ x
  else
    convert(AbstractMatrix,W) \ x
  end
end
function LinearAlgebra.mul!(Y::AbstractVecOrMat, W::WOperator, B::AbstractVecOrMat)
  if W._func_cache === nothing
    # Allocate cache only if needed
    W._func_cache = Vector{eltype(W)}(undef, size(Y, 1))
  end
  # Compute mass_matrix * B
  if isa(W.mass_matrix, UniformScaling)
    @. Y = W.mass_matrix.λ * B
  else
    mul!(Y, W.mass_matrix, B)
  end
  # Compute J * B
  mul!(W._func_cache, W.J, B)
  # Subtract result
  axpy!(-W.gamma, W._func_cache, Y)
end

function calc_W!(integrator, cache::StochasticDiffEqMutableCache, γdt, repeat_step)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack J,W = cache
    is_compos = is_composite(integrator.alg)
    alg = unwrap_alg(integrator, true)
    mass_matrix = integrator.f.mass_matrix
    nlcache = cache.nlsolve.cache

    new_W = true
    if has_invW(f)
      # skip calculation of inv(W) if step is repeated
      !repeat_step && f.invW(W,uprev,p,γdt,t) # W == inverse W
      is_compos && calc_J!(integrator, cache, true)
    elseif has_jac(f) && f.jac_prototype !== nothing
      # skip calculation of J if step is repeated
      if repeat_step || (!integrator.last_stepfail && nlcache.nl_iters == 1 && nlcache.ηold < alg.new_jac_conv_bound)
        new_jac = false
      else # Compute a new Jacobian
        new_jac = true
        DiffEqBase.update_coefficients!(W,uprev,p,t)
      end
      # skip calculation of W if step is repeated
      if !repeat_step && (integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps(typeof(integrator.t)))
        set_gamma!(W, γdt)
      else
        new_W = false
      end
    else
      # skip calculation of J if step is repeated
      if repeat_step || (!integrator.last_stepfail && nlcache.nl_iters == 1 && nlcache.ηold < alg.new_jac_conv_bound)
        new_jac = false
      else # Compute a new Jacobian
        new_jac = true
        calc_J!(integrator, cache, is_compos)
      end
      # skip calculation of W if step is repeated
      if !repeat_step && (integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps(typeof(integrator.t)))
        for j in 1:length(u), i in 1:length(u)
            @inbounds W[i,j] = mass_matrix[i,j]-γdt*J[i,j]
        end
      else
        new_W = false
      end
    end
    nlcache.new_W = new_W
    return nothing
  end
end

function calc_W!(integrator, cache::StochasticDiffEqConstantCache, γdt, repeat_step)
  @unpack t,uprev,p,f = integrator
  uf = cache.uf
  isarray = typeof(uprev) <: AbstractArray
  is_compos = is_composite(integrator.alg)
  mass_matrix = integrator.f.mass_matrix
  if has_jac(f)
    J = f.jac(uprev, p, t)
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(mass_matrix, γdt, J)
  else
    J = calc_J(integrator, cache, is_compos)
    W = mass_matrix - γdt*J
  end
  is_compos && (integrator.eigen_est = isarray ? opnorm(J, Inf) : J)
  J, W
end


text:
# These types are defined in grampc_macro.h, and are used as aliases
# here to allow for easier copy-pasting of the GRAMPC C structures into
# the Julia code.
const typeRNum      = Cdouble
const ctypeRNum     = Cdouble  # This is const in C, but Julia doesn't care about that
const typeInt       = Cint
const ctypeInt      = Cint     # This is const in C, but Julia doesn't care about that
const typeLInt      = Cint
const typeLogical   = Clong
const typeBoolean   = Cint
const typeChar      = Cchar
const typeUSERPARAM = Cvoid

# These status codes are in grampc_mess.h
@bitflag GRAMPCStatus::UInt32 begin
    NONE                            =    0
    GRADIENT_CONVERGED              =    1
    CONSTRAINTS_CONVERGED           =    2
    USED_LINESEARCH_MIN             =    4
    USED_LINESEARCH_MAX             =    8
    USED_LINESEARCH_INIT            =   16
    MULTIPLIER_UPDATED              =   32
    REACHED_MULTIPLIER_MAX          =   64
    REACHED_PENALTY_MAX             =  128
    INFEASIBLE                      =  256
    INTEGRATOR_INPUT_NOT_CONSISTENT =  512
    INTEGRATOR_REACHED_MAXSTEPS     = 1024
    INTEGRATOR_STEPS_TOO_SMALL      = 2048
    INTEGRATOR_MATRIX_IS_SINGULAR   = 4096
    INTEGRATOR_USED_H_MIN           = 8192
end

@enum PenaltyEstimationStatus begin
    SATURATED_PENALTY         = 0
    USED_ACTIVE_CONSTRAINTS   = 1
    USED_CONSTRAINT_TOLERANCE = 2
end

# These values come from the grampc_macro.h file

# Boolean map for the fields that are binary options
const boolMap = Dict{Int, Bool}(
    0 => false,
    1 => true )

# The time discretizations for the "TimeDiscretization" setting field
const timeDiscretization = Dict{Int, Symbol}(
    0 => :uniform,
    1 => :nonuniform )

const timeDiscretizationStrings = Dict{Symbol, String}(
    :uniform => "uniform",
    :nonuniform => "nonuniform" )

# The integrator options for the "Integrator" field
const dynamicsIntegrators = Dict{Int, Symbol}(
    0 => :euler,
    1 => :modeuler,
    2 => :heun,
    3 => :rodas,
    4 => :rk45 )

const dynamicsIntegratorsStrings = Dict{Symbol, String}(
    :euler => "euler",
    :modeuler => "modeuler",
    :heun => "heun",
    :rodas => "rodas",
    :rk45 => "ruku45" )

# The integrator options for the "IntegratorCost" field
const costIntegrators = Dict{Int, Symbol}(
    0 => :trapezoidal,
    1 => :simpson )

const costIntegratorsStrings = Dict{Symbol, String}(
    :trapezoidal => "trapezoidal",
    :simpson => "simpson" )

# The line searches available for the "LineSearchType" field
const lineSearchType = Dict{Int, Symbol}(
    0 => :adaptive,
    1 => :explicit_1,
    2 => :explicit_2 )

const lineSearchTypeStrings = Dict{Symbol, String}(
    :adaptive => "adaptive",
    :explicit_1 => "explicit1",
    :explicit_2 => "explicit2" )

# The method for handling the constraints for the "ConstraintsHandling" field
const constraintHandlingMethod = Dict{Int, Symbol}(
    0 => :external_penalty,
    1 => :augmented_lagrangian )

const constraintHandlingMethodStrings = Dict{Symbol, String}(
    :external_penalty => "extpen",
    :augmented_lagrangian => "auglag" )

# This dictionary maps the field to its setting strings
const optionSymbolStringDict = Dict{Symbol, Dict}(
    :TimeDiscretization => timeDiscretizationStrings,
    :IntegratorCost => costIntegratorsStrings,
    :Integrator => dynamicsIntegratorsStrings,
    :LineSearchType => lineSearchTypeStrings,
    :ConstraintsHandling => constraintHandlingMethodStrings )

const optionSymbolSymbolDict = Dict{Symbol, Dict}(
    :TimeDiscretization => timeDiscretization,
    :IntegratorCost => costIntegrators,
    :Integrator => dynamicsIntegrators,
    :LineSearchType => lineSearchType,
    :ConstraintsHandling => constraintHandlingMethod )


text:
function nonnaivebanks(mdl::Model)

    decision = view(mdl.DR,:,:)
    value = view(mdl.VF,:,:)
    default_id =  [0] #[-(1:size(mdl.VF,1))...]
    check = true

    while check

        tmp = LinearIndices( decision[:,1] )[findall( all(decision .∈ [default_id], dims = 2 ) )]
        check = !(all(tmp .∈ [default_id]))
        check && append!(default_id, tmp[(tmp .∈ [default_id]) .== false])

    end

    if (length(default_id) > 1)

        decision[default_id[2:end],:] .= 0
        value[default_id[2:end],:] .= mdl.VFalter[default_id[2:end],:]

    end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Structures related to the IERS EOP support.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export EOPData_IAU1980, EOPData_IAU2000A

"""
    EOPData_IAU1980{T}

EOP Data for IAU 1980.

!!! note
    Each field will be an `AbstractInterpolation` indexed by the Julian Day.
    Hence, if one wants to obtain, for example, the X component of the polar
    motion with respect to the crust at 19 June 2018, the following can be used:

        x[date_to_jd(2018, 19, 06, 0, 0, 0)]

# Fields

- `x, y`: Polar motion with respect to the crust [arcsec].
- `UT1_UTC`: Irregularities of the rotation angle [s].
- `LOD`: Length of day offset [ms].
- `dPsi, dEps`: Celestial pole offsets referred to the model IAU1980
    [milliarcsec].
- `*_err`: Errors in the components [same unit as the component].
"""
struct EOPData_IAU1980{T}
    x::T
    y::T
    UT1_UTC::T
    LOD::T
    dPsi::T
    dEps::T

    # Errors
    x_err::T
    y_err::T
    UT1_UTC_err::T
    LOD_err::T
    dPsi_err::T
    dEps_err::T
end

"""
    EOPData_IAU2000A{T}

EOP Data for IAU 2000A.

!!! note
    Each field will be an `AbstractInterpolation` indexed by the Julian Day.
    Hence, if one want to obtain, for example, the X component of the polar
    motion with respect to the crust at 19 June 2018, the following can be used:

        x[date_to_jd(2018, 19, 06, 0, 0, 0)]

# Fields

- `x, y`: Polar motion with respect to the crust [arcsec].
- `UT1_UTC`: Irregularities of the rotation angle [s].
- `LOD`: Length of day offset [ms].
- `dX, dY`: Celestial pole offsets referred to the model IAU2000A [milliarcsec].
- `*_err`: Errors in the components [same unit as the component].
"""
struct EOPData_IAU2000A{T}
    x::T
    y::T
    UT1_UTC::T
    LOD::T
    dX::T
    dY::T

    # Errors
    x_err::T
    y_err::T
    UT1_UTC_err::T
    LOD_err::T
    dX_err::T
    dY_err::T
end


text:
function create_testfuncs(g::OneDGrid{Tg,<:Array}) where Tg
    g.nx > 8 || error("nx must be > 8")
    m = 5
    k₀ = g.k[2] # Fundamental wavenumber

    # Analytic function
    φ = π/3
    f₁ = @. cos(m*k₀*g.x + φ)

    # Transform of f₁ three ways
     f₁h = fft(f₁)
    f₁hr = rfft(f₁)
    f₁hr_mul = zeros(Complex{eltype(g.x)}, g.nkr)
    mul!(f₁hr_mul, g.rfftplan, f₁)

    # Analytical values of the fft and rfft
     f₁h_analytical = zeros(Complex{Tg}, size(f₁h))
    f₁hr_analytical = zeros(Complex{Tg}, size(f₁hr))

    for i in 1:g.nk
      if abs(real(g.k[i])) == m*k₀
        f₁h_analytical[i] = -exp(sign(real(g.k[i]))*im*φ)*g.nx/2
      end
    end

    for i in 1:g.nkr
      if abs(real(g.k[i])) == m*k₀
        f₁hr_analytical[i] = -exp(sign(real(g.kr[i]))*im*φ)*g.nx/2
      end
    end

    f₁, f₁h, f₁hr, f₁hr_mul, f₁h_analytical, f₁hr_analytical
end

function create_testfuncs(g::TwoDGrid{Tg,<:Array}) where Tg
    g.nx > 8 || error("nx must be > 8")
    x, y = gridpoints(g)
    nx, ny = g.nx, g.ny
    m, n = 5, 2
    k₀, l₀ = g.k[2], g.l[2]

    # Analytic functions
    f₁ = @. cos(m*k₀*x) * cos(n*l₀*y)
    f₂ = @. sin(m*k₀*x + n*l₀*y)

     f₁h = fft(f₁)
     f₂h = fft(f₂)
    f₁hr = rfft(f₁)
    f₂hr = rfft(f₂)

    f₁hr_mul = zeros(Complex{Tg}, (g.nkr, g.nl))
    f₂hr_mul = zeros(Complex{Tg}, (g.nkr, g.nl))
    mul!(f₁hr_mul, g.rfftplan, f₁)
    mul!(f₂hr_mul, g.rfftplan, f₂)

    # Theoretical results
     f₁h_analytical = zeros(Complex{Tg}, size(f₁h))
     f₂h_analytical = zeros(Complex{Tg}, size(f₂h))
    f₁hr_analytical = zeros(Complex{Tg}, size(f₁hr))
    f₂hr_analytical = zeros(Complex{Tg}, size(f₂hr))

    for j in 1:g.nl, i in 1:g.nk
      if ( abs(g.k[i]) == m*k₀ && abs(g.l[j]) == n*l₀ )
        f₁h_analytical[i, j] = -nx*ny/4
      end
      if ( g.k[i] == m*k₀ && g.l[j] == n*l₀ )
        f₂h_analytical[i, j] = -nx*ny/2
      elseif ( g.k[i] == -m*k₀ && g.l[j] == -n*l₀ )
        f₂h_analytical[i, j] = nx*ny/2
      end
    end
    f₂h_analytical = -im*f₂h_analytical

    for j in 1:g.nl, i in 1:g.nkr
      if ( abs(g.kr[i])==m*k₀ && abs(g.l[j])==n*l₀ )
        f₁hr_analytical[i, j] = - nx*ny/4
      end
      if ( g.kr[i] == m*k₀ && g.l[j] == n*l₀ )
        f₂hr_analytical[i, j] = -nx*ny/2
      end
    end
    f₂hr_analytical = -im*f₂hr_analytical

    f₁, f₂, f₁h, f₂h, f₁hr, f₂hr, f₁hr_mul, f₂hr_mul, f₁h_analytical, f₁hr_analytical, f₂h_analytical, f₂hr_analytical
end

function create_testfuncs(g::ThreeDGrid{Tg,<:Array}) where Tg
    g.nx > 8 || error("nx must be > 8")
    x, y, z = gridpoints(g)
    nx, ny, nz = g.nx, g.ny, g.nz
    mx, my, mz = 5, 2, 3
    k₀, l₀, m₀ = g.k[2], g.l[2], g.m[2]

    # Analytic functions
    f₁ = @. cos(mx*k₀*x) * cos(my*l₀*y) * cos(mz*m₀*z)
    f₂ = @. sin(mx*k₀*x + my*l₀*y + mz*m₀*z)

     f₁h = fft(f₁)
     f₂h = fft(f₂)
    f₁hr = rfft(f₁)
    f₂hr = rfft(f₂)

    f₁hr_mul = zeros(Complex{Tg}, (g.nkr, g.nl, g.nm))
    f₂hr_mul = zeros(Complex{Tg}, (g.nkr, g.nl, g.nm))
    mul!(f₁hr_mul, g.rfftplan, f₁)
    mul!(f₂hr_mul, g.rfftplan, f₂)

    # Theoretical results
     f₁h_analytical = zeros(Complex{Tg}, size(f₁h))
     f₂h_analytical = zeros(Complex{Tg}, size(f₂h))
    f₁hr_analytical = zeros(Complex{Tg}, size(f₁hr))
    f₂hr_analytical = zeros(Complex{Tg}, size(f₂hr))

    for k in 1:g.nm, j in 1:g.nl, i in 1:g.nk
      if ( abs(g.k[i]) == mx*k₀ && abs(g.l[j]) == my*l₀ && abs(g.m[k]) == mz*m₀ )
        f₁h_analytical[i, j, k] = nx*ny*nz/8
      end
      if ( g.k[i] == mx*k₀ && g.l[j] == my*l₀ && g.m[k] == mz*m₀ )
        f₂h_analytical[i, j, k] = nx*ny*nz/2
      elseif ( g.k[i] == -mx*k₀ && g.l[j] == -my*l₀ && g.m[k] == -mz*m₀ )
        f₂h_analytical[i, j, k] = -nx*ny*nz/2
      end
    end
    f₂h_analytical = -im*f₂h_analytical;

    for k in 1:g.nm, j in 1:g.nl, i in 1:g.nkr
      if ( abs(g.kr[i]) == mx*k₀ && abs(g.l[j]) == my*l₀ && abs(g.m[k]) == mz*m₀ )
        f₁hr_analytical[i, j, k] = nx*ny*nz/8
      end
      if ( g.kr[i] == mx*k₀ && g.l[j] == my*l₀ && g.m[k] == mz*m₀ )
        f₂hr_analytical[i, j, k] = nx*ny*nz/2
      end
    end
    f₂hr_analytical = -im*f₂hr_analytical;

    f₁, f₂, f₁h, f₂h, f₁hr, f₂hr, f₁hr_mul, f₂hr_mul, f₁h_analytical, f₁hr_analytical, f₂h_analytical, f₂hr_analytical
end

if CUDA.has_cuda()
  function create_testfuncs(g::OneDGrid{Tg, <:CuArray}) where Tg
    cpugrid = OneDGrid(g.nx, g.Lx)
    out = create_testfuncs(cpugrid)
    return map(x->CuArray(x), out)
  end

  function create_testfuncs(g::TwoDGrid{Tg, <:CuArray}) where Tg
    cpugrid = TwoDGrid(g.nx, g.Lx, g.ny, g.Ly)
    out = create_testfuncs(cpugrid)
    return map(x->CuArray(x), out)
  end

  function create_testfuncs(g::ThreeDGrid{Tg, <:CuArray}) where Tg
    cpugrid = ThreeDGrid(g.nx, g.Lx, g.ny, g.Ly, g.nz, g.Lz)
    out = create_testfuncs(cpugrid)
    return map(x->CuArray(x), out)
  end
end


text:
function test_alias(from, to, ex)
  alias = AliasDeclaration(from => to)
  test_ex(alias, ex)
end

@testset "Alias declarations" begin
  test_alias(:vkCmdSetCullModeEXT, :vkCmdSetCullMode, :(const cmd_set_cull_mode_ext = cmd_set_cull_mode))
  test_alias(:VK_RENDERING_RESUMING_BIT_KHR, :VK_RENDERING_RESUMING_BIT, :(const RENDERING_RESUMING_BIT_KHR = RENDERING_RESUMING_BIT))
  test_alias(:VkPhysicalDeviceImageFormatInfo2KHR, :VkPhysicalDeviceImageFormatInfo2, :(const PhysicalDeviceImageFormatInfo2KHR = PhysicalDeviceImageFormatInfo2))
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using RainFARM, Compat
using Compat.Statistics
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here

print("Testing main RainFARM function and aggregation:\n")
print("-----------------------------------------------\n")
nt=2; nf=8; ns=64; nas=8;
prf=rand(ns,ns,nt);
prf[2,2,1]=NaN;
prl=agg(prf,nas,nt);
ww=ones(64,64);
print("Testing fglob=false, fsmooth=false\n")
pr=rainfarm(prl,1.7,nf,1.;fglob=false, fsmooth=false, verbose=false);
pra=agg(pr,nas,nt);
eps1=Statistics.mean((pra-prl).^2);
@test eps1 < 1e-20
print("Testing fglob=true, fsmooth=false, weights\n")
pr=rainfarm(prl,1.7,nf,ww;fglob=true, fsmooth=false, verbose=false);
pra=agg(pr,nas,nt);
eps1=(Statistics.mean(pra[:,:,1])-Statistics.mean(prl[:,:,1])).^2;
@test eps1 < 1e-20
print("Testing fglob=false, fsmooth=true\n")
pr=rainfarm(prl,1.7,nf,1.;fglob=false, fsmooth=true, verbose=false);
pra=agg(pr,nas,nt);
eps1=(Statistics.mean(pra[:,:,1])-Statistics.mean(prl[:,:,1])).^2;
@test eps1 < 0.005
print("\nTesting auxiliary functions:\n")
print("----------------------------\n")
print("Testing fitslopex\n")
fx=(1.:10.).^-2;
sx=fitslopex(fx);
eps1=abs(sx-1);
@test eps1 < 1e-8
print("Testing lon_lat_fine\n")
lon=0:7; lat=0:7; nf=2;
(lon_f,lat_f) = lon_lat_fine(lon, lat,nf)
lon_c=collect(-0.25:0.5:7.25);
eps1=sum((lon_c-lon_f).^2)
@test eps1 < 1e-20