text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Shell_Pdrop
	Shell_Pdrop()=begin
		new(
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Ideal Pressure Drop",
				:Symbol=>"\\Delta P_{ideal}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Cross Flow Pressure Drop",
				:Symbol=>"\\Delta P_{CrossFlow}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"End Zones Pressure Drop",
				:Symbol=>"\\Delta P_{EndZones}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Window Pressure Drop",
				:Symbol=>"\\Delta P_{Window}",
				:Default=>0.01,
				:Lower=>1e-10,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Pressure Drop",
				:Symbol=>"\\Delta P_{Nozzle\\_In}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle Pressure Drop",
				:Symbol=>"\\Delta P_{Nozzle\\_Out}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Total Pressure Drop",
				:Symbol=>"\\Delta P_{Total}",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa"
			)),
			fricfactor(Dict{Symbol,Any}(
				:Brief=>"Friction Factor",
				:Symbol=>"f_i",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle Velocity",
				:Symbol=>"V_{Nozzle\\_In}",
				:Default=>1,
				:Upper=>1e5,
				:Lower=>0
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle Velocity",
				:Symbol=>"V_{Nozzle\\_Out}",
				:Default=>1,
				:Upper=>1e5,
				:Lower=>0
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outlet Nozzle rho-V^2",
				:Default=>1,
				:Upper=>1e6,
				:Unit=>"kg/s^2/m"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inlet Nozzle rho-V^2",
				:Default=>1,
				:Upper=>1e6,
				:Unit=>"kg/s^2/m"
			)),
			[
				:(Total = CrossFlow+ EndZones + InletNozzle + OutletNozzle + Window),
			],
			[
				"Shell Side Total Pressure Drop",
			],
			[:Ideal,:CrossFlow,:EndZones,:Window,:InletNozzle,:OutletNozzle,:Total,:FricFactor,:Vnozzle_in,:Vnozzle_out,:RVsquare_out,:RVsquare_in,]
		)
	end
	Ideal::press_delta
	CrossFlow::press_delta
	EndZones::press_delta
	Window::press_delta
	InletNozzle::press_delta
	OutletNozzle::press_delta
	Total::press_delta
	FricFactor::fricfactor
	Vnozzle_in::velocity
	Vnozzle_out::velocity
	RVsquare_out::positive
	RVsquare_in::positive
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Shell_Pdrop
function setEquationFlow(in::Shell_Pdrop)
	addEquation(1)
end
function atributes(in::Shell_Pdrop,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Pressure drop and velocities in the shell side section of a shell and tube heat exchanger."
	drive!(fields,_)
	return fields
end
Shell_Pdrop(_::Dict{Symbol,Any})=begin
	newModel=Shell_Pdrop()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
mutable struct ApplicationState
    resolution::NTuple{2,Int}
    scale::NTuple{2,Int}
    position::Point{2,Int}
    noise::Matrix{Float64}
    haschanged::Bool
end

shader_file(filename) = joinpath(@__DIR__, "shaders", filename)

ApplicationState(resolution, scale, position) = ApplicationState(resolution, scale, position, zeros(resolution...), false)

haschanged(app::ApplicationState) = app.haschanged

function ApplicationState(position = (1920 / 2, 1080 / 2))
    app = ApplicationState((512, 512), (4, 4), position)
    update!(app)
    app
end

struct VulkanConfig{S}
    instance::Created{Instance,InstanceCreateInfo}
    device::Created{Device,DeviceCreateInfo}
    surface::Created{SurfaceKHR,S}
    disp::QueueDispatch
    command_pools::ThreadedCommandPool
end

struct Application{WM<:AbstractWindowManager}
    state::ApplicationState
    wm::WM
    gui::GUIManager{WM}
    vk::VulkanConfig
    gr::GUIRenderer
end

main_window(wm::WindowManager) = first(values(wm.impl.windows))

function add_widget!(app::Application, w::Widget, widget_dependencies::WidgetDependencies, callbacks=WindowCallbacks())
    insert!(app.gui.callbacks, w, callbacks)
    add_widget!(app.gr, w, widget_dependencies)
end

function update_vertex_buffer(widget::Widget, app::Application, data)
    upload_data(app.gr.render_set.render_infos[widget].shader_dependencies.vertex_buffer, data)
end

function update_index_buffer(widget::Widget, app::Application, data)
    upload_data(app.gr.render_set.render_infos[widget].shader_dependencies.index_buffer, data)
end

function add_perlin_image!(app::Application)
    add_widget!(
        app,
        ImageWidget(app.state.position, (512, 512), (1.0, 1.0)),
        WidgetDependencies(ShaderSpecification.(shader_file.(["texture_2d.vert", "texture_2d.frag"]), GLSL)),
        WidgetCallbacks(
            on_drag = (src_w::ImageWidget, src_ed::EventDetails, _, ed::EventDetails) -> begin
                Δloc = Point(ed.location) - Point(src_ed.location)
                new_img = @set src_w.center = src_w.center + Δloc
                update_vertex_buffer(src_w, app, new_img)
            end,
            on_drop = (src_w::ImageWidget, src_ed::EventDetails, _, dst_ed::EventDetails) -> begin
                Δloc = Point(dst_ed.location) - Point(src_ed.location)
                app.state.position = src_w.center + Δloc
                new_img = @set src_w.center = src_w.center + Δloc
                add_perlin_image!(app)
            end,
        )
    )
end

function Base.run(app::Application, mode::ExecutionMode = Synchronous(); render = true)
    if render
        device = app.vk.device
        fstate = FrameState(setup_swapchain(handle(device), handle(app.vk.surface))...)
        add_perlin_image!(app)
        image_info = ImageCreateInfo(
            IMAGE_TYPE_2D,
            FORMAT_R16G16B16A16_SFLOAT,
            Extent3D(app.state.resolution..., 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,
        )
        image = unwrap(create_image(device, image_info))
        memory = DeviceMemory(image, MEMORY_PROPERTY_DEVICE_LOCAL_BIT)
        view_info = ImageViewCreateInfo(
            image,
            IMAGE_VIEW_TYPE_2D,
            image_info.format,
            ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
            ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
        )
        view = unwrap(create_image_view(device, view_info))
        add_widget!(
            gr,
            :perlin,
            app.gui.widgets[:perlin],
            ShaderInfo(
                Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "texture_2d.vert"), FormatGLSL()), DescriptorBinding[]),
                Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "texture_2d.frag"), FormatGLSL()), [DescriptorBinding(DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, 0, 1)]),
            ),
            (
                (resource) -> begin
                    Extent3D(app.state.resolution..., 1) ≠ resource.image.info.extent
                end,
                (resource) -> begin
                    image_info = resource.image.info
                    image_info = @set image_info.extent = Extent3D(app.state.resolution..., 1)
                    image = unwrap(create_image(device, image_info))
                    memory = DeviceMemory(image, MEMORY_PROPERTY_DEVICE_LOCAL_BIT)
                    view_info = resource.view.info
                    view_info = @set view_info.image = image
                    view = unwrap(create_image_view(device, view_info))
                    (SampledImage(GPUResource(image, memory, image_info), GPUResource(view, nothing, view_info), resource.sampler),)
                end,
            ),
            SampledImage(
                GPUResource(image, memory, image_info),
                GPUResource(view, nothing, view_info),
                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,
                    get_physical_device_properties(device.physical_device).limits.max_sampler_anisotropy,
                    false,
                    COMPARE_OP_ALWAYS,
                    0,
                    0,
                    BORDER_COLOR_FLOAT_OPAQUE_BLACK,
                    false,
                ),
            ),
        )
        box1 = Box((1500., 540.), (100., 20.), RGB(0.,0.7,0.5))
        app.gui.widgets[:box1] = box1
        app.gui.callbacks[box1] = WidgetCallbacks(
            on_double_click = (_...) -> update!(app.state)
        )
        update_buffers!(rdr.gpu, device, :box1, box1)
        add_widget!(
            gr,
            :box1,
            box1,
            ShaderInfo(
                Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "box.vert"), FormatGLSL()), DescriptorBinding[]),
                Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "box.frag"), FormatGLSL()), DescriptorBinding[]),
            ),
        )
        # text1 = StaticText(
        #     text = "Hello World!",
        #     origin = (500., 500.)
        # )
        # add_widget!(
        #     gr,
        #     :text1,
        #     text1,
        #     ShaderInfo(
        #         Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "box.vert"), FormatGLSL()), DesciptorBinding[]),
        #         Shader(device, ShaderFile(joinpath(@__DIR__, "shaders", "glyph_minimal_manual.frag"), FormatGLSL()), [DescriptorBinding(DESCRIPTOR_TYPE_STORAGE_BUFFER, 0, 1)]),
        #     ),
        #     create_buffer_resource(device, text1.properties.indexing.curves) do (device, size)
        #         unwrap(create_buffer(device, size, BUFFER_USAGE_STORAGE_BUFFER_BIT, SHARING_MODE_EXCLUSIVE, [0]))
        #     end,
        # )
        transfer_texture!(gr, app.state)
        map_window(main_window(app.wm))
        run(app.gui, mode; on_iter_last = () -> begin
            @timeit to "Main loop" begin
                frame = rstate.frame
                if needs_resource_update(gr)
                    @timeit to "Update resources" begin
                        wait_hasrendered(frame)
                        update_resources(gr)
                    end
                end
                if haschanged(app.state)
                    @timeit to "Transfer texture" transfer_texture!(gr, app.state)
                    app.state.haschanged = false
                end
                @timeit to "Draw next frame" next_frame!(frame, rdr, app)
            end
        end)
        gpu = app.state.gpu
        GC.@preserve gpu rdr rstate device_wait_idle(rdr.device)
    else
        run(app.gui, mode)
    end
end

function on_button_pressed(details::EventDetails)
    x, y = details.location
    click = details.data.button
    state = details.data.state
    buttons_pressed = pressed_buttons(state)
    printed_state = isempty(buttons_pressed) ? "" : "with $(join(string.(buttons_pressed), ", ")) button$(length(buttons_pressed) > 1 ? "s" : "") held"
    @info "$click at $x, $y $printed_state"
end

function Application()
    app_state = ApplicationState()
    connection = Connection()
    win = XCBWindow(
        connection;
        x = 20,
        y = 20,
        width = 1920,
        height = 1080,
        border_width = 50,
        window_title = "Givre",
        icon_title = "Givre",
        attributes = [XCB.XCB_CW_BACK_PIXEL],
        values = [0],
        map = false
    )
    wm = XWindowManager(connection, [win])
    wwm = WindowManager(wm)
    gm = GUIManager(wwm)
    queue_config = dictionary([QUEUE_COMPUTE_BIT | QUEUE_GRAPHICS_BIT => 1])
    instance, device = init(;
        instance_extensions = ["VK_KHR_surface", "VK_KHR_xcb_surface"],
        device_extensions = ["VK_KHR_swapchain", "VK_KHR_synchronization2"],
        enabled_features = PhysicalDeviceFeatures(:sampler_anisotropy),
        queue_config = queue_config,
    )
    surface_ci = XcbSurfaceCreateInfoKHR(connection.h, main_window(wwm).id)
    surface = Created(unwrap(create_xcb_surface_khr(instance, surface_ci)), surface_ci)
    queue_dispatch = QueueDispatch(device, queue_config; surface)
    vk = VulkanConfig(instance, device, surface, queue_dispatch, ThreadedCommandPool(device, queue_dispatch, dictionary([1 => QUEUE_COMPUTE_BIT | QUEUE_GRAPHICS_BIT])))
    app = Application(app_state, wwm, gm, vk, GUIRenderer(device, queue_dispatch))

    _update! = app_state -> begin
        update!(app_state)
        app_state.haschanged = true
    end

    key_mappings = Dict{KeyCombination,Any}(
        key"ctrl+q" => (ev, _) -> throw(CloseWindow(ev.win, "Received closing request from user input")),
        key"s" => (_, _) -> begin
            app_state.scale = app_state.scale .+ 2
            _update!(app_state)
        end,
        key"j" => (_, _) -> begin
            app_state.resolution = app_state.resolution .+ 50
            _update!(app_state)
        end,
        key"k" => (_, _) -> begin
            app_state.resolution = app_state.resolution .- 50
            _update!(app_state)
        end,
    )

    function on_key_pressed(details::EventDetails)
        @unpack win, data = details
        @info keystroke_info(wm.keymap, details)
        @unpack key, modifiers = data
        kc = KeyCombination(key, modifiers)
        if haskey(key_mappings, kc)
            key_mappings[kc](details, app_state)
        end
    end

    set_callbacks!(wm, win, WindowCallbacks(;
        on_key_pressed,
        on_key_released = identity,
        on_pointer_move = identity,
        on_mouse_button_pressed = identity,
        on_mouse_button_released = identity,
    ))

    app
end

ext: jl


text:
# Python dependencies.
# Pairs are rc_key => library
PYDEPS = Dict(
    :scipy => "scipy.spatial",
    :pymatgen => "pymatgen.io.cif"
)

# Attempts to load a python dependency.
# Calling functions which require these dependencies without successfully loading them will throw exceptions.
function load_pydep(pydep::String)
    try
        return pyimport(pydep)
    catch
        return nothing
    end
end

# Checks for ability to load a python dependency.
function check_pydep(pydep::Pair)
    if isnothing(load_pydep(pydep[2]))
        @warn "Error loading $(pydep[1]). Some functionaltiy may be missing."
        return missing
    end
end

# Called at module initialization to spin up deps
function init_pydeps()
    for pydep ∈ PYDEPS
        rc[pydep[1]] = load_pydep(pydep[2])
    end
end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
import MapImages
import MapImages: grid_extent, MapImage
import Sario
import CSV
import PyPlot;
plt = PyPlot;
using ImageFiltering
import Images: imresize
using DataFramesMeta
using PyCall;
gps = pyimport("apertools.gps");
include("./map_plots.jl")

yrsymbols = Symbol.(collect(2009:2017))
rename_years(df) = rename(df, Dict(yr => Symbol(:y, yr) for yr in yrsymbols) )

barrel2gal = 42.0
gal2barrel = 1.0 / barrel2gal

dfmt = dateformat"y-m-d H:M:S.s"
readdf(fname) = coalesce.(CSV.read(fname, dateformat = dfmt), 0)

demrsc = Sario.load("dem_all.rsc")
water_prod = readdf("water_production.csv")
oil_prod = readdf("oil_production.csv")
gas_prod = rename_years(readdf("gas_production.csv"));
gas_inj = readdf("gas_injection.csv");
water_inj = readdf("water_injection.csv")

#### Clean known sketchy values
gas_prod = @where(gas_prod, (:y2015 .< 10000000) .| (:y2017 .< 40000000))

waterinjbins = bin_vals(demrsc, water_inj)

waterprodbins_small = coarse_bin_vals(demrsc, water_prod, digits = 1)
waterprodbins = imresize(waterprodbins_small, size(demrsc)...)
waterprodbins .*= (sum(bin_vals(demrsc, water_prod)) / sum(waterprodbins))

oilprodbins_small = coarse_bin_vals(demrsc, oil_prod, digits = 1)
oilprodbins = imresize(oilprodbins_small, size(demrsc)...)
oilprodbins .*= (sum(bin_vals(demrsc, oil_prod)) / sum(oilprodbins))

netbins = waterinjbins .- waterprodbins .- oilprodbins

eqs = readdf("texnet_events_latest.csv")
eq_bins =
    bin_vals(demrsc, eqs, loncol = :Longitude, latcol = :Latitude, valcol = :Magnitude)

eqs[:, :MagnitudeExp] = exp.(eqs[:, :Magnitude])
eqs[:, :Date] = eqs[:, Symbol("Origin Date")]
eqs19 = @where(eqs, :Date .< Date(2020,1,1))
eq_bins_logexp = log.(bin_vals(demrsc, eqs19, loncol = :Longitude, latcol = :Latitude, valcol = :MagnitudeExp, smoothing=2))
eq_bins_logexp[eq_bins_logexp .== -Inf] .= NaN;

fracs_all = readdf("FracDetails_Dates_TotalVolume.csv")
# frac_bins = bin_vals(demrsc, fracs, loncol=:mid_longitude, latcol=:mid_latitude, valcol=:total_volume_gals)
fracs = readdf("Frac_2015_simple.csv")
frac_bins = bin_vals(demrsc, fracs, loncol = :lon, latcol = :lat, valcol = :volume)

#im78 = MapImage("../images_path78/velocities_2018_linear_max800.h5", "velos_shifted/1");
binsblur = imfilter(netbins, Kernel.gaussian(2))
eqsblur = imfilter(eq_bins, Kernel.gaussian(1))
fracsblur = imfilter(frac_bins, Kernel.gaussian(2))
#
fig, axes = plt.subplots(2, 2, sharex = true, sharey = true)
#axim = axes[1, 1].imshow(im78, cmap="seismic", vmin=-20, vmax=20, extent=grid_extent(im78)); fig.colorbar(axim, ax=axes[1,1])
axim = axes[1, 2].imshow(
    binsblur,
    vmax = 1e6,
    vmin = -1e6,
    cmap = "seismic",
    extent = grid_extent(demrsc),
);
fig.colorbar(axim, ax = axes[1, 2]);
axim = axes[2, 1].imshow(fracsblur, vmax = 2e7, cmap = "Reds", extent = grid_extent(demrsc));
fig.colorbar(axim, ax = axes[2, 1]);
axim = axes[2, 2].imshow(eqsblur, vmax = 5, cmap = "Reds", extent = grid_extent(demrsc));
fig.colorbar(axim, ax = axes[2, 2]);
#

axes[1, 1].set_title("InSAR")
axes[1, 2].set_title("Water Inj. - Water Prod. - Oil Prod")
axes[2, 1].set_title("Fracs")
axes[2, 2].set_title("Eqs")
# axes[1, 2].imshow(netbins, vmax=1e6, vmin=-1e6, cmap="seismic")
#
oil_prod17 = @where(oil_prod, cols(Symbol(2017)) .> 0);
water_inj17 = @where(water_inj, cols(Symbol(2017)) .> 0);
water_prod17 = @where(water_prod, cols(Symbol(2017)) .> 0);

# @time frac_simple = DataFrame(load("FracFocus_OneLine_Final.xlsx", "simple_heel"));
# frac2016 = @where(frac_simple, cols(st) .> Date(2016,1,1), cols(st) .< (Dates.Year(1) + Date(2016,1,1)));


text:
using Test, Random, LSHFunctions, LinearAlgebra, SparseArrays

include(joinpath("..", "utils.jl"))

#==================
Tests
==================#
@testset "SignALSH tests" begin
    Random.seed!(RANDOM_SEED)

    @testset "Construct SignALSH" begin
        hashfn = SignALSH(; maxnorm=1)
        @test n_hashes(hashfn) == 1
        @test isa(hashfn, SignALSH{Float32})
        @test isa(hashfn, LSHFunctions.AsymmetricLSHFunction)
        @test hashtype(hashfn) == Bool
        @test similarity(hashfn) == inner_prod

        hashfn = SignALSH(32; maxnorm=1)
        @test n_hashes(hashfn) == 32

        hashfn = SignALSH(; dtype=Float64, maxnorm=1)
        @test isa(hashfn, SignALSH{Float64})

        # n_hashes must be positive
        @test_throws ErrorException SignALSH(-1; maxnorm=1)
        @test_throws ErrorException SignALSH( 0; maxnorm=1)

        # maxnorm must be specified and non-negative
        @test_throws ErrorException SignALSH()
        @test_throws ErrorException SignALSH(; maxnorm=-1)

        # m must be positive
        @test_throws ErrorException SignALSH(; m=-1)
        @test_throws ErrorException SignALSH(; m=0)
    end

    @test_skip @testset "Can hash inputs correctly with SignALSH" begin
        input_length = 5
        n_hashes = 8

        hashfn = SignALSH(n_hashes; m=3)
        simhash = SimHash([hashfn.coeff_A hashfn.coeff_B])

        # By default, hashfn should be using Float32 for its coefficients
        @test isa(hashfn, SignALSH{Float32})

        X = rand(Float32, input_length, 128)
        ihashes = index_hash(hashfn, X)
        qhashes = query_hash(hashfn, X)

        @test eltype(ihashes) == eltype(qhashes)
        @test hashtype(hashfn) == Bool

        # 1. Compute the indexing hashes manually
        norms = map(norm, eachcol(X))
        maxnorm = maximum(norms)
        norms .*= 1 / maxnorm
        Y = X ./ maxnorm

        @test all(map(norm,eachcol(Y)) .≤ 1)
        @test all(0 .≤ norms .≤ 1)

        norm_powers = -[norms.^2 norms.^4 norms.^8] .+ 1/2
        Px = [Y; norm_powers']

        @test simhash(Px) == ihashes

        # 2. Compute the querying hashes manually
        Qx = [X ./ norms'; zeros(1,128); zeros(1,128); zeros(1,128)]

        @test simhash(Qx) == qhashes
    end

    @testset "SignALSH can't hash inputs of norm > maxnorm" begin
        hashfn = SignALSH(; maxnorm=0)
        @test_throws ErrorException index_hash(hashfn, rand(4))
        @test_throws ErrorException query_hash(hashfn, rand(4))

        # Should have no issue if norm(x) == maxnorm
        @test index_hash(hashfn, zeros(4)) |> length == 1
        @test query_hash(hashfn, zeros(4)) |> length == 1
    end

    @testset "SignALSH generates collisions for large inner products" begin
        input_length = 5; n_hashes = 128;

        # Compare a random vector x against four other vectors:
        # a) 10 * x
        # b) x
        # c) A vector of all zeros
        # d) -x
        x = randn(input_length)
        x2, x3, x4 = 10*x, zero(x), -x

        maxnorm = (x, x2, x3, x4) .|> norm |> maximum
        hashfn = SignALSH(n_hashes; maxnorm=maxnorm)

        x_query_hashes = query_hash(hashfn, x)

        dataset = [x2 x x3 x4]
        p_hashes = index_hash(hashfn, dataset)

        # Each collection of hashes should be different from one another
        @test let result = true
            for (ii,jj) in Iterators.product(1:4, 1:4)
	            if ii != jj && p_hashes[:,ii] == p_hashes[:,jj]
		            result = false
		            break
	            end
            end
            result
        end

        # The number of collisions should be highest for x and 2*x, second-highest
        # for x and x, second-lowest for x and zeros, and lowest for x and -x
        n_collisions = [sum(x_query_hashes .== p) for p in eachcol(p_hashes)]
        @test n_collisions[1] > n_collisions[2] >
              n_collisions[3] > n_collisions[4]
    end

    @testset "Can hash sparse arrays" begin
        input_size = 100
        n_inputs = 150
        n_hashes = 2

        hashfn = SignALSH(n_hashes; maxnorm=4*input_size)
        x = sprandn(input_size, n_inputs, 0.2)

        # Mostly just need to test that the following lines don't crash
        ih = index_hash(hashfn, x)
        qh = query_hash(hashfn, x)

        @test size(ih) == (n_hashes, n_inputs)
        @test size(qh) == (n_hashes, n_inputs)
    end

    @testset "Can hash matrix adjoints" begin
        input_size = 100
        n_inputs = 150
        n_hashes = 2
        hashfn = SignALSH(n_hashes; maxnorm=4*input_size)

        ## Test 1: regular matrix adjoint
        x = randn(n_inputs, input_size)'
        @test index_hash(hashfn, x) == index_hash(hashfn, copy(x))
        @test query_hash(hashfn, x) == query_hash(hashfn, copy(x))

        ## Test 2: sparse matrix adjoint
        x = sprandn(n_inputs, input_size, 0.2)'
        @test index_hash(hashfn, x) == index_hash(hashfn, copy(x))
        @test query_hash(hashfn, x) == query_hash(hashfn, copy(x))
    end

    @testset "Hash inputs of different sizes" begin
        n_hashes = 42
        hashfn = SignALSH(n_hashes; maxnorm=100)

        @test size(hashfn.coeff_A) == (n_hashes, 0)

        index_hash(hashfn, rand(10))
        @test size(hashfn.coeff_A) == (n_hashes, 10)

        query_hash(hashfn, rand(20))
        @test size(hashfn.coeff_A) == (n_hashes, 20)

        index_hash(hashfn, rand(5))
        @test size(hashfn.coeff_A) == (n_hashes, 20)

        query_hash(hashfn, rand(15))
        @test size(hashfn.coeff_A) == (n_hashes, 20)

        index_hash(hashfn, rand(20))
        @test size(hashfn.coeff_A) == (n_hashes, 20)

        query_hash(hashfn, rand(20))
        @test size(hashfn.coeff_A) == (n_hashes, 20)
    end

    @testset "Hash inputs of different sizes with resize_pow2 = true" begin
        n_hashes = 25
        hashfn = SignALSH(n_hashes; maxnorm=100, resize_pow2=true)

        @test size(hashfn.coeff_A) == (n_hashes, 0)

        index_hash(hashfn, rand(4))
        @test size(hashfn.coeff_A) == (n_hashes, 4)

        query_hash(hashfn, rand(8))
        @test size(hashfn.coeff_A) == (n_hashes, 8)

        index_hash(hashfn, rand(9))
        @test size(hashfn.coeff_A) == (n_hashes, 16)

        query_hash(hashfn, rand(17))
        @test size(hashfn.coeff_A) == (n_hashes, 32)

        index_hash(hashfn, rand(15))
        @test size(hashfn.coeff_A) == (n_hashes, 32)

        query_hash(hashfn, rand(31))
        @test size(hashfn.coeff_A) == (n_hashes, 32)

        index_hash(hashfn, rand(32))
        @test size(hashfn.coeff_A) == (n_hashes, 32)

        query_hash(hashfn, rand(32))
        @test size(hashfn.coeff_A) == (n_hashes, 32)
    end
end


text:
using Discourse

function get_user(username::String)
    """
        Get User information for a specific user.
        ARGS:
            username::String -> username to return
        RETURNS:
            user_info::Dictionary
    """
    return _get("/users/$(username).json")["user"]

end


text:
# Vertex traits w.r.t whether size changes propagates
abstract type Direction end
"""
    Input

Represents the input direction, i.e coming from the output of another vertex.
"""
struct Input <: Direction end
"""
    Output

Represents the output direction, i.e coming from the input of another vertex.
"""
struct Output <: Direction end
"""
    Both

Represents both directions (`Input` and `Output`).
"""
struct Both <: Direction end

"""
    opposite(d::Direction)

Return the opposite direction of `d`.
"""
opposite(::Input) = Output()
opposite(::Output) = Input()
opposite(b::Both) = b

"""
    neighbours(d::Direction, v)

Return vertices connected to `v` in direction `d`.
"""
neighbours(::Input, v) = inputs(v)
neighbours(::Output, v) = outputs(v)
neighbours(::Both, v) = vcat(inputs(v), outputs(v))

"""
    ΔSizeGraph

Represents the information on how a size change will propagate as a `MetaDiGraph`.

Each vertex `i` represents a unique `AbstractVertex vi`.
Each edge `e` represents propagation of a size change between vertices `e.src` and `e.dst`.
Edge weights denoted by the symbol `:size` represents the size of the output sent between the vertices.

For the `AbstractVertex vi` associated with vertex `i` in the graph `g`, the following holds `g[i, :vertex] == vi` and `g[vi,:vertex] == i`)

For an edge `e` in graph `g`, the following holds:

If `get_prop(g, e, :direction)` is of type `Output` this means that `Δnout` of `e.dst` is called after processing `e.src`.

If `get_prop(g, e, :direction)` is of type `Input` this means that `Δnin` of `e.dst` is called after processing `e.src`.
"""
function ΔSizeGraph()
    g = MetaDiGraph(0, :size, -1)
    set_indexing_prop!(g, :vertex)
    return g
end

"""
    ΔninSizeGraph(v)

Return a `ΔSizeGraph` for the case when nin of `v` is changed, i.e when Δnin(v, Δ) is called.
"""
ΔninSizeGraph(v) = ΔSizeGraph(Input(), v)

"""
    ΔnoutSizeGraph(v)

Return a `ΔSizeGraph` for the case when nout of `v` is changed, i.e when Δnout(v, Δ) is called.
"""
ΔnoutSizeGraph(v) = ΔSizeGraph(Output(), v)

function ΔSizeGraph(d::Direction, v)
    g = ΔSizeGraph()
    set_prop!(g, :start, v => d)
    verts = all_in_Δsize_graph(trait(v), d, v, (g, v))
    return g
end

function all_in_Δsize_graph(v::AbstractVertex, d::Direction, (g,from)::Tuple{MetaDiGraph, AbstractVertex})
    has_edge(g, vertexind!(g, v), vertexind!(g, from)) && return g
    add_edge!(g, from, v, d) && return g
    all_in_Δsize_graph(trait(v), d, v, (g, v))
end

function LightGraphs.add_edge!(g::MetaDiGraph, src::AbstractVertex, dst::AbstractVertex, d::Direction)
    srcind = vertexind!(g, src)
    dstind = vertexind!(g, dst)
    visited(g, srcind, dstind, d) && return true
    add_edge!(g, srcind, dstind, Dict(:direction => d, g.weightfield => edgesize(d, src, dst)))
    return false
end

function visited(g, srcind, dstind, d)
    has_edge(g, srcind, dstind) && return true
    return visited_out(d, g, dstind)
end

function visited_out(::Output, g, dstind)
    for e in filter_edges(g, :direction, Output())
        e.dst == dstind && return true
    end
    return false
end
visited_out(::Input, g, dstind) = false

edgesize(::Input, src, dst) = nout(src)
edgesize(::Output, src, dst) = nout(dst)

vertexproplist(g::MetaDiGraph, prop::Symbol) = map(p -> p[:vertex], props.([g], vertices(g)))

function vertexind!(g::MetaDiGraph, v::AbstractVertex,)
    ind = indexin([v], vertexproplist(g, :vertex))[]
    if nothing == ind
        add_vertex!(g, :vertex, v)
    end
    return g[v, :vertex]
end

"""
    SizeDiGraph(cg::CompGraph)

Return `cg` as a `MetaDiGraph g`.

Each vertex `i` represents a unique `AbstractVertex vi`.

For the `AbstractVertex vi` associated with vertex `i` in the graph `g`, the following holds `g[i, :vertex] == vi` and `g[vi,:vertex] == i`)

Each edge `e` represents output from `e.src` which is input to `e.dst`.
Edge weights denoted by the symbol `:size` represents the size of the output sent between the vertices.

Note that the order of the input edges to a vertex matters (in case there are more than one) and is not encoded in `g`.
Instead, use `indexin(vi, inputs(v))` to find the index (or indices if input multiple times) of `vi` in `inputs(v)`.
"""
SizeDiGraph(cg::CompGraph) = SizeDiGraph(mapfoldl(v -> flatten(v), (vs1, vs2) -> unique(vcat(vs1, vs2)), cg.outputs))

"""
    SizeDiGraph(v::AbstractVertex)

Return a SizeDiGraph of all parents of v
"""
SizeDiGraph(v::AbstractVertex)= SizeDiGraph(flatten(v))

"""
    SizeDiGraph(vertices::AbstractArray{AbstractVertex,1})

Return a SizeDiGraph of all given vertices
"""
function SizeDiGraph(vertices::AbstractArray{<:AbstractVertex,1})
    g = MetaDiGraph(0,:size, 0)
    set_indexing_prop!(g, :vertex)
    add_vertices!(g, length(vertices))
    for (ind, v) in enumerate(vertices)
        set_prop!(g, ind, :vertex, v)
        for in_ind in indexin(inputs(v), vertices)
            add_edge!(g, in_ind, ind, :size, nout(vertices[in_ind]))
        end
    end
    return g
end

"""
    fullgraph(v::AbstractVertex)

Return a `SizeDiGraph` of all vertices in the same graph (or connected component) as `v`
"""
fullgraph(v::AbstractVertex) = SizeDiGraph(all_in_graph(v))

"""
    all_in_graph(v::AbstractVertex)

Return an array of vertices in the same graph (or connected component) as `v`
"""
function all_in_graph(v::AbstractVertex, visited = AbstractVertex[])
    v in visited && return visited
    push!(visited, v)
    foreach(vi -> all_in_graph(vi, visited), inputs(v))
    foreach(vo -> all_in_graph(vo, visited), outputs(v))
    return visited
end

"""
    all_in_Δsize_graph(v::AbstractVertex, d::Direction)

Return an array of vertices which will be affected if `v` changes size in direction `d`.
"""
function all_in_Δsize_graph(v::AbstractVertex, d::Direction, visited=[])
    (v, d) in visited && return visited
    push!(visited, (v, d))
    all_in_Δsize_graph(trait(v),d, v, visited)
    return unique(map(e -> e[1], visited))
end
function all_in_Δsize_graph(v::AbstractVertex, d::Both, visited=[])
    all_in_Δsize_graph(v, Input(), visited)
    foreach(vout -> all_in_Δsize_graph(vout, Input(), visited), outputs(v))
    return unique(map(e -> e[1], visited))
end

all_in_Δsize_graph(t::DecoratingTrait, d, v, visited) = all_in_Δsize_graph(base(t), d, v, visited)
function all_in_Δsize_graph(::Immutable, ::Direction, v, visited) end
all_in_Δsize_graph(::SizeAbsorb, d, v, visited) = foreach(vn -> all_in_Δsize_graph(vn, opposite(d), visited), neighbours(d, v))
function all_in_Δsize_graph(::SizeTransparent, d, v, visited)
    foreach(vin -> all_in_Δsize_graph(vin, Output(), visited), inputs(v))
    foreach(vout -> all_in_Δsize_graph(vout, Input(), visited), outputs(v))
end


text:
"""
LX_TOKENS

List of names of latex tokens. (See markdown tokens)
"""
const LX_TOKENS = (:LX_BRACE_OPEN, :LX_BRACE_CLOSE, :LX_COMMAND)

"""
$(TYPEDEF)

Structure to keep track of the definition of a latex command declared via a
`\\newcommand{\\name}[narg]{def}` or of an environment via
`\\newenvironment{name}[narg]{pre}{post}`.
The parametric type depends on the definition type, for a command it will be
a SubString (<: AbstractString), for an environment it will be a Pair of
SubString (corresponding to pre-post).
"""
struct LxDef{T}
    name::String
    narg::Int
    def ::T
    # location of the definition
    from::Int
    to  ::Int
end
# if offset unspecified, start from basically -∞ (configs etc)
function LxDef(name::AS, narg::Int, def)
    o = FD_ENV[:OFFSET_LXDEFS] += 5  # precise offset doesn't matter
    LxDef(string(name), narg, def, o, o + 3) # just forward a bit
end

from(lxd::LxDef) = lxd.from
to(lxd::LxDef)   = lxd.to

"""
pastdef(λ)

Convenience function to return a copy of a definition indicated as having
already been earlier in the context i.e.: earlier than any other definition
appearing in the current chunk.
"""
pastdef(λ::LxDef) = LxDef(λ.name, λ.narg, λ.def)

"""
$(TYPEDEF)

Super type for `LxCom` and `LxEnv`.
"""
abstract type LxObj <: AbstractBlock end

"""
$(TYPEDEF)

A `LxCom` has a similar content as a `Block`, with the addition of the definition and a vector of brace blocks.
"""
struct LxCom <: LxObj
    ss    ::SubString                 # \\command
    lxdef ::Union{Nothing,Ref{LxDef}} # definition of the command
    braces::Vector{OCBlock}           # relevant {...} with the command
end
LxCom(ss, def) = LxCom(ss, def, Vector{OCBlock}())

"""
$TYPEDEF

A `LxEnv` is similar to a `LxCom` but for an environment.

    `\\begin{aaa} ... \\end{aaa}`
    `\\begin{aaa}{opt1}{opt2} ... \\end{aaa}`
"""
struct LxEnv <: LxObj
    ss    ::SubString
    lxdef ::Union{Nothing,Ref{LxDef}}
    braces::Vector{OCBlock}
    ocpair::Pair{Token,Token}
end
LxEnv(ss, def, ocp) = LxCom(ss, def, Vector{OCBlock}(), ocp)

"""
$SIGNATURES

Content of an `LxEnv` block.

    `\\begin{aaa}{opt1} XXX \end{aaa}` --> ` XXX `
"""
function content(lxe::LxEnv)
    s = str(lxe.ss)
    cfrom = nextind(s, to(lxe.ocpair.first))
    if !isempty(lxe.braces)
        cfrom = nextind(s, to(lxe.braces[end]))
    end
    cto = prevind(s, from(lxe.ocpair.second))
    return subs(s, cfrom, cto)
end

"""
$SIGNATURES

For a given `LxObj`, retrieve the definition attached to the corresponding `LxDef` via the
reference.
"""
function getdef(lxo::LxObj)
    isnothing(lxo.lxdef) && return nothing
    return getindex(lxo.lxdef).def
end

"""
$SIGNATURES

For a given `LxObj`, retrieve the name of the object via the reference.
Example: `\\cite` --> `cite` or `\\begin{aaa}` --> `aaa`.
"""
function getname(lxc::LxCom)::String
    if isnothing(lxc.lxdef)
        s = String(lxc.ss)
        j = findfirst('{', s)
        return lxc.ss[2:prevind(s, j)]
    end
    return String(getindex(lxc.lxdef).name)[2:end]
end
getname(lxe::LxEnv) = envname(lxe.ocpair.first)
envname(lxe::LxEnv) = getname(lxe)


text:
for N in [Float64, Rational{Int}, Float32]
    # =====================================
    # Run decompose for different set types
    # =====================================

    function test_directions(set)
        res = σ(N[1, 0], set)[1] == one(N)
        res &= σ(N[0, 1], set)[2] == one(N)
        res &= σ(N[-1, 0], set)[1] == -one(N)
        res &= σ(N[0, -1], set)[2] == -one(N)
        return res
    end
    partition = [[1, 2], [3, 4], [5, 6]]
    b = BallInf(zeros(N, 6), one(N))
    d = decompose(b, partition, HPolygon)
    @test d.array[1] isa HPolygon && test_directions(d.array[1])
    d = decompose(b, partition, Hyperrectangle)
    @test d.array[1] isa Hyperrectangle && test_directions(d.array[1])

    d = decompose(b, partition, LinearMap)
    @test d.array[1] isa LinearMap && test_directions(d.array[1])

    d = decompose(b, [[i] for i in 1:6], Interval)
    @test d.array[1] isa Interval &&
        σ(N[1], d.array[1])[1] == one(N) && σ(N[-1], d.array[1])[1] == -one(N)

    # ===================
    # 1D/3D decomposition
    # ===================

    b1 = Ball1(zeros(N, 7), one(N))
    d = decompose(b1, [[i] for i in 1:7], Hyperrectangle)
    @test length(d.array) == 7
    d = decompose(b1, [[1], 2:3, 4:6, [7]], Hyperrectangle)
    @test length(d.array) == 4

    # ===================
    # template directions
    # ===================

    for dir in [BoxDirections{N}, OctDirections{N}, BoxDiagDirections{N}]
        d = decompose(b, partition, dir)
        @test d isa CartesianProductArray && array(d)[1] isa HPolytope
    end

    # ==========================
    # different options per block
    # ==========================

    block2oa1 = Dict(1 => Hyperrectangle, 2 => HPolygon, 3 => OctDirections{N},
                     4 => Interval)
    block2oa2 = [Hyperrectangle, HPolygon, OctDirections{N}, Interval]
    for block2oa in [block2oa1, block2oa2]  # both Dict and Vector work
        d = decompose(b1, [[1], 2:3, 4:6, [7]], block2oa)
        @test d isa CartesianProductArray && array(d)[1] isa Hyperrectangle &&
              array(d)[2] isa HPolygon && array(d)[3] isa HPolytope &&
              array(d)[4] isa Interval
    end

    # ==================
    # uniform block size
    # ==================

    d = decompose(b, Hyperrectangle; block_size=2)
    @test length(d.array) == 3
    d = decompose(b1, Hyperrectangle; block_size=3)
    @test length(d.array) == 3
end

# tests that do not work with Rational{Int}
for N in [Float64, Float32]
    # =============================
    # Check that Issue #43 is fixed
    # =============================

    X = CartesianProductArray([BallInf(N[0.767292, 0.936613], N(0.1)),
                               BallInf(N[0.734104, 0.87296], N(0.1))])
    A = N[1.92664 1.00674 1.0731 -0.995149;
          -2.05704 3.48059 0.0317863 1.83481;
          0.990993 -1.97754 0.754192 -0.807085;
          -2.43723 0.782825 -3.99255 3.93324]
    Ω0 = CH(X, A * X)
    dec = decompose(Ω0, [1:2, 3:4], HPolygon)
    dec1 = dec.array[1]

    @test dec1.constraints[1].b ≈ N(2.84042586)
    @test dec1.constraints[2].b ≈ N(4.04708832)
    @test dec1.constraints[3].b ≈ N(-0.667292)
    @test dec1.constraints[4].b ≈ N(-0.836613)

    # =====================
    # ε-close approximation
    # =====================

    function test_directions(set)
        res = σ(N[1, 0], set)[1] == one(N)
        res &= σ(N[0, 1], set)[2] == one(N)
        res &= σ(N[-1, 0], set)[1] == -one(N)
        res &= σ(N[0, -1], set)[2] == -one(N)
        return res
    end
    partition = [[1, 2], [3, 4], [5, 6]]
    b = BallInf(zeros(N, 6), one(N))
    d = decompose(b, partition, HPolygon => N(1e-2))
    @test d.array[1] isa HPolygon && test_directions(d.array[1])
    d = decompose(b, partition, N(1e-2))
    @test d.array[1] isa HPolygon && test_directions(d.array[1])
    d = decompose(b, partition, [N(1e-2), HPolygon => N(1e-2), N(1e-2)])
    @test d.array[1] isa HPolygon && test_directions(d.array[1])
end


text:
module TestRollCallEncoding
    using Base.Test
    import RollCallDataIO

    encoding1 = RollCallDataIO.RollCallEncoding()

    @test isa(encoding1, RollCallDataIO.RollCallEncoding)
    @test encoding1.yea === (1, 2, 3)
    @test encoding1.nay === (4, 5, 6)
    @test encoding1.missing === (7, 8, 9)
    @test encoding1.inactive === (0,)

    encoding2 = RollCallDataIO.RollCallEncoding(
        (1, 2),
        (3, 4),
        (5, 6),
        (7, 8, 9),
    )

    @test isa(encoding2, RollCallDataIO.RollCallEncoding)
    @test encoding2.yea === (1, 2)
    @test encoding2.nay === (3, 4)
    @test encoding2.missing === (5, 6)
    @test encoding2.inactive === (7, 8, 9)
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
function reinitModel!(
    OMDATA ,
    data   :: Dict;
    forcing :: Bool,
    thermal :: Bool,
)
    cvtsT = (a,) -> reshape(a, 1, size(a)...)
    SHF     = cvtsT(data["SHF"])
    SHF_QSW = cvtsT(data["SHF_QSW"])
    SFWF    = cvtsT(data["SFWF"])
    TAUX    = cvtsT(data["TAUX"])
    TAUY    = cvtsT(data["TAUY"])
    HMXL    = cvtsT(data["HBLT"])

    if forcing
        @. OMDATA.x2o["SWFLX"]      = - SHF_QSW
        @. OMDATA.x2o["NSWFLX"]     = - (SHF - SHF_QSW)
        @. OMDATA.x2o["VSFLX"]      = SFWF
        @. OMDATA.x2o["TAUX_east"]  = TAUX / 10.0
        @. OMDATA.x2o["TAUY_north"] = TAUY / 10.0

        @. OMDATA.mb.fi.HMXL = HMXL / 100.0
    end

    if thermal
        OMDATA.mb.fi.sv[:TEMP] .= data["TEMP"]
        OMDATA.mb.fi.sv[:SALT] .= data["SALT"]
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function loglikelihoods!(LL::AbstractMatrix, hmm::AbstractHMM{Univariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t, i] = logpdf(hmm.B[i], observations[t])
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::AbstractHMM{Multivariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t, i] = logpdf(hmm.B[i], view(observations, t, :))
    end
end

"""
    loglikelihoods(hmm, observations; robust) -> Matrix

Return the log-likelihood per-state and per-observation.

**Output**
- `Matrix{Float64}`: log-likelihoods matrix (`T x K`).

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000)
LL = likelihoods(hmm, y)
```
"""
function loglikelihoods(hmm::AbstractHMM, observations; logl = nothing, robust = false)
    (logl !== nothing) && deprecate_kwargs("logl")
    T, K = size(observations, 1), size(hmm, 1)
    LL = Matrix{Float64}(undef, T, K)

    loglikelihoods!(LL, hmm, observations)
    if robust
        replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    end
    LL
end


text:
module TestEconometrics

using Dates
using DataFrames
using DataArrays
using TimeData
using Base.Test
using EconDatasets
using TimeData
using Econometrics

#######################################
## create input and output test data ##
#######################################

## create test data
##-----------------

## Econometrics.createTestData()

## build R results
##----------------

## rscriptPath = joinpath(Pkg.dir("Econometrics"), "test/data/")

## comment the following line if you don't have docker
#run(`docker run -t --rm -v $rscriptPath:/home/docker/ juliafinmetrix/rfinm_deb bash R CMD BATCH --no-save --no-restore r_results.R`)

## build MATLAB results
##---------------------

# add MATLAB BATCH job

###############
## run tests ##
###############

my_tests = [
            ## distributions
            "t_loc_scale_test.jl",
            "nchisq_test.jl",

            ## models
            "model_test.jl",
            "normiid_test.jl",
            "tlsiid_test.jl",
            "garch_test.jl"
            ## "cir.jl",

            ## miscellaneous
            ## "bsOptions_test.jl",
            ## "returns.jl"
            ]

println("Running tests:")

for my_test in my_tests
    println("\n * $(my_test)")
    println("---------------------------------------")
    println("---------------------------------------")
    include(my_test)
end

end


text:
using HighestDensityRegions
using HighestDensityRegions:
    getindex_or_zero, interpolation_coefficient_and_index, interpolate_in
using Test

using Distributions
using DocStringExtensions

@testset "utilities" begin
    # getindex
    v = 1:3
    @test getindex_or_zero(v, 3) ≡ 3
    @test getindex_or_zero(v, 0) ≡ 0

    # interpolation index: arbitrary cumprobs
    c = [0.1, 0.4, 1.0]
    @test interpolation_coefficient_and_index(c, 0.1) == (1.0, 1)
    @test interpolation_coefficient_and_index(c, 0.2) == (1/3, 2)
    for _ in 1:20
        α = rand()
        γ, i = interpolation_coefficient_and_index(c, α)
        @test interpolate_in(c, γ, i) ≈ α
    end
    for (i, α) in enumerate(c)
        @test interpolation_coefficient_and_index(c, α) == (1.0, i)
    end

    # interpolation index: integer cumprobs
    N = 8
    cN = Vector((1:N) ./ N)
    for (i, α) in enumerate(cN)
        @test interpolation_coefficient_and_index(cN, α) == (1.0, i)
    end
end

"""
    $SIGNATURES

Calculate densities and counts using the cumulative density function for a
univariate distribution.
"""
function analytical_densities_and_counts(distribution, breaks)
    @assert issorted(breaks)
    counts = diff(cdf.(Ref(distribution), breaks))
    densities = counts ./ diff(breaks)
    densities, counts
end

"""
    $SIGNATURES

Calculate the empirical coverage probability of a (HDR) density threshold.
"""
coverage(threshold, densities, counts) =
    sum(counts[densities .≥ threshold]) / sum(counts)

"""
    $SIGNATURES

Test empirical coverage of HDR thresholds. `atol` is used for comparison, use
dense breaks for more precision.
"""
function test_coverage(αs, thresholds, densities, counts, atol = 0.01)
    for (i, α) in enumerate(αs)
        @test coverage(thresholds[i], densities, counts) ≈ α atol = atol
    end
end

@testset "normal thresholds with evenly spaced breaks" begin
    N = 1000
    breaks = range(-4; stop = 4, length = N)
    densities, counts =
        analytical_densities_and_counts(Normal(0, 1), breaks)
    αs = [0.25, 0.5, 0.75, 0.9]
    thresholds = hdr_thresholds(αs, densities, counts)
    test_coverage(αs, thresholds, densities, counts)
    thresholds_even = hdr_thresholds(αs, densities)
    @test thresholds_even ≈ thresholds
end

@testset "gamma thresholds with log breaks" begin
    breaks = exp.(range(0; stop = 2, length = 1000)) .- 1
    densities, counts = analytical_densities_and_counts(Gamma(1, 1), breaks)
    αs = [0.25, 0.5, 0.75, 0.9]
    thresholds = hdr_thresholds(αs, densities, counts)
    test_coverage(αs, thresholds, densities, counts)
end


text:
using .NumericalFluxes:
    CentralNumericalFluxHigherOrder, CentralNumericalFluxDivergence

struct DGModel{BL, G, NFND, NFD, GNF, AS, DS, HDS, D, DD, MD}
    balancelaw::BL
    grid::G
    numerical_flux_first_order::NFND
    numerical_flux_second_order::NFD
    numerical_flux_gradient::GNF
    state_auxiliary::AS
    state_gradient_flux::DS
    states_higher_order::HDS
    direction::D
    diffusion_direction::DD
    modeldata::MD
end
function DGModel(
    balancelaw,
    grid,
    numerical_flux_first_order,
    numerical_flux_second_order,
    numerical_flux_gradient;
    state_auxiliary = create_auxiliary_state(balancelaw, grid),
    state_gradient_flux = create_gradient_state(balancelaw, grid),
    states_higher_order = create_higher_order_states(balancelaw, grid),
    direction = EveryDirection(),
    diffusion_direction = direction,
    modeldata = nothing,
)
    DGModel(
        balancelaw,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        state_auxiliary,
        state_gradient_flux,
        states_higher_order,
        direction,
        diffusion_direction,
        modeldata,
    )
end

function (dg::DGModel)(
    tendency,
    state_conservative,
    ::Nothing,
    t;
    increment = false,
)

    balance_law = dg.balancelaw
    device = typeof(state_conservative.data) <: Array ? CPU() : CUDA()

    grid = dg.grid
    topology = grid.topology

    dim = dimensionality(grid)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dim == 2 ? 1 : Nq
    Nfp = Nq * Nqk
    nrealelem = length(topology.realelems)

    state_gradient_flux = dg.state_gradient_flux
    Qhypervisc_grad, Qhypervisc_div = dg.states_higher_order
    state_auxiliary = dg.state_auxiliary

    FT = eltype(state_conservative)
    num_state_gradient_flux = number_state_gradient_flux(balance_law, FT)
    nhyperviscstate = num_hyperdiffusive(balance_law, FT)

    Np = dofs_per_element(grid)

    workgroups_volume = (Nq, Nq, Nqk)
    ndrange_volume = (nrealelem * Nq, Nq, Nqk)
    workgroups_surface = Nfp
    ndrange_interior_surface = Nfp * length(grid.interiorelems)
    ndrange_exterior_surface = Nfp * length(grid.exteriorelems)

    communicate =
        !(isstacked(topology) && typeof(dg.direction) <: VerticalDirection)

    update_auxiliary_state!(
        dg,
        balance_law,
        state_conservative,
        t,
        dg.grid.topology.realelems,
    )

    if nhyperviscstate > 0
        hypervisc_indexmap = create_hypervisc_indexmap(balance_law)
    else
        hypervisc_indexmap = nothing
    end

    exchange_state_conservative = NoneEvent()
    exchange_state_gradient_flux = NoneEvent()
    exchange_Qhypervisc_grad = NoneEvent()
    exchange_Qhypervisc_div = NoneEvent()

    comp_stream = Event(device)

    ########################
    # Gradient Computation #
    ########################
    if communicate
        exchange_state_conservative = MPIStateArrays.begin_ghost_exchange!(
            state_conservative;
            dependencies = comp_stream,
        )
    end

    if num_state_gradient_flux > 0 || nhyperviscstate > 0

        comp_stream = volume_gradients!(device, (Nq, Nq))(
            balance_law,
            Val(dim),
            Val(N),
            dg.diffusion_direction,
            state_conservative.data,
            state_gradient_flux.data,
            Qhypervisc_grad.data,
            state_auxiliary.data,
            grid.vgeo,
            t,
            grid.D,
            hypervisc_indexmap,
            topology.realelems,
            ndrange = (Nq * nrealelem, Nq),
            dependencies = (comp_stream,),
        )

        comp_stream = interface_gradients!(device, workgroups_surface)(
            balance_law,
            Val(dim),
            Val(N),
            dg.diffusion_direction,
            dg.numerical_flux_gradient,
            state_conservative.data,
            state_gradient_flux.data,
            Qhypervisc_grad.data,
            state_auxiliary.data,
            grid.vgeo,
            grid.sgeo,
            t,
            grid.vmap⁻,
            grid.vmap⁺,
            grid.elemtobndy,
            hypervisc_indexmap,
            grid.interiorelems;
            ndrange = ndrange_interior_surface,
            dependencies = (comp_stream,),
        )

        if communicate
            exchange_state_conservative = MPIStateArrays.end_ghost_exchange!(
                state_conservative;
                dependencies = exchange_state_conservative,
            )

            # update_aux may start asynchronous work on the compute device and
            # we synchronize those here through a device event.
            wait(device, exchange_state_conservative)
            update_auxiliary_state!(
                dg,
                balance_law,
                state_conservative,
                t,
                dg.grid.topology.ghostelems,
            )
            exchange_state_conservative = Event(device)
        end

        comp_stream = interface_gradients!(device, workgroups_surface)(
            balance_law,
            Val(dim),
            Val(N),
            dg.diffusion_direction,
            dg.numerical_flux_gradient,
            state_conservative.data,
            state_gradient_flux.data,
            Qhypervisc_grad.data,
            state_auxiliary.data,
            grid.vgeo,
            grid.sgeo,
            t,
            grid.vmap⁻,
            grid.vmap⁺,
            grid.elemtobndy,
            hypervisc_indexmap,
            grid.exteriorelems;
            ndrange = ndrange_exterior_surface,
            dependencies = (comp_stream, exchange_state_conservative),
        )

        if communicate
            if num_state_gradient_flux > 0
                exchange_state_gradient_flux =
                    MPIStateArrays.begin_ghost_exchange!(
                        state_gradient_flux,
                        dependencies = comp_stream,
                    )
            end
            if nhyperviscstate > 0
                exchange_Qhypervisc_grad = MPIStateArrays.begin_ghost_exchange!(
                    Qhypervisc_grad,
                    dependencies = comp_stream,
                )
            end
        end

        if num_state_gradient_flux > 0
            # update_aux_diffusive may start asynchronous work on the compute device
            # and we synchronize those here through a device event.
            wait(device, comp_stream)
            update_auxiliary_state_gradient!(
                dg,
                balance_law,
                state_conservative,
                t,
                dg.grid.topology.realelems,
            )
            comp_stream = Event(device)
        end
    end

    if nhyperviscstate > 0
        #########################
        # Laplacian Computation #
        #########################

        comp_stream =
            volume_divergence_of_gradients!(device, workgroups_volume)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                grid.vgeo,
                grid.D,
                topology.realelems;
                ndrange = ndrange_volume,
                dependencies = (comp_stream,),
            )

        comp_stream =
            interface_divergence_of_gradients!(device, workgroups_surface)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                CentralNumericalFluxDivergence(),
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                grid.vgeo,
                grid.sgeo,
                grid.vmap⁻,
                grid.vmap⁺,
                grid.elemtobndy,
                grid.interiorelems;
                ndrange = ndrange_interior_surface,
                dependencies = (comp_stream,),
            )

        if communicate
            exchange_Qhypervisc_grad = MPIStateArrays.end_ghost_exchange!(
                Qhypervisc_grad,
                dependencies = exchange_Qhypervisc_grad,
            )
        end

        comp_stream =
            interface_divergence_of_gradients!(device, workgroups_surface)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                CentralNumericalFluxDivergence(),
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                grid.vgeo,
                grid.sgeo,
                grid.vmap⁻,
                grid.vmap⁺,
                grid.elemtobndy,
                grid.exteriorelems;
                ndrange = ndrange_exterior_surface,
                dependencies = (comp_stream, exchange_Qhypervisc_grad),
            )

        if communicate
            exchange_Qhypervisc_div = MPIStateArrays.begin_ghost_exchange!(
                Qhypervisc_div,
                dependencies = comp_stream,
            )
        end

        ####################################
        # Hyperdiffusive terms computation #
        ####################################

        comp_stream =
            volume_gradients_of_laplacians!(device, workgroups_volume)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                state_conservative.data,
                state_auxiliary.data,
                grid.vgeo,
                grid.ω,
                grid.D,
                topology.realelems,
                t;
                ndrange = ndrange_volume,
                dependencies = (comp_stream,),
            )

        comp_stream =
            interface_gradients_of_laplacians!(device, workgroups_surface)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                CentralNumericalFluxHigherOrder(),
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                state_conservative.data,
                state_auxiliary.data,
                grid.vgeo,
                grid.sgeo,
                grid.vmap⁻,
                grid.vmap⁺,
                grid.elemtobndy,
                grid.interiorelems,
                t;
                ndrange = ndrange_interior_surface,
                dependencies = (comp_stream,),
            )

        if communicate
            exchange_Qhypervisc_div = MPIStateArrays.end_ghost_exchange!(
                Qhypervisc_div,
                dependencies = exchange_Qhypervisc_div,
            )
        end

        comp_stream =
            interface_gradients_of_laplacians!(device, workgroups_surface)(
                balance_law,
                Val(dim),
                Val(N),
                dg.diffusion_direction,
                CentralNumericalFluxHigherOrder(),
                Qhypervisc_grad.data,
                Qhypervisc_div.data,
                state_conservative.data,
                state_auxiliary.data,
                grid.vgeo,
                grid.sgeo,
                grid.vmap⁻,
                grid.vmap⁺,
                grid.elemtobndy,
                grid.exteriorelems,
                t;
                ndrange = ndrange_exterior_surface,
                dependencies = (comp_stream, exchange_Qhypervisc_div),
            )

        if communicate
            exchange_Qhypervisc_grad = MPIStateArrays.begin_ghost_exchange!(
                Qhypervisc_grad,
                dependencies = comp_stream,
            )
        end
    end

    ###################
    # RHS Computation #
    ###################
    comp_stream = volume_tendency!(device, (Nq, Nq))(
        balance_law,
        Val(dim),
        Val(N),
        dg.direction,
        tendency.data,
        state_conservative.data,
        state_gradient_flux.data,
        Qhypervisc_grad.data,
        state_auxiliary.data,
        grid.vgeo,
        t,
        grid.ω,
        grid.D,
        topology.realelems,
        increment;
        ndrange = (nrealelem * Nq, Nq),
        dependencies = (comp_stream,),
    )

    comp_stream = interface_tendency!(device, workgroups_surface)(
        balance_law,
        Val(dim),
        Val(N),
        dg.direction,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        tendency.data,
        state_conservative.data,
        state_gradient_flux.data,
        Qhypervisc_grad.data,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.interiorelems;
        ndrange = ndrange_interior_surface,
        dependencies = (comp_stream,),
    )

    if communicate
        if num_state_gradient_flux > 0 || nhyperviscstate > 0
            if num_state_gradient_flux > 0
                exchange_state_gradient_flux =
                    MPIStateArrays.end_ghost_exchange!(
                        state_gradient_flux;
                        dependencies = exchange_state_gradient_flux,
                    )

                # update_aux_diffusive may start asynchronous work on the
                # compute device and we synchronize those here through a device
                # event.
                wait(device, exchange_state_gradient_flux)
                update_auxiliary_state_gradient!(
                    dg,
                    balance_law,
                    state_conservative,
                    t,
                    dg.grid.topology.ghostelems,
                )
                exchange_state_gradient_flux = Event(device)
            end
            if nhyperviscstate > 0
                exchange_Qhypervisc_grad = MPIStateArrays.end_ghost_exchange!(
                    Qhypervisc_grad;
                    dependencies = exchange_Qhypervisc_grad,
                )
            end
        else
            exchange_state_conservative = MPIStateArrays.end_ghost_exchange!(
                state_conservative;
                dependencies = exchange_state_conservative,
            )

            # update_aux may start asynchronous work on the compute device and
            # we synchronize those here through a device event.
            wait(device, exchange_state_conservative)
            update_auxiliary_state!(
                dg,
                balance_law,
                state_conservative,
                t,
                dg.grid.topology.ghostelems,
            )
            exchange_state_conservative = Event(device)
        end
    end

    comp_stream = interface_tendency!(device, workgroups_surface)(
        balance_law,
        Val(dim),
        Val(N),
        dg.direction,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        tendency.data,
        state_conservative.data,
        state_gradient_flux.data,
        Qhypervisc_grad.data,
        state_auxiliary.data,
        grid.vgeo,
        grid.sgeo,
        t,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        grid.exteriorelems;
        ndrange = ndrange_exterior_surface,
        dependencies = (
            comp_stream,
            exchange_state_conservative,
            exchange_state_gradient_flux,
            exchange_Qhypervisc_grad,
        ),
    )

    # The synchronization here through a device event prevents CuArray based and
    # other default stream kernels from launching before the work scheduled in
    # this function is finished.
    wait(device, comp_stream)
end

function init_ode_state(dg::DGModel, args...; init_on_cpu = false)
    device = arraytype(dg.grid) <: Array ? CPU() : CUDA()

    balance_law = dg.balancelaw
    grid = dg.grid

    state_conservative = create_conservative_state(balance_law, grid)

    topology = grid.topology
    Np = dofs_per_element(grid)

    state_auxiliary = dg.state_auxiliary
    dim = dimensionality(grid)
    N = polynomialorder(grid)
    nrealelem = length(topology.realelems)

    if !init_on_cpu
        event = Event(device)
        event = kernel_init_state_conservative!(device, min(Np, 1024))(
            balance_law,
            Val(dim),
            Val(N),
            state_conservative.data,
            state_auxiliary.data,
            grid.vgeo,
            topology.realelems,
            args...;
            ndrange = Np * nrealelem,
            dependencies = (event,),
        )
        wait(device, event)
    else
        h_state_conservative = similar(state_conservative, Array)
        h_state_auxiliary = similar(state_auxiliary, Array)
        h_state_auxiliary .= state_auxiliary
        event = kernel_init_state_conservative!(CPU(), Np)(
            balance_law,
            Val(dim),
            Val(N),
            h_state_conservative.data,
            h_state_auxiliary.data,
            Array(grid.vgeo),
            topology.realelems,
            args...;
            ndrange = Np * nrealelem,
        )
        wait(event) # XXX: This could be `wait(device, event)` once KA supports that.
        state_conservative .= h_state_conservative
    end

    event = Event(device)
    event = MPIStateArrays.begin_ghost_exchange!(
        state_conservative;
        dependencies = event,
    )
    event = MPIStateArrays.end_ghost_exchange!(
        state_conservative;
        dependencies = event,
    )
    wait(device, event)

    return state_conservative
end

function restart_ode_state(dg::DGModel, state_data; init_on_cpu = false)
    bl = dg.balancelaw
    grid = dg.grid

    state = create_state(bl, grid)
    state .= state_data

    device = arraytype(dg.grid) <: Array ? CPU() : CUDA()
    event = Event(device)
    event = MPIStateArrays.begin_ghost_exchange!(state; dependencies = event)
    event = MPIStateArrays.end_ghost_exchange!(state; dependencies = event)
    wait(device, event)

    return state
end

function restart_auxiliary_state(bl, grid, aux_data)
    state_auxiliary = create_auxiliary_state(bl, grid)
    state_auxiliary .= aux_data
    return state_auxiliary
end

# fallback
function update_auxiliary_state!(
    dg::DGModel,
    balance_law::BalanceLaw,
    state_conservative::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end

function update_auxiliary_state_gradient!(
    dg::DGModel,
    balance_law::BalanceLaw,
    state_conservative::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end

function indefinite_stack_integral!(
    dg::DGModel,
    m::BalanceLaw,
    state_conservative::MPIStateArray,
    state_auxiliary::MPIStateArray,
    t::Real,
    elems::UnitRange = dg.grid.topology.elems,
)

    device = typeof(state_conservative.data) <: Array ? CPU() : CUDA()

    grid = dg.grid
    topology = grid.topology

    dim = dimensionality(grid)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dim == 2 ? 1 : Nq

    FT = eltype(state_conservative)

    # do integrals
    nelem = length(elems)
    nvertelem = topology.stacksize
    horzelems = fld1(first(elems), nvertelem):fld1(last(elems), nvertelem)

    event = Event(device)
    event = kernel_indefinite_stack_integral!(device, (Nq, Nqk))(
        m,
        Val(dim),
        Val(N),
        Val(nvertelem),
        state_conservative.data,
        state_auxiliary.data,
        grid.vgeo,
        grid.Imat,
        horzelems;
        ndrange = (length(horzelems) * Nq, Nqk),
        dependencies = (event,),
    )
    wait(device, event)
end

function reverse_indefinite_stack_integral!(
    dg::DGModel,
    m::BalanceLaw,
    state_conservative::MPIStateArray,
    state_auxiliary::MPIStateArray,
    t::Real,
    elems::UnitRange = dg.grid.topology.elems,
)

    device = typeof(state_auxiliary.data) <: Array ? CPU() : CUDA()

    grid = dg.grid
    topology = grid.topology

    dim = dimensionality(grid)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dim == 2 ? 1 : Nq

    FT = eltype(state_auxiliary)

    # do integrals
    nelem = length(elems)
    nvertelem = topology.stacksize
    horzelems = fld1(first(elems), nvertelem):fld1(last(elems), nvertelem)

    event = Event(device)
    event = kernel_reverse_indefinite_stack_integral!(device, (Nq, Nqk))(
        m,
        Val(dim),
        Val(N),
        Val(nvertelem),
        state_conservative.data,
        state_auxiliary.data,
        horzelems;
        ndrange = (length(horzelems) * Nq, Nqk),
        dependencies = (event,),
    )
    wait(device, event)
end

function nodal_update_auxiliary_state!(
    f!,
    dg::DGModel,
    m::BalanceLaw,
    state_conservative::MPIStateArray,
    t::Real,
    elems::UnitRange = dg.grid.topology.realelems;
    diffusive = false,
)
    device = typeof(state_conservative.data) <: Array ? CPU() : CUDA()

    grid = dg.grid
    topology = grid.topology

    dim = dimensionality(grid)
    N = polynomialorder(grid)
    Nq = N + 1
    nelem = length(elems)

    Np = dofs_per_element(grid)

    nodal_update_auxiliary_state! =
        kernel_nodal_update_auxiliary_state!(device, min(Np, 1024))
    ### update state_auxiliary variables
    event = Event(device)
    if diffusive
        event = nodal_update_auxiliary_state!(
            m,
            Val(dim),
            Val(N),
            f!,
            state_conservative.data,
            dg.state_auxiliary.data,
            dg.state_gradient_flux.data,
            t,
            elems,
            grid.activedofs;
            ndrange = Np * nelem,
            dependencies = (event,),
        )
    else
        event = nodal_update_auxiliary_state!(
            m,
            Val(dim),
            Val(N),
            f!,
            state_conservative.data,
            dg.state_auxiliary.data,
            t,
            elems,
            grid.activedofs;
            ndrange = Np * nelem,
            dependencies = (event,),
        )
    end
    wait(device, event)
end

"""
    courant(local_courant::Function, dg::DGModel, m::BalanceLaw,
            state_conservative::MPIStateArray, direction=EveryDirection())
Returns the maximum of the evaluation of the function `local_courant`
pointwise throughout the domain.  The function `local_courant` is given an
approximation of the local node distance `Δx`.  The `direction` controls which
reference directions are considered when computing the minimum node distance
`Δx`.
An example `local_courant` function is
    function local_courant(m::AtmosModel, state_conservative::Vars, state_auxiliary::Vars,
                           diffusive::Vars, Δx)
      return Δt * cmax / Δx
    end
where `Δt` is the time step size and `cmax` is the maximum flow speed in the
model.
"""
function courant(
    local_courant::Function,
    dg::DGModel,
    m::BalanceLaw,
    state_conservative::MPIStateArray,
    Δt,
    simtime,
    direction = EveryDirection(),
)
    grid = dg.grid
    topology = grid.topology
    nrealelem = length(topology.realelems)

    if nrealelem > 0
        N = polynomialorder(grid)
        dim = dimensionality(grid)
        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq
        device = grid.vgeo isa Array ? CPU() : CUDA()
        pointwise_courant = similar(grid.vgeo, Nq^dim, nrealelem)
        event = Event(device)
        event = Grids.kernel_min_neighbor_distance!(device, (Nq, Nq, Nqk))(
            Val(N),
            Val(dim),
            direction,
            pointwise_courant,
            grid.vgeo,
            topology.realelems;
            ndrange = (nrealelem * Nq, Nq, Nqk),
            dependencies = (event,),
        )
        event = kernel_local_courant!(device, Nq * Nq * Nqk)(
            m,
            Val(dim),
            Val(N),
            pointwise_courant,
            local_courant,
            state_conservative.data,
            dg.state_auxiliary.data,
            dg.state_gradient_flux.data,
            topology.realelems,
            Δt,
            simtime,
            direction;
            ndrange = nrealelem * Nq * Nq * Nqk,
            dependencies = (event,),
        )
        wait(device, event)
        rank_courant_max = maximum(pointwise_courant)
    else
        rank_courant_max = typemin(eltype(state_conservative))
    end

    MPI.Allreduce(rank_courant_max, max, topology.mpicomm)
end

function copy_stack_field_down!(
    dg::DGModel,
    m::BalanceLaw,
    state_auxiliary::MPIStateArray,
    fldin,
    fldout,
    elems = topology.elems,
)

    device = typeof(state_auxiliary.data) <: Array ? CPU() : CUDA()

    grid = dg.grid
    topology = grid.topology

    dim = dimensionality(grid)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dim == 2 ? 1 : Nq

    # do integrals
    nelem = length(elems)
    nvertelem = topology.stacksize
    horzelems = fld1(first(elems), nvertelem):fld1(last(elems), nvertelem)

    event = Event(device)
    event = kernel_copy_stack_field_down!(device, (Nq, Nqk))(
        Val(dim),
        Val(N),
        Val(nvertelem),
        state_auxiliary.data,
        horzelems,
        Val(fldin),
        Val(fldout);
        ndrange = (length(horzelems) * Nq, Nqk),
        dependencies = (event,),
    )
    wait(device, event)
end

function MPIStateArrays.MPIStateArray(dg::DGModel)
    balance_law = dg.balancelaw
    grid = dg.grid

    state_conservative = create_conservative_state(balance_law, grid)

    return state_conservative
end

function create_hypervisc_indexmap(balance_law::BalanceLaw)
    # helper function
    _getvars(v, ::Type) = v
    function _getvars(v::Vars, ::Type{T}) where {T <: NamedTuple}
        fields = getproperty.(Ref(v), fieldnames(T))
        collect(Iterators.Flatten(_getvars.(fields, fieldtypes(T))))
    end

    gradvars = vars_state_gradient(balance_law, Int)
    gradlapvars = vars_gradient_laplacian(balance_law, Int)
    indices = Vars{gradvars}(1:varsize(gradvars))
    SVector{varsize(gradlapvars)}(_getvars(indices, gradlapvars))
end


text:
struct FBAOut{T<:Real} <: AbstractMetState
    v::Vector{T} # Flux vector solution
    obj_val::T # The value of the objective
    obj_ider # The used objective identifier
    sol # The LP solution object
end

# empty FBAOut
objval(out::FBAOut) = out.obj_val
objider(out::FBAOut) = out.obj_ider
get_lp_model(out::FBAOut) = out.sol


text:
using DataFrames
import CSV
f = readdir()
filter!(n -> startswith(n, "svd-pred"), f)
filter!(n -> endswith(n, "csv"), f)
predict = vcat(DataFrame.(CSV.File.(f))...)
CSV.write("predictions.csv", predict)

max_stars_repo_name: viralemergence/trefle

max_issues_repo_name: viralemergence/trefle

max_forks_repo_name: viralemergence/trefle


text:
# module Octo.Schema

export InvalidChangesetError
export validate_length

struct InvalidChangesetError <: Exception
    msg
end

function validates(M, nt::NamedTuple)
    Tname = Base.typename(M)
    if haskey(validation_models, Tname)
        validations = validation_models[Tname]
        validations((M, nt))
    end
end

function validate_length((M, nt), field::Symbol; kwargs...)
    if haskey(nt, field)
        value = nt[field]
        check = (min= (value, num) -> length(value) >= num,)
        for (key, num) in kwargs
            if haskey(check, key)
                if check[key](value, num)
                else
                    throw(InvalidChangesetError(string(field, " length must be ", key, " = ", num)))
                end
            end
        end
    end
end

# end module Octo.Schema


text:
#**************************************************************************************
# test_Calculations.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file contains all the required function to test calculate the Geometrical efficiency.
#
#**************************************************************************************

const absoluteTol1 = 1.0e-14 #1.0e-16
const absoluteTol2 = 1.0e-11 #1.0e-13

let poly1(z::Float64) = @evalpoly(z, 1.0, 2.0),
	poly2(z::Float64) = @evalpoly(z, 1.0, 2.0, 3.0)

	@testset "integrate start=$str, end=$nd" for
	str::Float64 = -30.0:5.0:30.0,
	nd::Float64  = -30.0:5.0:30.0
	#str === nd &&  continue

		#@test G.integrate(poly0, str, nd)[1] ≈ @evalpoly(nd, 0.0, 1.0) - @evalpoly(str, 0.0, 1.0) atol=absoluteTol2
		@test G.integrate(poly1, str, nd)[1] ≈ @evalpoly(nd, 0.0, 1.0, 1.0) - @evalpoly(str, 0.0, 1.0, 1.0) atol=absoluteTol2
		@test G.integrate(poly2, str, nd)[1] ≈ @evalpoly(nd, 0.0, 1.0, 1.0, 1.0) - @evalpoly(str, 0.0, 1.0, 1.0, 1.0) atol=absoluteTol2

		@test G.integrate(sin, str/10.0, nd/10.0)[1] ≈ cos(str/10.0) - cos(nd/10.0) atol=absoluteTol2
	end #testset_integrate
end #let

@testset "special cases" begin

	@debug("special case - point at the surface of cylindrical detector; very restrict test")
	@testset "cylindrical detector of cryRadius $cryRadius" for
	cryRadius::Float64    = 1.0:0.5:3.0,
	cryLegth::Float64 = 0.0:.5:2.0

	let acylDetector = CylDetector(cryRadius, cryLegth)

		@test geoEff(acylDetector, Point(0)) ≈ 0.5
		@test geoEff(acylDetector, Point(0, prevfloat(cryRadius)))  ≈ 0.5
		@test geoEff(acylDetector, Point(0, nextfloat(-cryRadius))) ≈ 0.5
		@test geoEff(acylDetector, Point(0,         cryRadius/2.0)) ≈ 0.5
		@test geoEff(acylDetector, Point(0,        -cryRadius/2.0)) ≈ 0.5

		@test geoEff(acylDetector, Point(eps())) ≈ 0.5
		@test 0.25 < geoEff(acylDetector, Point(eps(), prevfloat(cryRadius)))  <= 0.5	#  0.37485153429958035
		@test 0.25 < geoEff(acylDetector, Point(eps(), nextfloat(-cryRadius))) <= 0.5	#  0.3749998497981526
		@test geoEff(acylDetector, Point(eps(),          cryRadius/2.0)) ≈ 0.5
		@test geoEff(acylDetector, Point(eps(),         -cryRadius/2.0)) ≈ 0.5

		@test geoEff(acylDetector, Point(eps(cryRadius), cryRadius))  ≈ 0.25
		@test 0.0 < geoEff(acylDetector, Point(cryRadius, cryRadius)) < 0.25

		@test_throws G.NotImplementedError	geoEff(acylDetector, Point(eps(cryRadius), nextfloat(cryRadius)))
		@test_throws G.NotImplementedError	geoEff(acylDetector, Point(nextfloat(cryRadius), nextfloat(cryRadius)))
		@test_throws G.InValidGeometry		geoEff(acylDetector, Point(-1, 0))
		@test_throws G.InValidGeometry		geoEff(acylDetector, Point(-1, prevfloat(cryRadius)))
	end #let
	end #testset_cylindrical_detector

	@debug("special case - point at the surface of Borehole detector")
	@testset "Borehole detector of cryRadius $cryRadius and height $height, k $k" for
	cryRadius::Float64 = 1.0:0.5:11.0,
	height::Float64    = 1.0:0.5:11.0,
	k::Float64         = 1.1:0.5:11.0

	let holeradius::Float64 = cryRadius/k,		# k > 1
		aboreDetector = BoreDetector(cryRadius, height, holeradius)

		@test  0.0 < geoEff(aboreDetector, Point(0.0))        < 1.0
		@test  0.0 < geoEff(aboreDetector, Point(-0.1))       < 1.0 ### invert Detector
		@test  0.0 < geoEff(aboreDetector, Point(height/2.0))  < 1.0
		@test  0.0 < geoEff(aboreDetector, Point(-height/2.0)) < 1.0
		@test_skip  geoEff(aboreDetector, Point(height/2.0)) ≈ geoEff(aboreDetector, Point(-height/2.0)) #Test Failed

		@test   0.0 < geoEff(aboreDetector, Point(height))  < 1.0
		@test_skip   0.0 < geoEff(aboreDetector, Point(-height)) < 1.0	#StackOverflowError
		@test_skip   geoEff(aboreDetector, Point(height)) ≈ geoEff(aboreDetector, Point(-height)) #StackOverflowError

		@test  0.0 < geoEff(aboreDetector, Point(1.5*height))   < 1.0
		@test_skip   0.0 < geoEff(aboreDetector, Point(-1.5*height))  < 1.0 #StackOverflowError
		@test_skip   geoEff(aboreDetector, Point(1.5*height)) ≈ geoEff(aboreDetector, Point(-1.5*height)) #StackOverflowError
	end #let
	end #testset_Borehole_detector

	@debug("special case - point at the surface of well detector")
	@testset "Well detectors of cryRadius $cryRadius and height $height, k $k" for
	cryRadius::Float64 = 1.0:0.5:11.0,
	height::Float64    = 1.0:0.5:11.0,
	k::Float64         = 1.1:0.5:11.0

	let holeradius::Float64 = cryRadius/k,		# k > 1
		welldepth::Float64 = height/k,			# k > 1
		awellDetector = WellDetector(cryRadius, height, holeradius, welldepth)

		@test geoEff(awellDetector, Point(welldepth)) ≈ 0.5
		@test geoEff(awellDetector, Point(welldepth, prevfloat(holeradius))) ≈ 0.5
		@test geoEff(awellDetector, Point(welldepth, nextfloat(-holeradius))) ≈ 0.5
		@test geoEff(awellDetector, Point(welldepth, holeradius/2.0)) ≈ 0.5
		@test geoEff(awellDetector, Point(welldepth, -holeradius/2.0)) ≈ 0.5

		@test geoEff(awellDetector, Point(nextfloat(welldepth))) ≈ 0.5
		@test geoEff(awellDetector, Point(nextfloat(welldepth), prevfloat(holeradius))) ≈ 0.5
		@test geoEff(awellDetector, Point(nextfloat(welldepth), prevfloat(-holeradius))) ≈ 0.5
		@test geoEff(awellDetector, Point(nextfloat(welldepth),  holeradius/2.0)) ≈ 0.5
		@test geoEff(awellDetector, Point(nextfloat(welldepth), -holeradius/2.0)) ≈ 0.5

	end #let
	end #testset_Well_detector
end #testset_spectial_cases

@testset "scaling" begin

	@debug("scaling test - cylindrical detector with point source")
	@testset "[J=$j] test, Scaling $k at cryRadius=$cryRadius" for
	cryRadius = [1,2,3,4,5,6,7,8,9,10.1,10.5,10.6],
	j=2:5:100, 	# j > 1
	k=2:5:100

		acylDetector  = CylDetector(  cryRadius)
		acylDetectork = CylDetector(k*cryRadius)
		axPnt  = Point(  cryRadius/j); naxPnt  = Point(  cryRadius/j,   cryRadius/j)
		axPntk = Point(k*cryRadius/j); naxPntk = Point(k*cryRadius/j, k*cryRadius/j)

		@test geoEff(acylDetector, axPnt)  ≈ geoEff(acylDetectork, axPntk)		# axial point
		@test geoEff(acylDetector, naxPnt) ≈ geoEff(acylDetectork, naxPntk)		# non-axial point
	end #testset_cylindrical_detector

	@debug("scaling test - Borehole detector with point source")
	@testset "cryRadius=$cryRadius, holeRadius=$holeRadius" for
	cryRadius  = [1,2,3,4,5,6,7,8,9,10.1,10.5,10.6],
	holeRadius::Float64 = [1,2,3,4,5,6,7,8,9,10.1,10.5,10.6]/2.2
	holeRadius > cryRadius && continue

		for j=2:100, k=2:100
		aboreDetector  = BoreDetector(  cryRadius,   j,   holeRadius)
		aboreDetectork = BoreDetector(k*cryRadius, k*j, k*holeRadius)
		axPnt  = Point(  cryRadius/j); naxPnt  = Point(  cryRadius/j,   holeRadius/j)
		axPntk = Point(k*cryRadius/j); naxPntk = Point(k*cryRadius/j, k*holeRadius/j)

			@test geoEff(aboreDetector , axPnt)  ≈ geoEff(aboreDetectork, axPntk) atol= absoluteTol1	# axial point
			@test geoEff(aboreDetector , naxPnt) ≈ geoEff(aboreDetectork, naxPntk) atol= absoluteTol1	# non-axial point
		end #for
	end #testset_Borehole_detector

	@debug("scaling test - Well-type detector with point source")
	@testset "cryRadius=$cryRadius, holeRadius=$holeRadius" for
	cryRadius  = [1,2,3,4,5,6,7,8,9,10.1,10.5,10.6],
	holeRadius = [1,2,3,4,5,6,7,8,9,10.1,10.5,10.6]/2.2
	holeRadius > cryRadius && continue

		for j=2:5:100, k=2:5:100
			awellDetector  = WellDetector(  cryRadius,   j,   holeRadius,   j/2.0)
			awellDetectork = WellDetector(k*cryRadius, k*j, k*holeRadius, k*j/2.0)
			axPnt  = Point(  cryRadius/j); naxPnt  = Point(  cryRadius/j,   holeRadius/j)
			axPntk = Point(k*cryRadius/j); naxPntk = Point(k*cryRadius/j, k*holeRadius/j)

			@test geoEff(awellDetector , axPnt)  ≈ geoEff(awellDetectork, axPntk) atol= absoluteTol2	# axial point
			@test geoEff(awellDetector , naxPnt) ≈ geoEff(awellDetectork, naxPntk) atol= absoluteTol2	# non-axial point
		end #for
	end #testset_Well-type_detector
end #tesset_scalling

@testset "geoEff-relaxed`" begin

	@testset "geoEff on CylDetector(5,$cryLength)" for
	cryLength=[0,1,5,10]

		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1, 1))    < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1, 1//2)) < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1, pi))   < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1, 1.0))  < 0.5

		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1//2, 1))    < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1//2, 1//2)) < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1//2, pi))   < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1//2, 1.0))  < 0.5

		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),e, 1))    < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),e, 1//2)) < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),e, pi))   < 0.5 #
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),e, 1.0))  < 0.5

		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1.0, 1))    < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1.0, 1//2)) < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1.0, pi))   < 0.5
		@test 0.0 < geoEff(Detector(5,cryLength),(Point(1),1.0, 1.0))  < 0.5
     end #testset_CylDetector

	@testset "geoEff on WellDetector(5, 4, $holeRadius, $holeDepth)" for
	holeRadius = 3:0.5:4,
	holeDepth  = 0.1:1.0:3.1

	local (mim, mam) = holeDepth >= 1 ? (0.5, 1.0) : (0.0, 0.5)

		@test geoEff(Detector(5, 4, holeRadius, 1),(Point(1),0, 0)) ≈ 0.5 atol= absoluteTol1
		@test mim < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),0, 0))    < mam

		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1, 1))    < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1, 1//2)) < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1, pi))   < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1, 1.0))  < 1.0

		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1//2, 1))    < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1//2, 1//2)) < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1//2, pi))   < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1//2, 1.0))  < 1.0

		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),e, 1))    < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),e, 1//2)) < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),e, pi))   < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),e, 1.0))  < 1.0

		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1.0, 1))    < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1.0, 1//2)) < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1.0, pi))   < 1.0
		@test 0.0 < geoEff(Detector(5, 4, holeRadius, holeDepth),(Point(1),1.0, 1.0))  < 1.0
	end #testset_WellDetector
end #testset_geoEff_relaxed

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
using NeuroSignals
using Test

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


text:
function draw_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.draw_mode = getproperty(han.b["draw_button"],:active,Bool)

    nothing
end

function draw_start(han::Tracker_Handles,x,y)

    plot_image(han,han.current_frame')
    r = Gtk.getgc(han.c)
    Cairo.save(r)
    ctxcopy = copy(r)

    #if han.tracked[han.frame]
        new_whisker=Whisker1()

        push!(new_whisker.x,x)
        push!(new_whisker.y,y)
        push!(new_whisker.scores,1.0f0)
        push!(new_whisker.thick,1.0f0)

        new_whisker.len=1

        han.wt.whiskers=[new_whisker]
        han.woi_id = 1
    #end

    plot_whiskers(han)

    push!((han.c.mouse, :button1motion),  (c, event) -> draw_move(han, event.x, event.y, ctxcopy))
    push!((han.c.mouse, :motion), Gtk.default_mouse_cb)
    push!((han.c.mouse, :button1release), (c, event) -> draw_stop(han, event.x, event.y, ctxcopy))

    nothing
end

function draw_move(han::Tracker_Handles, x,y,ctxcopy)

    r=Gtk.getgc(han.c)

    han.wt.whiskers[han.woi_id].len+=1

    front_dist = (han.wt.whiskers[han.woi_id].x[1]-x)^2+(han.wt.whiskers[han.woi_id].y[1]-y)^2
    end_dist = (han.wt.whiskers[han.woi_id].x[end]-x)^2+(han.wt.whiskers[han.woi_id].y[end]-y)^2

    if end_dist<front_dist #drawing closer to end

        push!(han.wt.whiskers[han.woi_id].x,x)
        push!(han.wt.whiskers[han.woi_id].y,y)
        push!(han.wt.whiskers[han.woi_id].thick,1.0)
        push!(han.wt.whiskers[han.woi_id].scores,1.0)

    else
        unshift!(han.wt.whiskers[han.woi_id].x,x)
        unshift!(han.wt.whiskers[han.woi_id].y,y)
        unshift!(han.wt.whiskers[han.woi_id].thick,1.0)
        unshift!(han.wt.whiskers[han.woi_id].scores,1.0)
    end

    #redraw whisker
    set_source(r,ctxcopy)
    paint(r)

    plot_whiskers(han)

    nothing
end

function draw_stop(han::Tracker_Handles,x,y,ctxcopy)

    assign_woi(han)
    han.tracked[han.frame]=true

    pop!((han.c.mouse, :button1motion))
    pop!((han.c.mouse, :motion))
    pop!((han.c.mouse, :button1release))

    nothing
end


text:
"""
Parameterized lower triangular matrices.
"""
abstract ParamLowerTriangular{T,S<:AbstractMatrix} <: Base.LinAlg.AbstractTriangular{T,S}

"""
Parameterized lower triangular matrix in which each element of the lower triangle is a parameter
"""
immutable ColMajorLowerTriangular{T,S<:AbstractMatrix} <: ParamLowerTriangular{T,S}
    Lambda::LowerTriangular{T,S}
end

ColMajorLowerTriangular(typ,n::Integer) = ColMajorLowerTriangular(LowerTriangular(eye(typ,n)))

ColMajorLowerTriangular(n::Integer) = ColMajorLowerTriangular(LowerTriangular(eye(n)))

Base.convert(::Type{LowerTriangular},A::ColMajorLowerTriangular) = A.Lambda

Base.size(A::ColMajorLowerTriangular, args...) = size(A.Lambda, args...)

Base.size(A::ColMajorLowerTriangular) = size(A.Lambda)

Base.copy(A::ColMajorLowerTriangular) = ColMajorLowerTriangular(copy(A.Lambda))

Base.copy!(A::ColMajorLowerTriangular,B::ColMajorLowerTriangular) = (copy!(A.Lambda.data,B.Lambda.data);A)

Base.full(A::ColMajorLowerTriangular) = full(A.Lambda)

@inline nlower(n::Integer) = (n*(n+1))>>1

function Base.getindex{T}(A::ColMajorLowerTriangular{T},s::Symbol)
    s == :θ || throw(KeyError(s))
    Ad = A.Lambda.data
    n = size(Ad,1)
    res = Array(T,nlower(n))
    k = 0
    for j = 1:n, i in j:n
        @inbounds res[k += 1] = Ad[i,j]
    end
    res
end

Base.getindex(A::ColMajorLowerTriangular,i::Integer,j::Integer) = A.Lambda[i,j]

function Base.setindex!{T}(A::ColMajorLowerTriangular{T},v::StridedVector{T},s::Symbol)
    s == :θ || throw(KeyError(s))
    Ad = A.Lambda.data
    n = Base.LinAlg.chksquare(Ad)
    if length(v) ≠ nlower(n)
        throw(DimensionMismatch("length(v) = $(length(v)), should be $(nlower(n))"))
    end
    k = 0
    for j in 1:n, i in j:n
        Ad[i,j] = v[k += 1]
    end
    A
end

"""
lower bounds on the parameters
"""
function lowerbd{T}(A::ColMajorLowerTriangular{T})
    n = size(A.Lambda.data,1)
    res = fill(convert(T,-Inf),nlower(n))
    k = -n
    for j in n+1:-1:2
        res[k += j] = zero(T)
    end
    res
end

"""
length of the parameter vector for the term
"""
nθ(A::ColMajorLowerTriangular) = nlower(size(A.Lambda.data,1))

function Base.scale!(A::ColMajorLowerTriangular,B::HBlkDiag)
    Ba = B.arr
    r,s,k = size(Ba)
    Al = A.Lambda
    n = Base.LinAlg.chksquare(Al)
    n == r || throw(DimensionMismatch())
    if r == 1
        scale!(Ba,Al[1,1])
    else
        Ac_mul_B!(Al,reshape(Ba,(r,s*k)))
    end
    B
end

function Base.scale!{T}(A::ColMajorLowerTriangular{T},B::Diagonal{T})
    size(A,1) == 1 || throw(DimensionMismatch())
    scale!(A.Lambda.data[1,1],B.diag)
    B
end

LT(A::ReMat) = ColMajorLowerTriangular(eltype(A.z),1)

LT(A::VectorReMat) = (Az = A.z; ColMajorLowerTriangular(eltype(Az),size(Az,1)))

function Base.scale!{T}(A::ColMajorLowerTriangular{T},B::DenseVecOrMat{T})
    al = A.Lambda
    if (l = Base.LinAlg.chksquare(al)) == 1
        return scale!(al.data[1],B)
    end
    m,n = size(B,1),size(B,2)
    Ac_mul_B!(al,reshape(B,(l,div(m,l)*n)))
    B
end

function Base.scale!{T}(A::ColMajorLowerTriangular{T},B::SparseMatrixCSC{T})
    al = A.Lambda
    (l = Base.LinAlg.chksquare(al)) == 1 || error("Code not yet written")
    scale!(al[1],B.nzval)
    B
end

function Base.scale!{T}(A::SparseMatrixCSC{T},B::ColMajorLowerTriangular)
    bl = B.Lambda
    (l = Base.LinAlg.chksquare(bl)) == 1 || error("Code not yet written")
    scale!(A.nzval,bl[1])
    A
end

function Base.scale!{T}(A::Diagonal{T},B::ColMajorLowerTriangular{T})
    bl = B.Lambda
    if (l = Base.LinAlg.chksquare(bl)) ≠ 1
        throw(DimensionMismatch(
        "for scale!(A::Diagonal,B::ColMajorLowerTriangular) B must be 1×1"))
    end
    scale!(bl[1],A.diag)
    A
end

function Base.scale!{T}(A::HBlkDiag{T},B::ColMajorLowerTriangular{T})
    aa = A.arr
    r,s,k = size(aa)
    bl = B.Lambda
    for i in 1:k
        A_mul_B!(sub(aa,:,:,i),bl)
    end
    A
end

function Base.scale!{T}(A::StridedMatrix{T},B::ColMajorLowerTriangular{T})
    bl = B.Lambda
    l = Base.LinAlg.chksquare(bl)
    l == 1 && return scale!(A,bl.data[1])
    m,n = size(A)
    q,r = divrem(n,l)
    r == 0 || throw(DimensionMismatch("size(A,2) = $n must be a multiple of size(B,1) = $l"))
    for k in 0:(q-1)
        A_mul_B!(sub(A,:,k*l + (1:l)),bl)
    end
    A
end

immutable DiagonalLowerTriangular{T} <: ParamLowerTriangular{T}
    diag::Vector{T}
end

DiagonalLowerTriangular(typ,n::Integer) = DiagonalLowerTriangular(ones(typ,n))

DiagonalLowerTriangular(n::Integer) = DiagonalLowerTriangular(ones(n))

Base.full(A::DiagonalLowerTriangular) = Diagonal(A.diag)

Base.convert(::Type{LowerTriangular},A::DiagonalLowerTriangular) = LowerTriangular(full(A))

Base.size(A::DiagonalLowerTriangular) = (n=length(A.diag);(n,n))

function Base.getindex(A::DiagonalLowerTriangular,s::Symbol)
    s == :θ || throw(KeyError(s))
    copy(A.diag)
end

nθ(A::DiagonalLowerTriangular) = length(A.diag)

Base.getindex(A::DiagonalLowerTriangular,i::Integer,j::Integer) = full(A)[i,j]

function Base.setindex!{T}(A::DiagonalLowerTriangular{T},v::StridedVector{T},s::Symbol)
    s == :θ || throw(KeyError(s))
    copy!(A.diag,v)
    A
end

function lowerbd{T}(A::DiagonalLowerTriangular{T})
    zeros(A.diag)
end

function Base.scale!(A::DiagonalLowerTriangular,B::HBlkDiag)
    bb = B.arr
    r,s,k = size(bb)
    dd = A.diag
    r == length(dd) || throw(DimensionMismatch())
    if r == 1
        scale!(bb,dd[1])
    else
        scale!(dd,reshape(bb,(r,s*k)))
    end
    bb
end

function Base.scale!{T}(A::DiagonalLowerTriangular{T},B::DenseVecOrMat{T})
    m,n = size(B)
    dd = A.diag
    l = length(dd)
    scale!(dd,reshape(B,(l,div(m,l)*n)))
    B
end

function Base.scale!{T}(A::HBlkDiag,B::DiagonalLowerTriangular{T})
    aa = A.arr
    r,s,k = size(aa)
    for i in 1:k
        scale!(sub(aa,:,:,i),B.diag)
    end
    A
end

function Base.scale!{T}(A::StridedMatrix{T},B::DiagonalLowerTriangular{T})
    dd = B.diag
    l = length(dd)
    l == 1 && return scale!(A,dd[1])
    m,n = size(A)
    q,r = divrem(n,l)
    r == 0 || throw(DimensionMismatch("size(A,2) = $n must be a multiple of size(B,1) = $l"))
    for k in 0:(q-1)
        scale!(sub(A,:,k*l + (1:l)),dd)
    end
    A
end


text:
"""
S=[S1, S2, S3] = Magnitude of principal stresses
 Pp = Pore pressure
 mu = Coefficient of sliding friction
 a = trend of S1, except when S1 is vertical a = trend of SHmax minus 90 degrees
 b = -plunge of S1 (plunge is angle from horizontal)
 c = rake of S2, 0 if S1 or S3 is vertical, 90 if S2 is vertical
 faults = 2 column matrix with strike(0-360), dip(0-90)
          (use right-hand rule)), e.g faults=[208,60;36,62;...]
"""
function FaultStress(s,Pp,μ,a,b,c,faults)

n=size(faults,1)

str=faults[:,1].*(π./180)
dip=faults[:,2].*(π./180)
a=a*(π/180)
b=b*(π/180)
c=c*(π/180)

#Principal Stress Tensor
IM=Matrix{Float64}(I,3,3)
S=s.*IM

#Transformation from Principal Stress to Geografic Coordinates

R1=[cos(a)*cos(b) sin(a)*cos(b) -sin(b);
    cos(a)*sin(b)*sin(c)-sin(a)*cos(c) sin(a)*sin(b)*sin(c)+cos(a)*cos(c) cos(b)*sin(c);
    cos(a)*sin(b)*cos(c)+sin(a)*sin(c) sin(a)*sin(b)*cos(c)-cos(a)*sin(c) cos(b)*cos(c)]

Sg=R1'*S*R1

    Sn=zeros(n)
    τ=zeros(n)
    rake=zeros(n)
    for k=1:n

    # Transformation from geographic to fault coordinate system
    R2=[cos(str[k]) sin(str[k]) 0;
         sin(str[k])*cos(dip[k]) -cos(str[k])*cos(dip[k]) -sin(dip[k]);
         -sin(str[k])*sin(dip[k]) cos(str[k])*sin(dip[k]) -cos(dip[k])]

        Sf=R2*Sg*R2'

        Sn[k]=Sf[3,3]

        if Sf[3,2]>0.0
            rake[k]=atan(Sf[3,2]/Sf[3,1])

        elseif Sf[3,2]>0.0 && Sf[3,1]>0.0
            rake[k]=π-atan(Sf[3,2]/-Sf[3,1])
        else
            rake[k]=atan(-Sf[3,2]/-Sf[3,1])-π
        end

        R3=[cos(rake[k]) sin(rake[k]) 0;
            -sin(rake[k]) cos(rake[k]) 0;
            0 0 1]
        Sr=R3*Sf*R3'
        τ[k]=Sr[3,1]
    end
    τ=abs.(τ)
    σn=Sn.-Pp
    CFF=τ .- μ.*σn

    return σn, τ, CFF
end


text:
#@reexport using Compat
@reexport using RecipesBase
@reexport using ColorTypes
#using Compat: range

const mygreen = RGBA{Float64}(151/255,180/255,118/255,1)
const mygreen2 = RGBA{Float64}(113/255,161/255,103/255,1)
const myblue = RGBA{Float64}(74/255,144/255,226/255,1)

@recipe function plot(p::Polygon)
    z = [p.vert; p.vert[1]]
    linecolor --> mygreen
    fillrange --> 0
    fillcolor --> mygreen
    ratio := 1
    legend := :none
    grid := false
    x := real.(z)
    y := imag.(z)
    ()
end

@recipe function plot(m::ConformalMap)

  xmax = pop!(plotattributes, :xmax, nothing)
  xmaxc =  xmax == nothing ? 4.0 : xmax

  nspokes = pop!(plotattributes, :nspokes, nothing)
  nθ =  nspokes == nothing ? 20 : nspokes

  nrings = pop!(plotattributes, :nrings, nothing)
  nr =  nrings == nothing ? 10 : nrings

  rmax = sqrt(2)*xmaxc

  xmaxp = xmaxc*abs(m.ps.ccoeff[1])
  dxp = round(xmaxp)*0.5

  layout := (1,2)
  ratio := 1
  legend := :none
  xlims --> [(-xmaxc,xmaxc) (-xmaxp,xmaxp)]
  ylims --> [(-xmaxc,xmaxc) (-xmaxp,xmaxp)]
  framestyle := :frame
  grid := false

  # make the spokes
  dθ = 2π/nθ
  θg = 0
  nrg = 100
  rg = range(1,stop=rmax,length=nrg)
  for jθ in 1:nθ
    ζg = collect(rg*exp(im*θg))
    zg = m(ζg)

    @series begin
      subplot := 1
      color --> mygreen2

      real.(ζg), imag.(ζg)
    end
    @series begin
      subplot := 2
      color -->  mygreen2

      real.(zg), imag.(zg)
    end
    θg += dθ
  end

  # make the rings
  dr = (rmax-1)/nr
  rg = 1

  nθg = 241
  θg = range(0,stop=2π,length=nθg)
  for jr in 1:nr
    ζg = collect(rg*exp.(im*θg));
    zg = m(ζg)
    @series begin
      subplot := 1
      color -->  myblue

      real.(ζg), imag.(ζg)
    end
    @series begin
      subplot := 2
      color -->  myblue

      real.(zg), imag.(zg)
    end
    rg += dr
  end

end


text:
#= ATTENTION:
   your own folder is considered as the current working directory
   for running your solver
=#

include("../../LibSPP/librarySPP.jl")
include("./spp_dm1.jl")
include("../../../setSPP.jl")

function main()
    println("Etudiant(e)s : LE et OUSMAN")

    # Collecting the names of instances to solve located in the folder Data ----
    target = "../../Data"
    fnames = getfname(target)

    fres = splitdir(splitdir(pwd())[end-1])[end]
    io = open("../res/"*fres*".res", "w")
    for instance = 1:length(fnames)

        # Load one numerical instance ------------------------------------------
        C, A = loadSPP(string(target,"/",fnames[instance]))

        t1 = @elapsed x, zInit, Einit = GreedyConstruction(C, A)
        t2 = @elapsed xBest, zBest = GreedyImprovement(C, A, x, zInit, Einit)

        # Saving results -------------------------------------------------------
        println(io, fnames[instance], " ", zInit, " ", zBest, " ", t1, " ", t2, " ", t1+t2)
    end
    close(io)

end

main()

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
include("../src/HMat2d.jl");

function laplacian2d(nx,ny)
    M = speye(nx*ny);
    for x=1:nx
        for y=1:ny
            s = x+(y-1)*nx;
            M[s,s] = 4;
            if x > 1
                M[s,s-1] = -1;
            end
            if x < nx
                M[s,s+1] = -1;
            end
            if y > 1
                M[s,s-nx] = -1;
            end
            if y < ny
                M[s,s+nx] = -1;
            end
        end
    end
    return M;
end

EPS = 1e-6;
MaxRank = 2;
minn = 4;

n=64;
A = full(laplacian2d(n,n));
Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
A = inv(full(A[Z2C,Z2C]));
@time AH = HMat2dd2h(A, [n,n], [n,n], WEAK, [0,0], [0,0], 1, EPS, MaxRank, minn);

n=64;
BH = hcopy(AH);
CH = hmul(AH,BH);
@time CH = hmul(AH,BH);
#@profile CH = hmul(AH,BH);
#Profile.print()

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat