text:
using ModelicaScriptingTools
using Test: @testset, @test, @test_nowarn, @test_throws
using OMJulia: sendExpression
using DataFrames: Not, select!, DataFrame
using CSV
using Markdown
using Documenter

if !isdir("out")
    mkdir("out")
end

if !isdir("regRefData")
    mkdir("regRefData")
end

struct DummyPage
    source:: String
    workdir:: String
    mapping:: Dict
end

DummyPage(workdir) = DummyPage("dummydoc.md", workdir, Dict())

struct DummyInternal
    errors:: Array
end

struct DummyDocument
    internal:: DummyInternal
end

DummyDocument() = DummyDocument(DummyInternal([]))

@testset "MoST" begin
    @testset "moescape" begin
        @test raw"\\\"test\t\\data\?\\\"" == moescape("\"test\t\\data?\"")
    end
    @testset "mounescape" begin
        @test "\"test\t\\data?\"" == mounescape(raw"\\\"test\t\\data\?\\\"")
    end
    @testset "uniquehierarchy" begin
        @testset "simple" begin
            funcnames = [
                "a.b.c",
                "x.y.c",
                "ab.cd.ec",
                "ab.cd\$c"
            ]
            hier = uniquehierarchy(funcnames)
            expected = Dict(
                "a.b.c" => "b.c",
                "x.y.c" => "y.c",
                "ab.cd.ec" => "ec",
                "ab.cd\$c" => "cd.c"
            )
            @test expected == hier
        end
        @testset "string postfix != hierarchy postfix" begin
            funcnames = [
                "a.b.cd",
                "x.y.d"
            ]
            # these two names share a postfix (d) on the string level,
            # but they do not (cd vs d) on the hierarchical level
            # => last part of the hierarchy is sufficient for uniqueness
            hier = uniquehierarchy(funcnames)
            expected = Dict(
                "a.b.cd" => "cd",
                "x.y.d" => "d",
            )
            @test expected == hier
        end
    end
    @testset "replacefuncnames" begin
        input = """
        <math xmlns="http://www.w3.org/1998/Math/MathML">
        <mrow>
            <mrow><mi> foo</mi></mrow>
            <mo>≡</mo>
            <mrow>
                <mrow><mi>bla\$blubb.f</mi></mrow>
                <mo></mo>
                <mrow>
                    <mo>(</mo>
                    <mrow><mi> r</mi></mrow>
                    <mo>,</mo>
                    <mrow><mi> k</mi></mrow>
                    <mo>)</mo>
                </mrow>
            </mrow>
        </mrow>
        </math>"""
        replaced = replacefuncnames(input, Dict("bla\$blubb.f" => "f"))
        @test replace(input, "bla\$blubb.f" => "f") == replaced
    end
    @testset "Documenter.jl extension" begin
        @testset "DocExample" begin
            x = Markdown.parse("""
            ```@modelica
            %modeldir=res
            %libs=Modelica@3.2.3
            DocExample
            ```
            """).content[1]
            page = DummyPage(".")
            Documenter.Selectors.dispatch(
                Documenter.Expanders.ExpanderPipeline, x,
                page, DummyDocument()
            )
            result = page.mapping[x]
            @test result isa Documenter.Documents.MultiOutput
            @test length(result.content) == 6
            @test result.content[1] isa Markdown.Header{3}
            @test result.content[1].text == ["DocExample"]
            @test result.content[2] isa Documenter.Documents.RawHTML
            @test strip(result.content[2].code) == """

This is an example documentation for the DocExample class.
"""
            @test result.content[3] isa Markdown.Code
            @test replace(result.content[3].code, r"\s+" => "") == replace(read("res/DocExample.mo", String), r"\s+" => "")
            @test result.content[4] isa Documenter.Documents.RawHTML
            expected = """

<math xmlns="http://www.w3.org/1998/Math/MathML">
            <mrow><msup><mrow><mrow><mi> r
            </mi></mrow></mrow><mo>′</mo></msup><mo>≡</mo><mrow><mn> 1.0
            </mn><mo>/</mo><mrow><mrow><mi>g</mi></mrow><mo></mo><mrow><mo>(</mo><mrow><mi> foo
            </mi></mrow><mo>)</mo></mrow></mrow></mrow></mrow>
            </math>
<math xmlns="http://www.w3.org/1998/Math/MathML">
            <mrow><mrow><mi> foo
            </mi></mrow><mo>≡</mo><mrow><mrow><mi>f</mi></mrow><mo></mo><mrow><mo>(</mo><mrow><mi> r
            </mi></mrow><mo>,</mo><mrow><mi> k
            </mi></mrow><mo>)</mo></mrow></mrow></mrow>
            </math>
"""
            @test replace(result.content[4].code, r"\s+" => "") == replace(expected, r"\s+" => "")
            @test result.content[5] isa Markdown.MD
            @test occursin("```modelica\nfunction f", string(result.content[5]))
            @test occursin("```modelica\nfunction g", string(result.content[5]))
            variables14 = Markdown.parse("""
            | name | unit |                  label | value |
            | ----:| ----:| ----------------------:| -----:|
            |    r |  "V" |         some potential |   0.0 |
            |  foo |      | second sample variable |       |
            |    k |      |         some parameter |   2.0 |
            """)
            variables16 = Markdown.parse("""
            | name | unit |                  label | value |
            | ----:| ----:| ----------------------:| -----:|
            |  foo |      | second sample variable |       |
            |    r |  "V" |         some potential |   0.0 |
            |    k |      |         some parameter |   2.0 |
            """)
            @test result.content[6] isa Markdown.MD
            @test result.content[6] in [variables14, variables16]
        end
        @testset "Example (model without functions)" begin
            # only checks if models without functions pose errors
            x = Markdown.parse("""
            ```@modelica
            %modeldir=res
            %libs=Modelica@3.2.3
            Example
            ```
            """).content[1]
            page = DummyPage(".")
            Documenter.Selectors.dispatch(
                Documenter.Expanders.ExpanderPipeline, x,
                page, DummyDocument()
            )
            result = page.mapping[x]
        end
    end
    withOMC("out", "res") do omc
        installAndLoad(omc, "Modelica"; version="3.2.3")
        @testset "getVersion" begin
            major, minor, patch = getVersion(omc)
            # just test that version is sensible (and we have correct types)
            @test major >= 0
            @test minor >= 0
            @test patch >= 0
        end
        @testset "loadModel" begin
            mopath = sendExpression(omc, "getModelicaPath()")
            @testset "load existing correct model" begin
                @test_nowarn loadModel(omc, "Example")
            end
            @testset "load existing model with syntax error" begin
                modelfile = joinpath(pwd(), "res/SyntaxError.mo")
                try
                    loadModel(omc, "SyntaxError")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg == "Could not load SyntaxError"
                    @test occursin("Missing token: SEMICOLON", actual.omc)
                    @test occursin("$modelfile:6:3-6:3", actual.omc)
                end
            end
            @testset "load existing model with instantiation error" begin
                modelfile = joinpath(pwd(), "res/UndefinedVariable.mo")
                try
                    loadModel(omc, "UndefinedVariable")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg == "Model check of UndefinedVariable failed"
                    @test occursin("Error: Variable r not found in scope UndefinedVariable", actual.omc)
                    @test occursin("$modelfile:3:3-3:13", actual.omc)
                end
            end
            @testset "load non-existent model" begin
                try
                    loadModel(omc, "DoesNotExist")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg == "Could not load DoesNotExist"
                    @test occursin("Error: Failed to load package DoesNotExist (default)", actual.omc)
                end
            end
        end
        @testset "filename" begin
            @testset "Example" begin
                loadModel(omc, "Example")
                @test endswith(ModelicaScriptingTools.filename(omc, "Example"), "Example.mo")
            end
            @testset "<interactive>" begin
                sendExpression(omc, "model Foo end Foo;")
                @test "<interactive>" == ModelicaScriptingTools.filename(omc, "Foo")
            end
        end
        @testset "getSimulationSettings" begin
            loadModel(omc, "Example")
            @testset "read from model file" begin
                res = Dict(
                    "startTime" => 0.0, "stopTime" => 5.0,
                    "numberOfIntervals" => 50, "outputFormat" => "csv",
                    "variableFilter" => "sub\\.alias",
                    "tolerance" => 1.0e-6
                )
                @test res == getSimulationSettings(omc, "Example")
            end
            @testset "use override" begin
                @test 1 == getSimulationSettings(omc, "Example"; override=Dict("startTime" => 1))["startTime"]
            end
            @testset "nonexistant model" begin
                @test_throws MoSTError getSimulationSettings(omc, "DoesNotExist")
            end
        end
        @testset "simulate" begin
            mopath = sendExpression(omc, "getModelicaPath()")
            @testset "simulate correct model" begin
                loadModel(omc, "Example")
                @test_nowarn simulate(omc, "Example")
            end
            @testset "simulate model with arithmetic error" begin
                loadModel(omc, "ArithmeticError")
                try
                    simulate(omc, "ArithmeticError")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg == "Simulation of ArithmeticError failed"
                    @test occursin("Simulation execution failed for model: ArithmeticError", actual.omc)
                    @test occursin("division by zero at time 1.0", actual.omc)
                end
            end
            @testset "simulate model with initialization warning" begin
                loadModel(omc, "MissingInitialValue")
                try
                    simulate(omc, "MissingInitialValue")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg == "Simulation of MissingInitialValue failed"
                    @test occursin("Warning: The initial conditions are not fully specified.", actual.omc)
                end
            end
            @testset "simulate model with inconsistent units" begin
                try
                    # OpenModelica 1.16 error occurs here
                    loadModel(omc, "InconsistentUnits")
                    # OpenModelica 1.15 / 1.14 error occurs here
                    simulate(omc, "InconsistentUnits")
                catch actual
                    @test isa(actual, MoSTError)
                    @test actual.msg in [
                        "Simulation of InconsistentUnits failed",
                        "Model InconsistentUnits could not be instantiated"
                    ]
                    @test occursin("Warning: The following equation is INCONSISTENT", actual.omc)
                    @test occursin("sub.alias = r", actual.omc)
                    @test occursin("sub-expression \"r\" has unit \"A\"", actual.omc)
                    @test occursin("sub-expression \"sub.alias\" has unit \"V\"", actual.omc)
                end
            end
        end
        @testset "regressionTest" begin
            @testset "helper function resamplequi" begin
                testdata = DataFrame(
                    "time" => range(0, 1; length=10),
                    "x" => range(4, 10; length=10)
                )
                expected = DataFrame(
                    "time" => [0, 0.5, 1],
                    "x" => [testdata[1, "x"], testdata[6, "x"], testdata[10, "x"]]
                )
                @test expected == ModelicaScriptingTools.resamplequi(testdata, 3)
            end
            # we can only test correct regression test here
            @testset "regression tests of correct model" begin
                # setup simulation and reference data
                loadModel(omc, "Example")
                simulate(omc, "Example")
                cp("out/Example_res.csv", "regRefData/Example_res.csv"; force=true)
                regressionTest(omc, "Example", "regRefData"; relTol=1e-3, variableFilter="sub\\.alias")
            end
            @testset "regression test with missing alias in reference" begin
                # setup simulation and reference data
                loadModel(omc, "Example")
                simulate(omc, "Example")
                cp("out/Example_res.csv", "regRefData/Example_res.csv"; force=true)
                # remove column r from reference data, which should not have been selected in the first place
                csvfile = "regRefData/Example_res.csv"
                data = DataFrame(CSV.File(csvfile))
                select!(data, Not(:r))
                CSV.write(csvfile, data)
                regressionTest(omc, "Example", "regRefData"; relTol=1e-3, variableFilter="sub\\.alias")
            end
            @testset "regression test with ouputFormat=mat" begin
                # setup simulation and reference data
                loadModel(omc, "Example")
                simulate(omc, "Example", getSimulationSettings(omc, "Example"; override=Dict("outputFormat" => "mat")))
                cp("out/Example_res.mat", "regRefData/Example_res.mat"; force=true)
                regressionTest(omc, "Example", "regRefData"; relTol=1e-3, variableFilter="sub\\.alias", outputFormat="mat")
            end
        end
        @testset "testmodel" begin
            loadModel(omc, "Example")
            simulate(omc, "Example", getSimulationSettings(omc, "Example"))
            cp("out/Example_res.csv", "regRefData/Example_res.csv"; force=true)
            testmodel(omc, "Example", regRelTol=1e-4)
        end
        @testset "getDocAnnotation" begin
            loadModel(omc, "Example")
            expected = "\n
This is an example documentation for the Example class.
\n      "
            @test expected == getDocAnnotation(omc, "Example")
        end
        @testset "getcode" begin
            loadModel(omc, "Example.ExSub")
            expected = """model ExSub
              Modelica.SIunits.Voltage alias;
            end ExSub;"""
            @test expected == getcode(omc, "Example.ExSub")
        end
        @testset "getequations" begin
            loadModel(omc, "Example")
            expected = [
                "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n<mrow><msup><mrow><mrow><mi> r \n</mi></mrow></mrow><mo>′</mo></msup><mo>≡</mo><mn> 1.0 \n</mn></mrow>\n</math>"
            ]
            @test expected == getequations(omc, "Example")
        end
        @testset "deprefix" begin
            @testset "FunctionNames" begin
                loadModel(omc, "FunctionNames")
                vars = getvariables(omc, "FunctionNames")
                adict = ModelicaScriptingTools.aliasdict(vars)
                eqs = getequations(omc, "FunctionNames")
                funcs = getfunctions(omc, "FunctionNames")
                (funcdict, funcs) = uniquefunctions(funcs)
                eqs = [replacefuncnames(e, funcdict) for e in eqs]
                funcdict = uniquehierarchy(funcs[1:end, 1])
                eqs = [replacefuncnames(e, funcdict) for e in eqs]
                prefixes = [commonhierarchy(e, adict) for e in eqs]
                de = [deprefix(e, p) for (e, p) in zip(eqs, prefixes)]
                @test ["x", "f", "_b", "f", "_b"] == findidentifiers(de[1])
                @test ["_b", "f", "x", "g", "x"] == findidentifiers(de[2])
            end
        end
        @testset "getvariables" begin
            loadModel(omc, "Example")
            expected = [
                Dict(
                    "label" => "", "name" => "r",
                    "variability" => "continuousState", "unit" => "\"V\"",
                    "initial" => "0.0", "type" => "Real",
                    "quantity" => "\"ElectricPotential\"",
                    "bindExpression" => "", "aliasof" => ""
                ),
                Dict(
                    "label" => "", "name" => "sub.alias",
                    "variability" => "continuous", "unit" => "\"V\"",
                    "initial" => "", "type" => "Real",
                    "quantity" => "\"ElectricPotential\"",
                    "bindExpression" => "r", "aliasof" => "r"
                )
            ]
            @test expected == getvariables(omc, "Example")
        end
        @testset "getfunctions" begin
            @testset "DocExample" begin
                loadModel(omc, "DocExample")
                expected = [
                    "DocExample.f" "function(x :: Real * y :: Real) => Real" "function DocExample.f\n  input Real x;\n  input Real y;\n  output Real res;\nalgorithm\n  res := x ^ y + y;\nend DocExample.f;\n\n\n";
                    "DocExample.g" "function(x :: Real) => Real" "function DocExample.g\n  input Real x;\n  output Real y;\nalgorithm\n  y := 2.0 * x;\nend DocExample.g;\n\n\n"
                ]
                # seems like OpenModelica < 1.17.0 changes order of terms in equations
                if getVersion(omc) >= Tuple([1, 17, 0])
                    expected[2, 3] = replace(expected[2, 3], "y := 2.0 * x" => "y := x * 2.0")
                end
                @test expected == getfunctions(omc, "DocExample")
            end
            @testset "FunctionNames" begin
                loadModel(omc, "FunctionNames")
                funcs = getfunctions(omc, "FunctionNames")
                expected = [
                    "FunctionNames.f", "FunctionNames.f2", "FunctionNames.Submodel\$sm.f",
                    "FunctionNames.Submodel\$sm.g", "FunctionNames.Submodel\$sm.f.inf",
                    "FunctionNames.f.inf", "FunctionNames.f2.inf"
                ]
                # OpenModelica > 1.17.0 seems to drop the subclass names from
                # function names
                # Also it does no longer treat FunctionNames.sm.f.inf as
                # separate function.
                if getVersion(omc) >= Tuple([1, 17, 0])
                    expected = [
                        "FunctionNames.f", "FunctionNames.f2", "FunctionNames.f.inf",
                        "FunctionNames.f2.inf", "FunctionNames.sm.f", "FunctionNames.sm.g"
                    ]
                end
                @test expected == funcs[:, 1]
            end
        end
        @testset "funclist" begin
            loadModel(omc, "FunctionNames")
            funcs = getfunctions(omc, "FunctionNames")
            funclist = functionlist(funcs)
            expectedstr = """Functions:

            ```modelica
            function inf
              input Real x;
              input Real a = 1.0;
              output Real y;
            algorithm
              y := x + a;
            end inf;
            ```

            ```modelica
            function g
              input Real x;
              output Real y;
            algorithm
              y := 1.0 + x;
            end g;
            ```

            ```modelica
            function f
              input Real x1;
              input Real x2 = 1.0;
              output Real y;
            algorithm
              y := 2.0 * x1 + inf(x2, 1.0);
            end f;
            ```
            """
            if getVersion(omc) >= Tuple([1, 17, 0])
                expectedstr = replace(expectedstr, "y := 1.0 + x" => "y := x + 1.0")
                # FIXME OpenModelica 1.17.0-dev unfortunately uses x2 = 0.0,
                # which is wrong => if this is fixed, the following line must
                # be removed
                expectedstr = replace(expectedstr, "input Real x2 = 1.0" => "input Real x2 = 0.0")
            end
            expected = Markdown.parse(expectedstr)
            @test expected == funclist
        end
    end
end

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
export Mie

"""
    Mie(; m, n, cutoff, nl_only, lorentz_mixing, force_units, energy_units, skip_shortcut)

The Mie generalized interaction.
When `m` equals 6 and `n` equals 12 this is equivalent to the Lennard-Jones interaction.
"""
struct Mie{S, C, T, F, E} <: GeneralInteraction
    m::T
    n::T
    cutoff::C
    nl_only::Bool
    lorentz_mixing::Bool
    force_units::F
    energy_units::E
    mn_fac::T
end

function Mie(;
                m,
                n,
                cutoff=NoCutoff(),
                nl_only=false,
                lorentz_mixing=true,
                force_units=u"kJ * mol^-1 * nm^-1",
                energy_units=u"kJ * mol^-1",
                skip_shortcut=false)
    mn_fac = convert(typeof(m), (n / (n - m)) * (n / m) ^ (m / (n - m)))
    return Mie{skip_shortcut, typeof(cutoff), typeof(m), typeof(force_units), typeof(energy_units)}(
        m, n, cutoff, nl_only, lorentz_mixing, force_units, energy_units, mn_fac)
end

@fastmath @inbounds function force(inter::Mie{S, C, T},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C, T}
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)
    r = √r2

    if !S && iszero(atom_i.σ) || iszero(atom_j.σ)
        return ustrip.(zero(coord_i)) * inter.force_units
    end

    # Lorentz-Berthelot mixing rules use the arithmetic average for σ
    # Otherwise use the geometric average
    σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    m = inter.m
    n = inter.n
    # Derivative obtained via wolfram
    const_mn = inter.mn_fac * ϵ
    σ_r = σ / r
    params = (m, n, σ_r, const_mn)

    if cutoff_points(C) == 0
        f = force_divr_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        f = force_divr_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        if r2 < cutoff.activation_dist
            f = force_divr_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_divr_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_divr_nocutoff(::Mie, r2, invr2, (m, n, σ_r, const_mn))
    return -const_mn / r2 * (m * σ_r ^ m - n * σ_r ^ n)
end

@inline @inbounds function potential_energy(inter::Mie{S, C, T},
                                            coord_i,
                                            coord_j,
                                            atom_i,
                                            atom_j,
                                            box_size) where {S, C, T}
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)
    r = √r2

    if !S && iszero(atom_i.σ) || iszero(atom_j.σ)
        return ustrip(zero(box_size[1])) * inter.energy_units
    end

    σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    m = inter.m
    n = inter.n
    const_mn = inter.mn_fac * ϵ
    σ_r = σ / r
    params = (m, n, σ_r, const_mn)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        if r2 < cutoff.activation_dist
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::Mie, r2, invr2, (m, n, σ_r, const_mn))
    return const_mn * (σ_r ^ m - σ_r ^ m)
end

ext: jl


text:
module TestFileMgrPages
  using Base.Test
  using HelloDB.Pages
  using HelloDB.Blocks
  using HelloDB
  using TestHelpers: set_test_db

  function test_getint()
    page = Page(false)
    write(page.contents, 5, 2, 6)
    @test getint(page, 1 * sizeof(Int)) == 2
  end

  function test_setint()
    page = Page(false)
    write(page.contents, 0, 0, 0, 0)
    setint(page, 2, 1)
    seek(page.contents, 2)
    @test read(page.contents, Int) == 1
  end

  function test_getstring()
    page = Page(false)
    write(page.contents, 0x00, 5, 'P', 'e', 'n', 'i', 's', 0x00)
    @test getstring(page, 1) == "Penis"
  end

  function test_setstring()
    page = Page(false)
    write(page.contents, 0x00, 0x00)
    setstring(page, 2, "boobs")
    @test getstring(page, 2) == "boobs"
    @test page.contents.size == 2 + length("boobs") + sizeof(Int)
  end

  function test_read()
    set_test_db()
    filename = DB.filemgr.dbdir * "/test_read"

    f = open(filename, "w+")
    write(f, 2, 'h', 'i')
    close(f)

    block = Block(filename, 0)
    page = Page()
    read(page, block)
    @test getstring(page, 0) == "hi"
    dropdb(DB)
  end

  function test_write()
    set_test_db()
    filename = "test_write"

    block1 = Block(filename, 0)
    block2 = Block(filename, 1)
    block3 = Block(filename, 4)

    page1 = Page()
    setstring(page1, 0, "STRING ONE")

    page = Page()
    setstring(page, 0, "STRING TWO")

    write(page1, block1)
    write(page, block2)
    write(page1, block3)

    readpage1 = Page()
    readpage2 = Page()
    readpage3 = Page()

    read(readpage1, block1)
    read(readpage2, block2)
    read(readpage3, block3)

    @test getstring(readpage1, 0) == "STRING ONE"
    @test getstring(readpage2, 0) == "STRING TWO"
    @test getstring(readpage3, 0) == "STRING ONE"
    dropdb(DB)
  end

  function test_append()
    set_test_db()
    filename = "test_append"

    block1 = Block(filename, 0)
    block2 = Block(filename, 2)

    page1 = Page()
    setstring(page1, 0, "ABC")

    newpage = Page()
    setstring(newpage, 0, "DEF")

    write(page1, block1)
    write(page1, block2)
    newblock = append(newpage, filename)

    readpage1 = Page()
    readpage2 = Page()
    readpage3 = Page()

    read(readpage1, block1)
    read(readpage2, block2)
    read(readpage3, newblock)

    @test getstring(readpage1, 0) == "ABC"
    @test getstring(readpage2, 0) == "ABC"
    @test getstring(readpage3, 0) == "DEF"
    @test newblock.number == 3

    dropdb(DB)
  end

  @test string_size(2) == sizeof(Int) + (2) # TODO; not sizeof(Char) == 4

  test_getint()
  test_setint()
  test_getstring()
  test_setstring()
  test_read()
  test_write()
  test_append()
end


text:
module PlotsTests

include("imgcomp.jl")

# don't actually show the plots
srand(1234)
default(show=false, reuse=true)
img_eps = isinteractive() ? 1e-2 : 10e-2

# facts("Gadfly") do
#     @fact gadfly() --> Plots.GadflyBackend()
#     @fact backend() --> Plots.GadflyBackend()
#
#     @fact typeof(plot(1:10)) --> Plots.Plot{Plots.GadflyBackend}
#     @fact plot(Int[1,2,3], rand(3)) --> not(nothing)
#     @fact plot(sort(rand(10)), rand(Int, 10, 3)) --> not(nothing)
#     @fact plot!(rand(10,3), rand(10,3)) --> not(nothing)
#
#     image_comparison_facts(:gadfly, skip=[4,6,23,24,27], eps=img_eps)
# end

facts("PyPlot") do
    @fact pyplot() --> Plots.PyPlotBackend()
    @fact backend() --> Plots.PyPlotBackend()

    image_comparison_facts(:pyplot, skip=[25,30], eps=img_eps)
end

facts("GR") do
    @fact gr() --> Plots.GRBackend()
    @fact backend() --> Plots.GRBackend()

    if is_linux() && isinteractive()
        image_comparison_facts(:gr, skip=[2,25,30], eps=img_eps)
    end
end

facts("Plotly") do
    @fact plotly() --> Plots.PlotlyBackend()
    @fact backend() --> Plots.PlotlyBackend()

    # # until png generation is reliable on OSX, just test on linux
    # @static is_linux() && image_comparison_facts(:plotly, only=[1,3,4,7,8,9,10,11,12,14,15,20,22,23,27], eps=img_eps)
end

# facts("Immerse") do
#     @fact immerse() --> Plots.ImmerseBackend()
#     @fact backend() --> Plots.ImmerseBackend()
#
#     # as long as we can plot anything without error, it should be the same as Gadfly
#     image_comparison_facts(:immerse, only=[1], eps=img_eps)
# end

# facts("PlotlyJS") do
#     @fact plotlyjs() --> Plots.PlotlyJSBackend()
#     @fact backend() --> Plots.PlotlyJSBackend()
#
#     # as long as we can plot anything without error, it should be the same as Plotly
#     image_comparison_facts(:plotlyjs, only=[1], eps=img_eps)
# end

facts("UnicodePlots") do
    @fact unicodeplots() --> Plots.UnicodePlotsBackend()
    @fact backend() --> Plots.UnicodePlotsBackend()

    # lets just make sure it runs without error
    @fact isa(plot(rand(10)), Plot) --> true
end

facts("Axes") do
    p = plot()
    axis = p.subplots[1][:xaxis]
    @fact typeof(axis) --> Axis
    @fact Plots.discrete_value!(axis, "HI") --> (0.5, 1)
    @fact Plots.discrete_value!(axis, :yo) --> (1.5, 2)
    @fact extrema(axis) --> (0.5,1.5)
    @fact axis[:discrete_map] --> Dict{Any,Any}(:yo  => 2, "HI" => 1)

    Plots.discrete_value!(axis, ["x$i" for i=1:5])
    Plots.discrete_value!(axis, ["x$i" for i=0:2])
    @fact extrema(axis) --> (0.5, 7.5)
end

# tests for preprocessing recipes

# facts("recipes") do

    # user recipe

    # type T end
    # @recipe function f(::T)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     bg := :yellow
    #     rand(10)
    # end
    # plot(T())

    # plot recipe

    # @recipe function f(::Type{Val{:hiplt}},plt::Plot)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     t := :path
    #     bg:=:green
    #     ()
    # end
    # plot(rand(10),t=:hiplt)

    # series recipe

    # @recipe function f(::Type{Val{:hi}},x,y,z)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     t := :path
    #     ()
    # end
    # plot(rand(10),t=:hiplt)

# end

FactCheck.exitstatus()
end # module


text:
using LinearAlgebra
using Test

include("../src/threedQuadratureRule.jl")

#  Tests a quadrature formula for tetrahedron.  Rules must produce the
#  stated accuracy for monomials when compared to the exact integral
#  over the unit tetrahedron:
#
#      int(int(int( monomial, z,0,1-x-y), y,0,1-x), x,0,1)
#

  # Test degree 0 monomials
  function pass0(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,1)
    err[1] = sum(w)-1.0/6.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 1 monomials
  function pass1(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,3)
    err[1] = dot(r,w)-1.0/24.0
    err[2] = dot(s,w)-1.0/24.0
    err[3] = dot(t,w)-1.0/24.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 2 monomials
  function pass2(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,6)
    err[1] = dot(r.*r,w)-1.0/60.0
    err[2] = dot(r.*s,w)-1.0/120.0
    err[3] = dot(r.*t,w)-1.0/120.0
    err[4] = dot(s.*s,w)-1.0/60.0
    err[5] = dot(s.*t,w)-1.0/120.0
    err[6] = dot(t.*t,w)-1.0/60.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 3 monomials
  function pass3(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,10)
    err[ 1] = dot(r.*r.*r,w)-1.0/120.0
    err[ 2] = dot(r.*r.*s,w)-1.0/360.0
    err[ 3] = dot(r.*r.*t,w)-1.0/360.0
    err[ 4] = dot(r.*s.*s,w)-1.0/360.0
    err[ 5] = dot(r.*s.*t,w)-1.0/720.0
    err[ 6] = dot(r.*t.*t,w)-1.0/360.0
    err[ 7] = dot(s.*s.*s,w)-1.0/120.0
    err[ 8] = dot(s.*s.*t,w)-1.0/360.0
    err[ 9] = dot(s.*t.*t,w)-1.0/360.0
    err[10] = dot(t.*t.*t,w)-1.0/120.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 4 monomials
  function pass4(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,15)
    err[ 1] = dot(r.*r.*r.*r,w)-1.0/210.0
    err[ 2] = dot(r.*r.*r.*s,w)-1.0/840.0
    err[ 3] = dot(r.*r.*r.*t,w)-1.0/840.0
    err[ 4] = dot(r.*r.*s.*s,w)-1.0/1260.0
    err[ 5] = dot(r.*r.*s.*t,w)-1.0/2520.0
    err[ 6] = dot(r.*r.*t.*t,w)-1.0/1260.0
    err[ 7] = dot(r.*s.*s.*s,w)-1.0/840.0
    err[ 8] = dot(r.*s.*s.*t,w)-1.0/2520.0
    err[ 9] = dot(r.*s.*t.*t,w)-1.0/2520.0
    err[10] = dot(r.*t.*t.*t,w)-1.0/840.0
    err[11] = dot(s.*s.*s.*s,w)-1.0/210.0
    err[12] = dot(s.*s.*s.*t,w)-1.0/840.0
    err[13] = dot(s.*s.*t.*t,w)-1.0/1260.0
    err[14] = dot(s.*t.*t.*t,w)-1.0/840.0
    err[15] = dot(t.*t.*t.*t,w)-1.0/210.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 5 monomials
  function pass5(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,21)
    err[ 1] = dot(r.*r.*r.*r.*r,w)-1.0/336.0
    err[ 2] = dot(r.*r.*r.*r.*s,w)-1.0/1680.0
    err[ 3] = dot(r.*r.*r.*r.*t,w)-1.0/1680.0
    err[ 4] = dot(r.*r.*r.*s.*s,w)-1.0/3360.0
    err[ 5] = dot(r.*r.*r.*s.*t,w)-1.0/6720.0
    err[ 6] = dot(r.*r.*r.*t.*t,w)-1.0/3360.0
    err[ 7] = dot(r.*r.*s.*s.*s,w)-1.0/3360.0
    err[ 8] = dot(r.*r.*s.*s.*t,w)-1.0/10080.0
    err[ 9] = dot(r.*r.*s.*t.*t,w)-1.0/10080.0
    err[10] = dot(r.*r.*t.*t.*t,w)-1.0/3360.0
    err[11] = dot(r.*s.*s.*s.*s,w)-1.0/1680.0
    err[12] = dot(r.*s.*s.*s.*t,w)-1.0/6720.0
    err[13] = dot(r.*s.*s.*t.*t,w)-1.0/10080.0
    err[14] = dot(r.*s.*t.*t.*t,w)-1.0/6720.0
    err[15] = dot(r.*t.*t.*t.*t,w)-1.0/1680.0
    err[16] = dot(s.*s.*s.*s.*s,w)-1.0/336.0
    err[17] = dot(s.*s.*s.*s.*t,w)-1.0/1680.0
    err[18] = dot(s.*s.*s.*t.*t,w)-1.0/3360.0
    err[19] = dot(s.*s.*t.*t.*t,w)-1.0/3360.0
    err[20] = dot(s.*t.*t.*t.*t,w)-1.0/1680.0
    err[21] = dot(t.*t.*t.*t.*t,w)-1.0/336.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 6 monomials
  function pass6(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,28)
    err[ 1] = dot(r.*r.*r.*r.*r.*r,w)-1.0/504.0
    err[ 2] = dot(r.*r.*r.*r.*r.*s,w)-1.0/3024.0
    err[ 3] = dot(r.*r.*r.*r.*r.*t,w)-1.0/3024.0
    err[ 4] = dot(r.*r.*r.*r.*s.*s,w)-1.0/7560.0
    err[ 5] = dot(r.*r.*r.*r.*s.*t,w)-1.0/15120.0
    err[ 6] = dot(r.*r.*r.*r.*t.*t,w)-1.0/7560.0
    err[ 7] = dot(r.*r.*r.*s.*s.*s,w)-1.0/10080.0
    err[ 8] = dot(r.*r.*r.*s.*s.*t,w)-1.0/30240.0
    err[ 9] = dot(r.*r.*r.*s.*t.*t,w)-1.0/30240.0
    err[10] = dot(r.*r.*r.*t.*t.*t,w)-1.0/10080.0
    err[11] = dot(r.*r.*s.*s.*s.*s,w)-1.0/7560.0
    err[12] = dot(r.*r.*s.*s.*s.*t,w)-1.0/30240.0
    err[13] = dot(r.*r.*s.*s.*t.*t,w)-1.0/45360.0
    err[14] = dot(r.*r.*s.*t.*t.*t,w)-1.0/30240.0
    err[15] = dot(r.*r.*t.*t.*t.*t,w)-1.0/7560.0
    err[16] = dot(r.*s.*s.*s.*s.*s,w)-1.0/3024.0
    err[17] = dot(r.*s.*s.*s.*s.*t,w)-1.0/15120.0
    err[18] = dot(r.*s.*s.*s.*t.*t,w)-1.0/30240.0
    err[19] = dot(r.*s.*s.*t.*t.*t,w)-1.0/30240.0
    err[20] = dot(r.*s.*t.*t.*t.*t,w)-1.0/15120.0
    err[21] = dot(r.*t.*t.*t.*t.*t,w)-1.0/3024.0
    err[22] = dot(s.*s.*s.*s.*s.*s,w)-1.0/504.0
    err[23] = dot(s.*s.*s.*s.*s.*t,w)-1.0/3024.0
    err[24] = dot(s.*s.*s.*s.*t.*t,w)-1.0/7560.0
    err[25] = dot(s.*s.*s.*t.*t.*t,w)-1.0/10080.0
    err[26] = dot(s.*s.*t.*t.*t.*t,w)-1.0/7560.0
    err[27] = dot(s.*t.*t.*t.*t.*t,w)-1.0/3024.0
    err[28] = dot(t.*t.*t.*t.*t.*t,w)-1.0/504.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 7 monomials
  function pass7(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,36)
    err[ 1] = dot(r.*r.*r.*r.*r.*r.*r,w)-1.0/720.0
    err[ 2] = dot(r.*r.*r.*r.*r.*r.*s,w)-1.0/5040.0
    err[ 3] = dot(r.*r.*r.*r.*r.*r.*t,w)-1.0/5040.0
    err[ 4] = dot(r.*r.*r.*r.*r.*s.*s,w)-1.0/15120.0
    err[ 5] = dot(r.*r.*r.*r.*r.*s.*t,w)-1.0/30240.0
    err[ 6] = dot(r.*r.*r.*r.*r.*t.*t,w)-1.0/15120.0
    err[ 7] = dot(r.*r.*r.*r.*s.*s.*s,w)-1.0/25200.0
    err[ 8] = dot(r.*r.*r.*r.*s.*s.*t,w)-1.0/75600.0
    err[ 9] = dot(r.*r.*r.*r.*s.*t.*t,w)-1.0/75600.0
    err[10] = dot(r.*r.*r.*r.*t.*t.*t,w)-1.0/25200.0
    err[11] = dot(r.*r.*r.*s.*s.*s.*s,w)-1.0/25200.0
    err[12] = dot(r.*r.*r.*s.*s.*s.*t,w)-1.0/100800.0
    err[13] = dot(r.*r.*r.*s.*s.*t.*t,w)-1.0/151200.0
    err[14] = dot(r.*r.*r.*s.*t.*t.*t,w)-1.0/100800.0
    err[15] = dot(r.*r.*r.*t.*t.*t.*t,w)-1.0/25200.0
    err[16] = dot(r.*r.*s.*s.*s.*s.*s,w)-1.0/15120.0
    err[17] = dot(r.*r.*s.*s.*s.*s.*t,w)-1.0/75600.0
    err[18] = dot(r.*r.*s.*s.*s.*t.*t,w)-1.0/151200.0
    err[19] = dot(r.*r.*s.*s.*t.*t.*t,w)-1.0/151200.0
    err[20] = dot(r.*r.*s.*t.*t.*t.*t,w)-1.0/75600.0
    err[21] = dot(r.*r.*t.*t.*t.*t.*t,w)-1.0/15120.0
    err[22] = dot(r.*s.*s.*s.*s.*s.*s,w)-1.0/5040.0
    err[23] = dot(r.*s.*s.*s.*s.*s.*t,w)-1.0/30240.0
    err[24] = dot(r.*s.*s.*s.*s.*t.*t,w)-1.0/75600.0
    err[25] = dot(r.*s.*s.*s.*t.*t.*t,w)-1.0/100800.0
    err[26] = dot(r.*s.*s.*t.*t.*t.*t,w)-1.0/75600.0
    err[27] = dot(r.*s.*t.*t.*t.*t.*t,w)-1.0/30240.0
    err[28] = dot(r.*t.*t.*t.*t.*t.*t,w)-1.0/5040.0
    err[29] = dot(s.*s.*s.*s.*s.*s.*s,w)-1.0/720.0
    err[30] = dot(s.*s.*s.*s.*s.*s.*t,w)-1.0/5040.0
    err[31] = dot(s.*s.*s.*s.*s.*t.*t,w)-1.0/15120.0
    err[32] = dot(s.*s.*s.*s.*t.*t.*t,w)-1.0/25200.0
    err[33] = dot(s.*s.*s.*t.*t.*t.*t,w)-1.0/25200.0
    err[34] = dot(s.*s.*t.*t.*t.*t.*t,w)-1.0/15120.0
    err[35] = dot(s.*t.*t.*t.*t.*t.*t,w)-1.0/5040.0
    err[36] = dot(t.*t.*t.*t.*t.*t.*t,w)-1.0/720.0
    return norm(err,Inf)<1e-15
  end

  # Test degree 8 monomials
  function pass8(r::Vector{Float64},s::Vector{Float64},t::Vector{Float64},w::Vector{Float64})
    err = zeros(Float64,45)
    err[ 1] = dot(r.*r.*r.*r.*r.*r.*r.*r,w)-1.0/990.0
    err[ 2] = dot(r.*r.*r.*r.*r.*r.*r.*s,w)-1.0/7920.0
    err[ 3] = dot(r.*r.*r.*r.*r.*r.*r.*t,w)-1.0/7920.0
    err[ 4] = dot(r.*r.*r.*r.*r.*r.*s.*s,w)-1.0/27720.0
    err[ 5] = dot(r.*r.*r.*r.*r.*r.*s.*t,w)-1.0/55440.0
    err[ 6] = dot(r.*r.*r.*r.*r.*r.*t.*t,w)-1.0/27720.0
    err[ 7] = dot(r.*r.*r.*r.*r.*s.*s.*s,w)-1.0/55440.0
    err[ 8] = dot(r.*r.*r.*r.*r.*s.*s.*t,w)-1.0/166320.0
    err[ 9] = dot(r.*r.*r.*r.*r.*s.*t.*t,w)-1.0/166320.0
    err[10] = dot(r.*r.*r.*r.*r.*t.*t.*t,w)-1.0/55440.0
    err[11] = dot(r.*r.*r.*r.*s.*s.*s.*s,w)-1.0/69300.0
    err[12] = dot(r.*r.*r.*r.*s.*s.*s.*t,w)-1.0/277200.0
    err[13] = dot(r.*r.*r.*r.*s.*s.*t.*t,w)-1.0/415800.0
    err[14] = dot(r.*r.*r.*r.*s.*t.*t.*t,w)-1.0/277200.0
    err[15] = dot(r.*r.*r.*r.*t.*t.*t.*t,w)-1.0/69300.0
    err[16] = dot(r.*r.*r.*s.*s.*s.*s.*s,w)-1.0/55440.0
    err[17] = dot(r.*r.*r.*s.*s.*s.*s.*t,w)-1.0/277200.0
    err[18] = dot(r.*r.*r.*s.*s.*s.*t.*t,w)-1.0/554400.0
    err[19] = dot(r.*r.*r.*s.*s.*t.*t.*t,w)-1.0/554400.0
    err[20] = dot(r.*r.*r.*s.*t.*t.*t.*t,w)-1.0/277200.0
    err[21] = dot(r.*r.*r.*t.*t.*t.*t.*t,w)-1.0/55440.0
    err[22] = dot(r.*r.*s.*s.*s.*s.*s.*s,w)-1.0/27720.0
    err[23] = dot(r.*r.*s.*s.*s.*s.*s.*t,w)-1.0/166320.0
    err[24] = dot(r.*r.*s.*s.*s.*s.*t.*t,w)-1.0/415800.0
    err[25] = dot(r.*r.*s.*s.*s.*t.*t.*t,w)-1.0/554400.0
    err[26] = dot(r.*r.*s.*s.*t.*t.*t.*t,w)-1.0/415800.0
    err[27] = dot(r.*r.*s.*t.*t.*t.*t.*t,w)-1.0/166320.0
    err[28] = dot(r.*r.*t.*t.*t.*t.*t.*t,w)-1.0/27720.0
    err[29] = dot(r.*s.*s.*s.*s.*s.*s.*s,w)-1.0/7920.0
    err[30] = dot(r.*s.*s.*s.*s.*s.*s.*t,w)-1.0/55440.0
    err[31] = dot(r.*s.*s.*s.*s.*s.*t.*t,w)-1.0/166320.0
    err[32] = dot(r.*s.*s.*s.*s.*t.*t.*t,w)-1.0/277200.0
    err[33] = dot(r.*s.*s.*s.*t.*t.*t.*t,w)-1.0/277200.0
    err[34] = dot(r.*s.*s.*t.*t.*t.*t.*t,w)-1.0/166320.0
    err[35] = dot(r.*s.*t.*t.*t.*t.*t.*t,w)-1.0/55440.0
    err[36] = dot(r.*t.*t.*t.*t.*t.*t.*t,w)-1.0/7920.0
    err[37] = dot(s.*s.*s.*s.*s.*s.*s.*s,w)-1.0/990.0
    err[38] = dot(s.*s.*s.*s.*s.*s.*s.*t,w)-1.0/7920.0
    err[39] = dot(s.*s.*s.*s.*s.*s.*t.*t,w)-1.0/27720.0
    err[40] = dot(s.*s.*s.*s.*s.*t.*t.*t,w)-1.0/55440.0
    err[41] = dot(s.*s.*s.*s.*t.*t.*t.*t,w)-1.0/69300.0
    err[42] = dot(s.*s.*s.*t.*t.*t.*t.*t,w)-1.0/55440.0
    err[43] = dot(s.*s.*t.*t.*t.*t.*t.*t,w)-1.0/27720.0
    err[44] = dot(s.*t.*t.*t.*t.*t.*t.*t,w)-1.0/7920.0
    err[45] = dot(t.*t.*t.*t.*t.*t.*t.*t,w)-1.0/990.0
    return norm(err,Inf)<1e-15
  end

  r,s,t,w = threedQuadratureRule(1)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)

#  r,s,t,w = threedQuadratureRule(4)
#  @test pass0(r,s,t,w)
#  @test pass1(r,s,t,w)
#  @test pass2(r,s,t,w)

  r,s,t,w = threedQuadratureRule(5)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)

  r,s,t,w = threedQuadratureRule(11)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)
  @test pass4(r,s,t,w)

  r,s,t,w = threedQuadratureRule(15)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)
  @test pass4(r,s,t,w)
  @test pass5(r,s,t,w)

  r,s,t,w = threedQuadratureRule(24)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)
  @test pass4(r,s,t,w)
  @test pass5(r,s,t,w)
  @test pass6(r,s,t,w)

  r,s,t,w = threedQuadratureRule(31)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)
  @test pass4(r,s,t,w)
  @test pass5(r,s,t,w)
  @test pass6(r,s,t,w)
  @test pass7(r,s,t,w)

  r,s,t,w = threedQuadratureRule(45)
  @test pass0(r,s,t,w)
  @test pass1(r,s,t,w)
  @test pass2(r,s,t,w)
  @test pass3(r,s,t,w)
  @test pass4(r,s,t,w)
  @test pass5(r,s,t,w)
  @test pass6(r,s,t,w)
  @test pass7(r,s,t,w)
  @test pass8(r,s,t,w)


text:
# (A D⁻¹ Aᵀ + δI) Δy = A D⁻¹ ξ₁ + ξ₂
# where D = s_l (x - lvar)⁻¹ + s_u (uvar - x)⁻¹ + ρI,
# and the right hand side of K2 is rhs = [ξ₁]
#                                        [ξ₂]
export K1KrylovParams

"""
Type to use the K1 formulation with a Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
The outer constructor

    K1KrylovParams(; uplo = :L, kmethod = :cg, preconditioner = :Identity,
                   ratol = 1.0e-10, rrtol = 1.0e-10)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:cg`
- `:minres`
- `:minres_qlp`

"""
struct K1KrylovParams <: SolverParams
  uplo::Symbol
  kmethod::Symbol
  preconditioner::Symbol
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ_min::Float64
  δ_min::Float64
end

function K1KrylovParams(;
  uplo::Symbol = :L,
  kmethod::Symbol = :cg,
  preconditioner::Symbol = :Identity,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ_min::T = 1e3 * sqrt(eps()),
  δ_min::T = 1e4 * sqrt(eps()),
) where {T <: Real}
  return K1KrylovParams(
    uplo,
    kmethod,
    preconditioner,
    atol0,
    rtol0,
    atol_min,
    rtol_min,
    ρ_min,
    δ_min,
  )
end

mutable struct PreallocatedDataK1Krylov{T <: Real, S, L <: LinearOperator, Ksol <: KrylovSolver} <:
               PreallocatedDataNormalKrylov{T, S}
  D::S
  rhs::S
  regu::Regularization{T}
  δv::Vector{T}
  K::L # augmented matrix (LinearOperator)
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function opK1prod!(
  res::AbstractVector{T},
  D::AbstractVector{T},
  A::Union{AbstractMatrix{T}, AbstractLinearOperator{T}},
  δv::AbstractVector{T},
  v::AbstractVector{T},
  vtmp::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  if uplo == :U
    mul!(vtmp, A, v)
    mul!(res, A', vtmp ./ D, α, β)
    res .+= (α * δv[1]) .* v
  else
    mul!(vtmp, A', v)
    mul!(res, A, vtmp ./ D, α, β)
    res .+= (α * δv[1]) .* v
  end
end

function PreallocatedData(
  sp::K1KrylovParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}
  D = similar(fd.c, id.nvar)
  # init Regularization values
  if iconf.mode == :mono
    regu =
      Regularization(T(sqrt(eps()) * 1e5), T(sqrt(eps()) * 1e5), T(sp.ρ_min), T(sp.δ_min), :classic)
    # Regularization(T(0.), T(0.), T(sp.ρ_min), T(sp.δ_min), :classic)
    D .= T(1.0e0) / 2
  else
    regu = Regularization(
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps(T)) * 1e0),
      T(sqrt(eps(T)) * 1e0),
      :classic,
    )
    D .= T(1.0e-2)
  end

  δv = [regu.δ] # put it in a Vector so that we can modify it without modifying opK2prod!
  K = LinearOperator(
    T,
    id.ncon,
    id.ncon,
    true,
    true,
    (res, v, α, β) -> opK1prod!(res, D, fd.A, δv, v, similar(fd.c), α, β, fd.uplo),
  )

  rhs = similar(fd.c, id.ncon)
  kstring = string(sp.kmethod)
  KS = eval(KSolver(sp.kmethod))(K, rhs)

  return PreallocatedDataK1Krylov(
    D,
    rhs,
    regu,
    δv,
    K, #K
    KS,
    sp.atol0,
    sp.rtol0,
    sp.atol_min,
    sp.rtol_min,
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK1Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  pad.rhs .= @views dd[(id.nvar + 1):end]
  if fd.uplo == :U
    @views mul!(pad.rhs, fd.A', dd[1:(id.nvar)] ./ pad.D, one(T), one(T))
  else
    @views mul!(pad.rhs, fd.A, dd[1:(id.nvar)] ./ pad.D, one(T), one(T))
  end
  rhsNorm = kscale!(pad.rhs)
  pad.K.nprod = 0
  ksolve!(pad.KS, pad.K, pad.rhs, I(id.ncon), verbose = 0, atol = pad.atol, rtol = pad.rtol)
  update_kresiduals_history!(res, pad.K, pad.KS.x, pad.rhs)
  kunscale!(pad.KS.x, rhsNorm)

  if fd.uplo == :U
    @views mul!(dd[1:(id.nvar)], fd.A, pad.KS.x, one(T), -one(T))
  else
    @views mul!(dd[1:(id.nvar)], fd.A', pad.KS.x, one(T), -one(T))
  end
  dd[1:(id.nvar)] ./= pad.D
  dd[(id.nvar + 1):end] .= pad.KS.x
  return 0
end

function update_pad!(
  pad::PreallocatedDataK1Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  pad.δv[1] = pad.regu.δ

  pad.D .= pad.regu.ρ
  pad.D[id.ilow] .+= pt.s_l ./ itd.x_m_lvar
  pad.D[id.iupp] .+= pt.s_u ./ itd.uvar_m_x
  pad.δv[1] = pad.regu.δ

  # update_preconditioner!(pad.pdat, pad, itd, pt, id, fd, cnts)

  return 0
end


text:
include("1.jl")
include("11.jl")
function mark_krossper(r::Robot)
    for side in (Nord,Sud,Ost,West)
        putmarkholeper(r,side)
        backingper(r, inverse(side))
    end
    putmarker!(r)
end

function putmarkholeper(r::Robot,side::HorizonSide)
    y = 0
    while (y != 1)
        y = disexit(r,side)
        putmarker!(r)
    end
end

function backingper(r::Robot,side::HorizonSide)
    while ismarker(r)
        disexit(r,side)
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
__precompile__()

"""
This module orchestrates a GWAS analysis.
"""
module MendelGWAS
#
# Required OpenMendel packages and modules.
#
using MendelBase
using MendelPlots
using SnpArrays
#
# Required external modules.
#
using CSV
using DataFrames
using Distributions
using GLM
using LinearAlgebra
using Missings
using Printf
using StatsBase
using StatsModels

export GWAS

"""
This is the wrapper function for the GWAS analysis option.
"""
function GWAS(control_file = ""; args...)
  #
  # Print the logo. Store the initial directory.
  #
  print(" \n \n")
  println("     Welcome to OpenMendel's")
  println("      GWAS analysis option")
  print(" \n \n")
  println("Reading the data.\n")
  initial_directory = pwd()
  #
  # The user specifies the analysis to perform via a set of keywords.
  # Start the keywords at their default values.
  #
  keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
  #
  # Keywords unique to this analysis should be first defined here
  # by setting their default values using the format:
  # keyword["some_keyword_name"] = default_value
  #
  keyword["distribution"] = "" # Binomial(), Gamma(), Normal(), Poisson(), etc.
  keyword["link"] = ""         # LogitLink(), IdentityLink(), LogLink(), etc.
  keyword["lrt_threshold"] = -1.0 # Default changed below to 0.05 / SNPs
  keyword["maf_threshold"] = 0.01
  keyword["min_success_rate_per_sample"] = 0.98
  keyword["min_success_rate_per_snp"] = 0.98
  keyword["manhattan_plot_file"] = ""
  keyword["output_table"] = "" # Table of all SNPs and their p-values
  keyword["qq_plot_file"] = ""
  keyword["regression"] = ""   # linear, logistic, or poisson
  keyword["regression_formula"] = ""
  keyword["signif_threshold"] = -1.0 # Default changed below to 0.05 / SNPs
  keyword["pcs"] = 0      # Number of Principal Components to include in model.
  #
  # Process the run-time user-specified keywords that will control the analysis.
  # This will also initialize the random number generator.
  #
  process_keywords!(keyword, control_file, args)
  #
  # Check that the correct analysis option was specified.
  #
  lc_analysis_option = lowercase(keyword["analysis_option"])
  if (lc_analysis_option != "" &&
      lc_analysis_option != "gwas")
     throw(ArgumentError(
       "An incorrect analysis option was specified.\n \n"))
  end
  keyword["analysis_option"] = "GWAS"
  #
  # Read the genetic data from the external files named in the keywords.
  #
  (pedigree, person, nuclear_family, locus, snpdata,
    locus_frame, phenotype_frame, person_frame, snp_definition_frame) =
    read_external_data_files(keyword)
  #
  # Check if SNP data were read.
  #
  if snpdata.snps == 0
    println(" \n\nERROR: This analysis requires SNP data and none were read!\n")
  else
    #
    # If the Significance Threshold has not been set by the user,
    # set it to the Bonferroni correction: 0.05 / SNPs.
    # If the LRT Threshold has not been set by the user,
    # set it to the Significance Threshold.
    #
    if keyword["signif_threshold"] < 0.0
      keyword["signif_threshold"] = 0.05 / snpdata.snps
    end
    if keyword["lrt_threshold"] < 0.0
      keyword["lrt_threshold"] = keyword["signif_threshold"]
    end
    #
    # If principal components are requested to be included in the model,
    # then call a function that will use the PCA routine in SnpArrays,
    # and will add these PCs to the person_frame. Once in the pedigree frame
    # these PCs will be included in the regression procedure.
    # Each PC will be named PCn where n is its number, e.g., "PC1".
    #
##!!  if keyword["pcs"] > 0
##!!    add_pcs!(person_frame, snpdata, keyword)
##!!  end
    #
    # Execute the specified analysis.
    #
    println(" \nAnalyzing the data.\n")
    execution_error = gwas_option(person, snpdata, person_frame, keyword)
    if execution_error
      println(" \n \nERROR: Mendel terminated prematurely!\n")
    else
      println(" \n \nMendel's analysis is finished.\n")
    end
  end
  #
  # Finish up by closing, and thus flushing, any output files.
  # Return to the initial directory.
  #
  close(keyword["output_unit"])
  cd(initial_directory)
  return nothing
end # function GWAS

"""
This function performs GWAS on a set of traits.
"""
function gwas_option(person::Person, snpdata::SnpDataStruct,
  person_frame::DataFrame, keyword::Dict{AbstractString, Any})

  TAB_CHAR :: Char = Char(9)

  people = person.people
  snps = snpdata.snps
  io = keyword["output_unit"]
  lrt_threshold = keyword["lrt_threshold"]
  maf_threshold = keyword["maf_threshold"]
  signif_threshold = keyword["signif_threshold"]
  min_success_rate_per_sample = keyword["min_success_rate_per_sample"]
  min_success_rate_per_snp = keyword["min_success_rate_per_snp"]
  #
  # Recognize the three basic GWAS regression models:
  # linear, logistic, and poisson.
  # For these three we will use fast internal regression code,
  # unless an interaction term is detected.
  #
  regression_type = lowercase(keyword["regression"])
  distribution_family = keyword["distribution"]
  link = keyword["link"]
  if regression_type == "linear" ||
     (distribution_family == Normal() && link == IdentityLink())
    keyword["distribution"] = Normal()
    keyword["link"] = IdentityLink()
    keyword["regression"] = "linear"
    fast_method = true
  elseif regression_type == "logistic" ||
     (distribution_family == Binomial() && link == LogitLink())
    keyword["distribution"] = Binomial()
    keyword["link"] = LogitLink()
    keyword["regression"] = "logistic"
    fast_method = true
  elseif regression_type == "poisson" ||
     (distribution_family == Poisson() && link == LogLink())
    keyword["distribution"] = Poisson()
    keyword["link"] = LogLink()
    keyword["regression"] = "poisson"
    fast_method = true
  elseif regression_type == "" && distribution_family == "" && link == ""
    throw(ArgumentError(
      "No regression type has been defined for this analysis.\n" *
      "Set the keyword regression to either:\n" *
      "  linear (for usual quantitative traits),\n" *
      "  logistic (for usual qualitative traits),\n" *
      "  poisson, or blank (for unusual analyses).\n" *
      "If keyword regression is not assigned, then the keywords\n" *
      "'distribution' and 'link' must be assigned names of functions.\n \n"))
  else
    fast_method = false
  end
  regression_type = lowercase(keyword["regression"])
  distribution_family = keyword["distribution"]
  link = keyword["link"]
  if regression_type != "" && regression_type != "linear" &&
     regression_type != "logistic" && regression_type != "poisson"
    throw(ArgumentError(
     "The keyword regression (currently: $regression_type) must be assigned\n" *
     "the value 'linear', 'logistic', 'poisson', or blank.\n \n"))
  end
  #
  # Retrieve and error check the regression formula.
  #
  regression_formula = keyword["regression_formula"]
  if regression_formula == ""
    throw(ArgumentError(
      "The keyword regression_formula appears blank.\n" *
      "A regression formula must be provided.\n \n"))
  end
  if !occursin("~", regression_formula)
    throw(ArgumentError(
      "The value of the keyword regression_formula ('$regression_formula')\n" *
      "does not contain the required '~' that should separate\n" *
      "the trait and the predictors.\n \n"))
  end
  side = split(regression_formula, "~")
  side[1] = strip(side[1])
  side[2] = strip(side[2])
  if side[1] == ""
    throw(ArgumentError(
      "The left hand side of the formula specified in\n" *
      "the keyword regression_formula appears blank.\n" *
      "The left hand side should contain only the name of the trait field\n" *
      "in the Pedigree data. If the trait field is unnamed, use 'Trait'.\n \n"))
  end
  entry_separators = [' ', TAB_CHAR, ',', ';', '+', '*', '&']
  if something(findfirst((in)(entry_separators), side[1]), 0) > 0
    lhs_string = side[1]
    throw(ArgumentError(
      "The left hand side ('$lhs_string') of the formula specified in\n" *
      "the keyword regression_formula appears to have multiple entries.\n" *
      "The left hand side should contain only the name of the trait field\n" *
      "in the Pedigree data. If the trait field is unnamed, use 'Trait'.\n \n"))
  end
  if side[2] == "" || side[2] == "0" || side[2] == "-1"
    println("WARNING: Since the right hand side of the provided")
    println("regression formula ('$regression_formula') is empty")
    println("(that is, has neither intercept nor predictor),")
    println("the intercept has been added to the regression.\n ")
    side[2] = "1"
  end
  lhs = Meta.parse(side[1])
  rhs = Meta.parse(side[2])
  if !(lhs in names(person_frame))
    lhs_string = string(lhs)
    throw(ArgumentError(
      "The field named on the left hand side of the formula specified in\n" *
      "the keyword regression_formula (currently: '$lhs_string')\n" *
      "is not in the Pedigree data file.\n" *
      "The left hand side should contain only the name of the trait field\n" *
      "in the Pedigree data. If the trait field is unnamed, use 'Trait'.\n \n"))
  end
  #
  # If the regression formula includes an interaction term,
  # do not use the fast internal regression code.
  #
  interaction_symbols = ['*', '&']
  if something(findfirst((in)(interaction_symbols), string(rhs)), 0) > 0
    fast_method = false
  end
  #
  # Create a list of the names of the columns used in the formula.
  # NB: Since this procedure only looks for string matches
  # it will be flummoxed by, for example, a column named "l",
  # if the "log" function is used in the formula.
  #
  persondf_colname = names(person_frame)
  persondf_colnames = size(persondf_colname, 1)
  indicat = falses(persondf_colnames)
  for i = 1:persondf_colnames
    namestring = string(persondf_colname[i])
    if namestring == string(lhs) ||
       first(something(findfirst(namestring, string(rhs)), 0:-1)) > 0
      indicat[i] = true
    end
  end
  formula_colname = fill(:placeholder, sum(indicat))
  j = 0
  for i = 1:persondf_colnames
    if indicat[i]
      j = j + 1
      formula_colname[j] = persondf_colname[i]
    end
  end
  if j != sum(indicat)
    throw(ArgumentError("Inconsistency in formula column count./n /n"))
  end
  #
  # Change sex designations to 1.0 (for females) and -1.0 (for males).
  # Since the field :Sex may have type String,
  # create a new field of type Float64 that replaces :Sex.
  # Note that unrecognized sexes are labeled female.
  #
  if first(something(findfirst("Sex", string(rhs)), 0:-1)) > 0 &&
     in(:Sex, names(person_frame))
    person_frame[!, :NumericSex] = ones(people)
    for i = 1:people
      s = person_frame[i, :Sex]
      if !isa(Meta.parse(string(s), raise=false), Number); s = lowercase(s); end
      if !(s in keyword["female"]); person_frame[i, :NumericSex] = -1.0; end
    end
    names_list = names(person_frame)
    deleteat!(names_list, findall((in)([:Sex]), names_list))
    select!(person_frame, names_list)
    rename!(person_frame, :NumericSex => :Sex)
  end
  #
  # For Logistic regression make sure the cases are 1.0,
  # non-cases are 0.0, and missing data is "missing".
  # Again, since the trait field may be of type String,
  # create a new field of type Float64 that replaces it.
  #
  case_label = keyword["affected_designator"]
  if regression_type == "logistic" && case_label != ""
    person_frame[!, :NumericTrait] = zeros(people)
    for i = 1:people
      s = strip(string(person_frame[i, lhs]))
      if s == "" || s == "NaN" || s == "missing" || s == "NA"
        person_frame[i, :NumericTrait] = missing
      elseif s == case_label
        person_frame[i, :NumericTrait] = 1.0
      end
    end
    names_list = names(person_frame)
    deleteat!(names_list, findall((in)([lhs]), names_list))
    select!(person_frame, names_list)
    rename!(person_frame, :NumericTrait => lhs)
  end
  #
  # To filter the SNP data, first find the SNPs and samples
  # that surpass the requested success rates and MAF.
  # Note that the sample_mask is on the people in entry order!
  #
  sample_mask, snp_mask = SnpArrays.filter(snpdata.snpmatrix,
    min_success_rate_per_row = min_success_rate_per_snp,
    min_success_rate_per_col = min_success_rate_per_sample,
    min_maf = maf_threshold)
  #
  # Create a model frame, reordering the individuals (i.e., the rows)
  # by their entry order to align with the corresponding SNP data.
  #
  model_frame = deepcopy(person_frame)
  sort!(model_frame, :EntryOrder)
  #
  # Restrict the model_frame to the columns used in the regression formula.
  # (Perhaps wait on this until we know it's working well.)
  #
##  select!(model_frame, formula_colname)
  #
  # Remove the individuals with low genotyping success rates
  # using the sample_mask from SnpArrays.
  # Note that now rows = people - too.few.genotypes.
  #
  model_frame = model_frame[sample_mask, :]
  #
  # Next, drop all individuals with missing values among the variables
  # used in the regression formula. Keep track of the corresponding mask
  # to use with the SNP data. Thus, model_frame has no missing data.
  # Note that the completeness_mask should be applied after the sample_mask.
  # Note that now rows = people - too.few.genotypes - incomplete.predictors.
  #
  completeness_mask = completecases(model_frame, formula_colname)
  dropmissing!(model_frame, formula_colname, disallowmissing=true)
  cases = size(model_frame, 1) # also = sum(completeness_mask)
##
##   model = ModelFrame(fm, person_frame[sample_mask, :])
##   #
##   # To ensure that the trait and SNPs occur in the same order, sort the
##   # the model dataframe by the entry-order column of the pedigree frame.
##   #
##   model.df[!, :EntryOrder] = person_frame[sample_mask, :EntryOrder]
##   sort!(model.df, :EntryOrder)
##   names_list = names(model.df)
##   deleteat!(names_list, findall((in)([:EntryOrder]), names_list))
##   select!(model.df, names_list)
##   #
##   # Find which predictors have complete data.
##   # Note that the completeness_mask should be applied after the sample_mask.
##   #
##   completeness_mask = completecases(model.df)
  #
  # First consider the base model with no SNPs included.
  # If using one of the standard regression models without interactions,
  # then use the fast regression code in OpenMendel's general utilities file.
  # Otherwise, for general distributions, use the code in the GLM package.
  #
  if fast_method
    #
    # Create a model matrix using the regression formula and the model_frame.
    # The model_frame is just a subset and rearrangement of the person_frame
    # (a subset due to minimum genotype success rates, and a rearrangement
    # due to the need to match the original SNP Array data order).
    # The model matrix has columns for the intercept, any interaction terms,
    # and any variable functions (such as "log(BMI)").
    # The functions and macros used here are in the StatsModels package.
    # Here y is the vector of responses and X is the matrix of predictors.
    #
    fm = @eval(@formula($lhs ~ $rhs))
##    fm = @eval(@formula($side[1] ~ $side[2]))
    schem = schema(fm, model_frame)
    fs = apply_schema(fm, schem)
    response, X = modelcols(fs, model_frame)
    predictors = size(X, 2)
    y = zeros(cases)
    y[1:end] = model_frame[!, lhs]
    predictor_temp = coefnames(fs)[2]
    if isa(predictor_temp, Array)
      predictor_name = predictor_temp
      if predictors != size(predictor_name, 1)
        throw(ArgumentError("Inconsistency in number of predictors.\n \n"))
      end
    elseif isa(predictor_temp, String) && predictors == 1
      predictor_name = [predictor_temp]
    else
      throw(ArgumentError("Inconsistency in number and type of predictors.\n \n"))
    end
##    #
##    # Create the model matrix, initially still with any missing predictor values.
##    # Copy only the complete rows into design matrix X and response vector y,
##    # thus they have size = people - too.few.genotypes - incomplete.predictors.
##    #
##    modelmatrx = ModelMatrix(model)
##    cases = sum(completeness_mask)
##    predictors = size(modelmatrx.m, 2)
##    X = zeros(cases, predictors)
##    base_estimate = zeros(predictors)
##    for j = 1:predictors
##      X[:, j] = modelmatrx.m[completeness_mask, j]
##    end
##    y = zeros(cases)
##    y[1:end] = model.df[completeness_mask, lhs]
    #
    # Estimate parameters under the base model.
    # Then, for linear regression, record the vector of residuals.
    #
    base_estimate = zeros(predictors)
    (base_estimate, base_loglikelihood) = fast_regress(X, y, regression_type)
    if regression_type == "linear"
      base_residual = y - (X * base_estimate)
    end
    #
    # Output the results of the base model.
    #
    println(io, "\nOutput for MendelGWAS\n")
    println(io, "Using regression type: ", regression_type)
    println(io, "with the canonical link function.")
    println(io, "\n \nResults for Base Model:\n  ", fm, "\n")
    println(io, "Base components' effect estimates: ")
    for j = 1:predictors
      println(io, "   ", string(predictor_name[j]), " : ",
        round(base_estimate[j], sigdigits = 6))
    end
    println(io, "Base model loglikelihood: ",
      round(base_loglikelihood, sigdigits = 8))
    println(io, " ")
  else
    #
    # For other distributions or when including interactions,
    # let the GLM package estimate parameters. Output results.
    # Note that model_frame only contains complete cases.
    #
    fm = @eval(@formula($lhs ~ $rhs))
    base_model = glm(fm, model_frame, distribution_family, link)
    println(io, "\nOutput for MendelGWAS\n")
    if regression_type == ""
      println(io, "Using distribution type: ", distribution_family)
      println(io, "and link function: ", link)
    else
      println(io, "Using regression type: ", regression_type)
      println(io, "with the canonical link function.")
    end
    println(io, "\n \nResults for Base Model:\n  ", base_model.mf.f, "\n")
    print(io, coeftable(base_model))
    println(io, "\n \n")
  end
  flush(io)
  #
  # Now consider the alternative model with a SNP included.
  # Add a column to the design matrix to hold the SNP dosages.
  # Change the regression formula to include the SNP.
  #
  if fast_method; X = [X zeros(cases)]; end
  rhs = Meta.parse(side[2] * " + SNP")
  fm = @eval(@formula($lhs ~ $rhs))
  #
  # Analyze the SNP predictors one by one.
  #
  dosage = zeros(sum(sample_mask))
  pvalue = ones(snps)
  beta_snp = fill(NaN, snps)
  if fast_method
    alt_estimate = zeros(predictors+1)
    if regression_type == "linear"
      alt_residual = zeros(cases)
    end
  end
  skipped_snps = 0
  signif_snps = 0
  chr_max_bp = zeros(Integer, 50) # Assumes number of chromosomes <= 50.
  current_chr_number = 1
  current_chr = snpdata.chromosome[1]
  for snp = 1:snps
    #
    # Find maximum basepair listed for each chromosome.
    # Note that this section assumes SNPs on the same chromosome
    # are listed in a contigous block.
    #
    if snpdata.chromosome[snp] != current_chr
      current_chr = snpdata.chromosome[snp]
      current_chr_number = current_chr_number + 1
    end
    if !ismissing(snpdata.basepairs[snp])
      chr_max_bp[current_chr_number] =
        max(chr_max_bp[current_chr_number], snpdata.basepairs[snp])
    end
    #
    # Skip analysis of SNPs with MAF or genotype success rate
    # below the specified thresholds.
    #
    if !snp_mask[snp]
      skipped_snps = skipped_snps + 1
      continue
    end
    #
    # Copy the filtered SNP genotypes into a dosage vector,
    # allowing missing genotypes to be simplistically imputed based on MAF.
    #
    copyto!(dosage, view(snpdata.snpmatrix, sample_mask, snp); impute = true)
    #
    # For the three basic regression types, analyze the alternative model
    # using internal score test code. If the score test p-value
    # is below the specified threshold, carry out a likelihood ratio test.
    #
    if fast_method
##      copyto!(X[:, end], dosage[completeness_mask])
##      copyto!(X[:, end], view(dosage, completeness_mask))
##      X[:, end] = deepcopy(dosage[completeness_mask])
##      X[:, end] .= view(dosage, completeness_mask)
      X[:, end] = dosage[completeness_mask]
##if snp == 1; println(" X[:, end] = ", typeof(X), " :: ", X[:, end]); end
      alt_estimate = [base_estimate; 0.0]
      score_test = fast_score_test(X, y, alt_estimate, regression_type)
      pvalue[snp] = ccdf(Chisq(1), score_test)
      if pvalue[snp] < lrt_threshold
        (alt_estimate, alt_loglikelihood) = fast_regress(X, y, regression_type)
        lrt = 2.0 * (alt_loglikelihood - base_loglikelihood)
        pvalue[snp] = ccdf(Chisq(1), lrt)
	beta_snp[snp] = alt_estimate[end]
      end
      #
      # Record the vector of residuals under this alternative model:
      # alt_residual = y - (X * alt_estimate)
      #
      if regression_type == "linear"
        mul!(alt_residual, X, alt_estimate) # alt_residual = X * alt_estimate
        alt_residual .= y .- alt_residual
      end
    #
    # For other distributions, analyze the alternative model using GLM package.
    # Add the (abridged) SNP data to the design matrix and then fit the model.
    #
    else
      model_frame[!, :SNP] = dosage[completeness_mask]
      snp_model = fit(GeneralizedLinearModel, fm, model_frame,
        distribution_family, link)
      #
      # Find which row in the output table holds the SNP predictor results.
      # Store the p-value and beta value for the SNP.
      #
      snp_predictor = findfirst(isequal("SNP"), coeftable(snp_model).rownms)
      if snp_predictor == 0 || snp_predictor == nothing
       throw(ArgumentError("SNP predictor not found in alternative model.\n\n"))
      end
      pvalue[snp] = coeftable(snp_model).cols[4][snp_predictor]
      beta_snp[snp] = coeftable(snp_model).cols[1][snp_predictor]
    end
    #
    # Output regression results for potentially significant SNPs.
    #
    if pvalue[snp] < signif_threshold
      signif_snps = signif_snps + 1
      if signif_snps == 1
        println(io, "\n\nSummary of Significant Results for Alternative Model:")
        if fast_method
          println(io, "  ", fm)
        else
          println(io, "  ", snp_model.mf.f)
        end
        println(io, "using p-value significance threshold: ", signif_threshold)
      end
      println(io, "\nResults for SNP ", snpdata.snpid[snp])
      println(io, " on chromosome ", snpdata.chromosome[snp],
        " at basepair ", snpdata.basepairs[snp])
      println(io, "SNP p-value: ", round(pvalue[snp], sigdigits = 4))
      println(io, "Minor allele frequency: ",
        round(snpdata.maf[snp], sigdigits = 4))
      if uppercase(snpdata.chromosome[snp]) == "X"
        hw = xlinked_hardy_weinberg_test(dosage, person.male)
      else
        hw = hardy_weinberg_test(dosage)
      end
      println(io, "Hardy-Weinberg p-value: ", round(hw, sigdigits = 4))
      if fast_method && pvalue[snp] < lrt_threshold
        println(io, "SNP effect estimate: ",
          round(alt_estimate[end], sigdigits = 4))
        println(io, "SNP model loglikelihood: ",
          round(alt_loglikelihood, sigdigits = 8))
      elseif fast_method
        println(io, "SNP effect estimate: ",
          round(alt_estimate[end], sigdigits = 4))
      else
        println(io, "\n", coeftable(snp_model), "\n")
      end
      #
      # For linear models, output the proportion of the base model's variance
      # that is explained by including this SNP in the model.
      #
      if fast_method && regression_type == "linear"
        variance_explained = 1 - (norm(alt_residual)^2 / norm(base_residual)^2)
        println(io, "Proportion of base-model variance explained: ",
          round(variance_explained, sigdigits = 4))
      end
    end
##if snp == 1; return; end
  end
  #
  # Report if there were no significant SNPs.
  #
  if signif_snps == 0
    println(io, "\n\nSummary of Significant Results for Alternative Model:")
    if fast_method
      println(io, "  ", fm)
    else
      println(io, "  ", snp_model.mf.f)
    end
    println(io, "using p-value significance threshold: ", signif_threshold)
    println(io, "\nNo SNPs passed the significance threshold.\n")
  end
  #
  # Output false discovery rates.
  #
  fdr = [0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90]
  (number_passing, threshold) =
    simes_fdr(pvalue[snp_mask], fdr, snps - skipped_snps)
  println(io, " \n \n ")
  println(io, "        P-value   Number of Passing")
  println(io, "FDR    Threshold     Predictors \n")
  for i = 1:length(fdr)
   @printf(io, "%4.2f   %8.5f   %9i\n", fdr[i], threshold[i], number_passing[i])
  end
  println(io, " ")
  flush(io)
  #
  # Create a new dataframe that will hold the data to output.
  #
  output_frame = DataFrame(SNP = snpdata.snpid,
    Chromosome = snpdata.chromosome,
    BasePair = snpdata.basepairs,
    Pvalue = pvalue,
    NegLog10Pvalue = -log10.(pvalue),
    Beta = beta_snp)
  #
  # If requested, output a full listing of the p-values in csv format.
  #
  table_file = string(keyword["output_table"])
  if table_file != ""
    CSV.write(table_file, output_frame;
      writeheader = true, delim = keyword["output_field_separator"],
      missingstring = keyword["output_missing_value"])
    println(" \nTable of p-values for all SNPs is in file: ", table_file)
  end
  #
  # If requested, output a Manhattan plot or QQ plot.
  #
  manhat_plot_file = string(keyword["manhattan_plot_file"])
  qq_plot_file = string(keyword["qq_plot_file"])
  if manhat_plot_file != "" || qq_plot_file != ""
    println(" \nPlots being drawn using MendelPlots package.")
  end
  if manhat_plot_file != ""
    manhattan(output_frame, outfile = manhat_plot_file,
##      linecolor = colorant"black",
      chrvar = "Chromosome", posvar = "BasePair", pvalvar = "Pvalue")
    println(" \nManhattan plot of GWAS results is in file: ", manhat_plot_file)
  end
  if qq_plot_file != ""
    qq(output_frame, outfile = qq_plot_file, pvalvar = "Pvalue")
    println(" \nQQ plot of GWAS results is in file: ", qq_plot_file)
  end
  return execution_error = false
end # function gwas_option

##!!"""
##!!Add principal components into a pedigree frame and regression formula.
##!!The PCs are calculated via pca() from the SnpArrays module.
##!!"""
##!!function add_pcs!(person_frame::DataFrame, snpdata::SnpDataStruct,
##!!  keyword::Dict{AbstractString, Any})
##!!
##!!  pcs = keyword["pcs"]
##!!  regress_form = keyword["regression_formula"]
##!!  #
##!!  # Perform PCA on the SNP data.
##!!  #
##!!  pcscore, pcloading, pcvariance = pca(snpdata.snpmatrix, pcs)
##!!  #
##!!  # Include the new covariates in the pedigree dataframe.
##!!  # NB: this process is *slow* for a large number of PCs!
##!!  #
##!!  for i = 1:pcs
##!!    person_frame[!, Symbol("PC$i")] = zscore(pcscore[:,i])
##!!    regress_form = regress_form * " + PC$i"
##!!  end
##!!  #
##!!  # Save the expanded regression formula to the keyword dictionary.
##!!  #
##!!  keyword["regression_formula"] = regress_form
##!!  return nothing
##!!end # function add_pcs!
#
# Method to obtain path to this package's data files
# so they can be used in the documentation and testing routines.
# For example, datadir("Control file.txt") will return
# "/path/to/package/data/Control file.txt"
#
datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)

end # module MendelGWAS


text:
"""
    StaticPokemon(types, stats, fastMove, chargedMoves)

Struct for holding the values associated with the mons that do not change
throughout the battle: types, stats, and moves. Note that like moves, this
struct is agnostic to the actual identity/dex/species of the mon.
"""
struct StaticPokemon
    primary_type::Int8
    secondary_type::Int8
    stats::Stats
    fastMove::FastMove
    charged_move_1::ChargedMove
    charged_move_2::ChargedMove
end

"""
    StaticPokemon(i; league = "great", cup = "all",
        custom_moveset = ["none"], custom_stats = ())

Construct a StaticPokemon from the index of a mon within its rankings
(optionally specified). Other optional inputs are a custom moveset or IVs.
"""
function StaticPokemon(i::Int64; league::String = "great", cup = "all",
  custom_moveset = ["none"], custom_stats = ())
    rankings = get_rankings(cup, league = league)
    gmid = get_gamemaster_mon_id(rankings[i]["speciesId"])
    gm = gamemaster["pokemon"][gmid]
    types = typings[convert(Array{String}, gm["types"])[1]],
        typings[convert(Array{String}, gm["types"])[2]]
    cp_limit = get_cp_limit(league)
    if custom_stats != ()
        level, atk, def, hp = parse.(Int8, custom_stats)
        if level == 0
            function get_cp(lvl)
                attack = (atk + gm["baseStats"]["atk"]) * cpm[lvl]
                defense = (def + gm["baseStats"]["def"]) * cpm[lvl]
                hitpoints = floor(Int16,
                    (hp + gm["baseStats"]["hp"]) * cpm[lvl])
                cp = floor(max(10,
                    (attack * sqrt(defense) * sqrt(hitpoints)) / 10.0))
                return cp
            end
            level = (1:0.5:40)[findfirst(x -> get_cp(x) > cp_limit,
                1:0.5:40) - 1]
        end
    elseif league == "master"
        level, atk, def, hp = 40, 15, 15, 15
    else
        level = gm["defaultIVs"]["cp$(cp_limit)"][1]
        atk = gm["defaultIVs"]["cp$(cp_limit)"][2]
        def = gm["defaultIVs"]["cp$(cp_limit)"][3]
        hp = gm["defaultIVs"]["cp$(cp_limit)"][4]
    end
    attack = haskey(gm, "tags") && "shadow" in gm["tags"] ? floor(UInt16,
        (6/5 * atk + gm["baseStats"]["atk"]) * cpm[level] * 100) :
        floor(UInt16, (atk + gm["baseStats"]["atk"]) * cpm[level] * 100)
    defense = haskey(gm, "tags") && "shadow" in gm["tags"] ? floor(UInt16,
        (5/6 * def + gm["baseStats"]["def"]) * cpm[level] * 100) :
        floor(UInt16, (def + gm["baseStats"]["def"]) * cpm[level] * 100)
    hitpoints = floor(Int16, (hp + gm["baseStats"]["hp"]) * cpm[level])
    stats = Stats(attack, defense, hitpoints)
    if haskey(rankings[i], "moveStr")
        moves = parse.(Ref(Int64), split(rankings[i]["moveStr"], "-"))
        fastMovesAvailable = gm["fastMoves"]
        sort!(fastMovesAvailable)
        fastMoveGm = gamemaster["moves"][get_gamemaster_move_id(
            fastMovesAvailable[moves[1]+1],)]
        fastMove = Move(fastMoveGm, types)
        chargedMovesAvailable = gm["chargedMoves"]
        if haskey(gm, "tags") &&
           "shadoweligible" in gm["tags"] && gm["level25CP"] < cp_limit
            push!(chargedMovesAvailable, "RETURN")
        elseif haskey(gm, "tags") && "shadow" in gm["tags"]
            push!(chargedMovesAvailable, "FRUSTRATION")
        end
        sort!(chargedMovesAvailable)
        chargedMove1 = Move(gamemaster["moves"][get_gamemaster_move_id(
            chargedMovesAvailable[moves[2]])], types)
        chargedMove2 = Move(gamemaster["moves"][get_gamemaster_move_id(
            chargedMovesAvailable[moves[3]])], types)
    else
        moveset = custom_moveset == ["none"] ?
            rankings[i]["moveset"] : custom_moveset
        fastMove = FastMove(moveset[1]::String, types)
        chargedMove1 = ChargedMove(moveset[2]::String, types)
        chargedMove2 = ChargedMove(moveset[3]::String, types)
    end
    return StaticPokemon(
        types[1],
        types[2],
        stats,
        fastMove,
        chargedMove1.energy > chargedMove2.energy ? chargedMove2 : chargedMove1,
        chargedMove1.energy > chargedMove2.energy ? chargedMove1 : chargedMove2
    )
end

"""
    StaticPokemon(mon; league = "great", cup = "all")

Construct a StaticPokemon from the name of the pokemon, and the meta it is
within. Movesets and IVs can also be specified by comma-separating the string
being passed in.
"""
function StaticPokemon(mon::String; league = "great", cup = "all")
    if occursin(",", mon)
        mon_arr = split(mon, ",")
        if length(mon_arr) == 4
            return StaticPokemon(get_rankings_mon_id(convert(String,
                mon_arr[1]), league = league, cup = cup), league = league,
                cup = cup, custom_moveset = convert.(String, mon_arr[2:4]))
        elseif length(mon_arr) == 7
            return StaticPokemon(get_rankings_mon_id(convert(String,
                mon_arr[1])), league = league, cup = cup,
                custom_moveset = convert.(String, mon_arr[2:4]),
                custom_stats = ("0", mon_arr[5], mon_arr[6], mon_arr[7]))
        elseif length(mon_arr) == 8
            return StaticPokemon(get_rankings_mon_id(convert(String,
                mon_arr[1])), league = league, cup = cup,
                custom_moveset = convert.(String, mon_arr[2:4]),
                custom_stats = (mon_arr[5], mon_arr[6], mon_arr[7], mon_arr[8]))
        end
    else
        return StaticPokemon(get_rankings_mon_id(
            mon, league = league, cup = cup),
            league = league, cup = cup)
    end
end

"""
    DynamicPokemon(hp, energy)

Struct for holding the values associated with the mons that change throughout
the battle: current hp and energy. Note that like moves, this struct is
agnostic to the actual identity/dex/species of the mon.
"""
struct DynamicPokemon
    data::UInt16    # hp (initially hitpoints stat) and energy (initially 0)
end

"""
    DynamicPokemon(mon)

Construct a starting DynamicPokemon from a StaticPokemon. This is just setting
the starting hp of the mon to the stat value, and the energy to zero.
"""

function get_energy(p::DynamicPokemon)
    return p.data >> 0x0009
end

function get_hp(p::DynamicPokemon)
    return p.data % 0x0200
end

function add_energy(p::DynamicPokemon, e::Int8)
    return DynamicPokemon(min(0x0064, get_energy(p) + e) << 9 + get_hp(p))
end

function subtract_energy(p::DynamicPokemon, e::Int8)
    curr_e = get_energy(p)
    return DynamicPokemon((curr_e - min(curr_e, e)) << 9 + get_hp(p))
end

function damage(p::DynamicPokemon, d::UInt16)
    return DynamicPokemon(p.data - min(get_hp(p), d))
end

function DynamicPokemon(mon::StaticPokemon)
    DynamicPokemon(mon.stats.hitpoints)
end


text:
# # Simple diffusion example
#
# This script provides our simplest example of Oceananigans.jl functionality:
# the diffusion of a one-dimensional Gaussian. This example demonstrates
#
#   * how to load `Oceananigans.jl`;
#   * how to instantiate an `Oceananigans.jl` model;
#   * how to create simple `Oceananigans.jl` output;
#   * how to set an initial condition with a function;
#   * how to time-step a model forward, and finally
#   * how to look at results.
#
# ## Using `Oceananigans.jl`
#
# To use `Oceananigans.jl` after it has been installed, we bring
# `Oceananigans.jl` functions and names into our 'namespace' by writing

using Oceananigans

# In addition, we import the submodule `Grids`.

using Oceananigans.Grids

# ## Instantiating and configuring a model
#
# To begin using Oceananigans, we instantiate an incompressible model by calling
# the `IncompressibleModel` constructor:

model = IncompressibleModel(
       grid = RegularCartesianGrid(size = (1, 1, 128), x = (0, 1), y = (0, 1), z = (-0.5, 0.5)),
    closure = IsotropicDiffusivity(κ = 1.0)
)
nothing # hide

# The keyword arguments `grid` and `closure` indicate that
# our model grid is Cartesian with uniform grid spacing, that our diffusive
# stress and tracer fluxes are determined by diffusion with a constant
# diffusivity `κ` (note that we do not use viscosity in this example).

# Note that by default, a `Model` has no-flux boundary condition on all
# variables. Next, we set an initial condition on our "passive tracer",
# temperature. Our objective is to observe the diffusion of a Gaussian.

## Build a Gaussian initial condition function with width `δ`:
δ = 0.1
Tᵢ(x, y, z) = exp(-z^2 / (2δ^2))

## Set `model.tracers.T` to the function `Tᵢ`:
set!(model, T=Tᵢ)

# ## Running your first `Model`
#
# Finally, we time-step the model forward using the function
# `time_step!`, with a time-step size that ensures numerical stability.

## Time-scale for diffusion across a grid cell
cell_diffusion_time_scale = model.grid.Δz^2 / model.closure.κ.T

## We create a `Simulation` which will handle time stepping the model. It will
## execute `Nt` time steps with step size `Δt` using a second-order Adams-Bashforth method.
simulation = Simulation(model, Δt = 0.1 * cell_diffusion_time_scale, stop_iteration = 1000)

run!(simulation)

# ## Visualizing the results
#
# We use `Plots.jl` to look at the results. Fields are
# stored as 3D arrays in Oceananigans so we plot `interior(T)[1, 1, :]`
# which will return a 1D array.

using Plots, Printf

## A convenient function for generating a label with the current model time
tracer_label(time) = @sprintf("t = %.3f", time)

## Plot initial condition
T = model.tracers.T

z = znodes(T)[:]

p = plot(Tᵢ.(0, 0, z), z, linewidth=2, label="t = 0",
         xlabel="Temperature", ylabel="z")

## Plot current solution
plot!(p, interior(T)[1, 1, :], z, linewidth=2, label=tracer_label(model.clock.time))

# Interesting! Next, we add an output writer that saves the temperature field
# in JLD2 files, and run the simulation for longer so that we can animate the results.

using Oceananigans.OutputWriters: JLD2OutputWriter

simulation.output_writers[:temperature] =
    JLD2OutputWriter(model, model.tracers, prefix = "one_dimensional_diffusion",
                     iteration_interval = 100, force = true)

## Run simulation for 10,000 more iterations
simulation.stop_iteration += 10000

run!(simulation)
nothing

# Finally, we animate the results by opening the JLD2 file, extract the
# iterations we ended up saving at, and plot the evolution of the
# temperature profile in a loop over the iterations.

using JLD2

file = jldopen(simulation.output_writers[:temperature].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)

    T = file["timeseries/T/$iter"][1, 1, :]
    t = file["timeseries/t/$iter"]

    plot(T, z, linewidth=2, title=tracer_label(t),
         label="", xlabel="Temperature", ylabel="z", xlims=(0, 1))
end

mp4(anim, "one_dimensional_diffusion.mp4", fps = 15) # hide


text:
# included in runtests.jl
using Test
import Sockets: TCPSocket
import Random: randstring
using WebSockets
import WebSockets: maskswitch!,
    write_fragment,
    read_frame,
    is_control_frame,
    handle_control_frame,
    locked_write,
    codeDesc
include("logformat.jl")

"""
The dummy websocket don't use TCP. Close won't work, but we can manipulate the contents
using otherwise the same functions as TCP sockets.
"""
dummyws(server::Bool)  = WebSocket(Base.BufferStream(), server)
io = IOBuffer()

# maskswitch
empty1 = UInt8[]
empty2 = UInt8[]
@test length(maskswitch!(empty1)) == 4
@test empty1 == empty2
# Test most basic frame, length <126

for len = [8, 125], fin=[true, false], clientwriting = [false, true]

    op = (rand(UInt8) & 0b1111)
    test_str = randstring(len)
    # maskswitch two times with same key == unmasked
    maskunmask = copy(codeunits(test_str))
    mskkey = maskswitch!(maskunmask)
    maskswitch!(maskunmask, mskkey)
    @test maskunmask == codeunits(test_str)

    # websocket fragment as Base.CodeUnits{UInt8,String}
    # for client writing, the data is masked and the mask is contained in the frame.
    # for server writing, the data is not masked, and the header is four bytes shorter.
    write_fragment(io, fin, op, clientwriting, copy(codeunits(test_str)))
    # test that the original input string was not masked.
    @test maskunmask == codeunits(test_str)
    frame = take!(io)
    # Check the frame header
    # Last frame bit
    @test bitstring(frame[1]) == (fin ? "1" : "0") * "000" * bitstring(op)[end-3:end]
    # payload length bit
    @test frame[2] & 0b0111_1111 == len
    # ismasked bit
    hasmsk = (frame[2] & 0b1000_0000) >>> 7 != 0
    @test hasmsk  == clientwriting
    # payload data
    if hasmsk
        framedata = copy(frame[7:end])
        maskswitch!(framedata, frame[3:6])
    else
        framedata = frame[3:end]
    end

    @test framedata == codeunits(test_str)

    # Test for WebSocketError when reading
    #  masked frame-> websocket|server
    #  unmasked frame -> websocket|client

    # Let's pretend TCP has moved our frame into the peer websocket
    receivingws = dummyws(!clientwriting)
    write(receivingws.socket, frame)
    @test_throws WebSockets.WebSocketError read_frame(receivingws)
    close(receivingws.socket)

    # Let's pretend receivingws didn't error like it should, but
    # echoed our message back with identical masking.
    dws = dummyws(clientwriting)
    @test dws.server == clientwriting
    write(dws.socket, frame)
    # read the frame back, now represented as a WebSocketFragment

    frag_back = read_frame(dws)
    close(dws.socket)
    @test frag_back.is_last == fin
    @test frag_back.rsv1 == false
    @test frag_back.rsv2 == false
    @test frag_back.rsv3 == false
    @test frag_back.opcode == op
    @test frag_back.is_masked == clientwriting
    @test frag_back.payload_len == len
    maskkey = UInt8[]
    if clientwriting
        maskkey = frame[3:6]
    end
    @test frag_back.maskkey == maskkey
    # the WebSocketFragment stores the data after unmasking
    @test codeunits(test_str) == frag_back.data
end

# Test length 126 or more

for len = 126:129, fin=[true, false], clientwriting = [false, true]
    op = 0b1111
    test_str = randstring(len)
    write_fragment(io, fin, op, clientwriting, copy(codeunits(test_str)))
    frame = take!(io)

    @test bitstring(frame[1]) == (fin ? "1" : "0") * "000" * bitstring(op)[end-3:end]
    @test frame[2] & 0b0111_1111 == 126
    @test bitstring(frame[4])*bitstring(frame[3]) == bitstring(hton(UInt16(len)))

    dws = dummyws(clientwriting)
    write(dws.socket, frame)
    frag_back = read_frame(dws)
    close(dws.socket)

    @test frag_back.is_last == fin
    @test frag_back.rsv1 == false
    @test frag_back.rsv2 == false
    @test frag_back.rsv3 == false
    @test frag_back.opcode == op
    @test frag_back.is_masked == clientwriting
    @test frag_back.payload_len == len
    @test test_str == String(frag_back.data)
end

# Test length typemax(UInt16) + 1

for clientwriting = [false, true]
    len = typemax(UInt16) +1
    op = 0b1111
    fin = true

    test_str = randstring(len)
    write_fragment(io, fin, op, clientwriting, copy(codeunits(test_str)))
    frame = take!(io)

    dws = dummyws(clientwriting)
    write(dws.socket, frame)
    frag_back = read_frame(dws)
    close(dws.socket)

    @test frag_back.is_last == fin
    @test frag_back.rsv1 == false
    @test frag_back.rsv2 == false
    @test frag_back.rsv3 == false
    @test frag_back.opcode == op
    @test frag_back.is_masked == clientwriting
    @test frag_back.payload_len == len
    @test test_str == String(frag_back.data)
end

# Test unknown opcodes

for op in 0xB:0xF
    clientwriting = false
    len = 10
    fin = true

    test_str = randstring(len)
    write_fragment(io, fin, op, clientwriting, copy(codeunits(test_str)))
    frame = take!(io)

    dws = dummyws(clientwriting)
    write(dws.socket, frame)
    frag_back = read_frame(dws)

    @test is_control_frame(frag_back)
    thiserror = ArgumentError("")
    try
        handle_control_frame(dws, frag_back)
    catch err
        thiserror = err
    end
    @test typeof(thiserror) <: ErrorException
    @test thiserror.msg == " while handle_control_frame(ws|client, wsf): Unknown opcode $op"

    close(dws.socket)
end

# Test multi-frame message

for clientwriting = [false, true]

    op = WebSockets.OPCODE_TEXT
    full_str = "123456"
    first_str = "123"
    second_str = "456"
    fin = false
    dws = dummyws(clientwriting)
    write_fragment(dws.socket, fin, op, clientwriting, copy(codeunits(first_str)))
    fin = true
    write_fragment(dws.socket, fin, op, clientwriting, copy(codeunits(second_str)))

    @test read(dws) == codeunits(full_str)
end

# Test read(ws) bad mask error handling

@info "Provoking close handshake from protocol error without a peer. Waits 10s, 'a reasonable time'."
for clientwriting in [false, true]
    op = WebSockets.OPCODE_TEXT
    test_str = "123456"
    fin = true
    write_fragment(io, fin, op, clientwriting, copy(codeunits(test_str)))
    frame = take!(io)
    # let's put this frame on the same kind of socket, and then read it as if it came from the peer
    # This will provoke a close handshake, but since there is no peer it times out.
    dws = dummyws(!clientwriting)
    write(dws.socket, frame)
    thiserror = ArgumentError("")
    try
        read(dws)
    catch err
        thiserror = err
    end
    @test typeof(thiserror) <: WebSocketClosedError
    expmsg = " while read(ws|$(dws.server ? "server" : "client")) WebSocket|$(dws.server ? "server" : "client") cannot handle incoming messages with$(dws.server ? "out" : "") mask. Ref. rcf6455 5.3 - Performed closing handshake."
    @test thiserror.message ==  expmsg
    @test !isopen(dws)
    close(dws.socket)

    # simple close frame
    clientwriting = false
end

# Close frame with no reason.
for clientwriting = [false, true]
    fin = true
    op = WebSockets.OPCODE_CLOSE
    write_fragment(io, fin, op, clientwriting, UInt8[])
    frame = take!(io)
    len = 0
    # Check the frame header
    # Last frame bit
    @test bitstring(frame[1]) == (fin ? "1" : "0") * "000" * bitstring(op)[end-3:end]
    # payload length bit
    @test frame[2] & 0b0111_1111 == len
    # ismasked bit
    hasmsk = (frame[2] & 0b1000_0000) >>> 7 != 0
    @test hasmsk  == clientwriting
    # the peer of the writer is
    dws = dummyws(clientwriting)
    write(dws.socket, frame)
    frag_back = read_frame(dws)
    @test frag_back.is_last == fin
    @test frag_back.opcode == op
    @test frag_back.is_masked == clientwriting
    @test frag_back.payload_len == len
    @test is_control_frame(frag_back)
end

# Close with no reason
for clientwriting in [false, true]
    op = WebSockets.OPCODE_CLOSE
    fin = true
    thisws = dummyws(!clientwriting)
    locked_write(thisws.socket, true, op, !thisws.server, UInt8[])
    close(thisws.socket)
    frame = read(thisws.socket)
    peerws = dummyws(clientwriting)
    write(peerws.socket, frame)
    close(peerws.socket)
    wsf = read_frame(peerws)
    @test is_control_frame(wsf)
    @test wsf.opcode == WebSockets.OPCODE_CLOSE
    @test wsf.payload_len == 0
end

# Close with status number
for clientwriting in [false, true], statusnumber in keys(codeDesc)
    op = WebSockets.OPCODE_CLOSE
    freereason = ""
    fin = true
    thisws = dummyws(!clientwriting)
    statuscode = reinterpret(UInt8, [hton(UInt16(statusnumber))])
    locked_write(thisws.socket, true, op, !thisws.server, statuscode)
    close(thisws.socket)
    frame = read(thisws.socket)

    # Check the frame header
    # Last frame bit
    @test bitstring(frame[1]) == (fin ? "1" : "0") * "000" * bitstring(op)[end-3:end]
    # payload length bit
    @test frame[2] & 0b0111_1111 == 2
    # ismasked bit
    hasmsk = (frame[2] & 0b1000_0000) >>> 7 != 0
    @test hasmsk  == clientwriting
    # the peer of the writer is
    peerws = dummyws(clientwriting)
    write(peerws.socket, frame)
    close(peerws.socket)
    wsf = read_frame(peerws)
    @test is_control_frame(wsf)
    @test wsf.opcode == WebSockets.OPCODE_CLOSE
    @test wsf.payload_len == 2
    scode = Int(reinterpret(UInt16, reverse(wsf.data))[1])
     @test scode == statusnumber
    reason = string(scode) * ":" * get(codeDesc, scode, "")
end

# Close with status number and freereason
for clientwriting in [false, true], statusnumber in keys(codeDesc)
    freereason = "q.e.d"
    op = WebSockets.OPCODE_CLOSE
    fin = true
    thisws = dummyws(!clientwriting)
    statuscode = vcat(reinterpret(UInt8, [hton(UInt16(statusnumber))]),
                codeunits(freereason))
    locked_write(thisws.socket, true, op, !thisws.server, copy(statuscode))
    close(thisws.socket)
    frame = read(thisws.socket)

    # Check the frame header
    # Last frame bit
    @test bitstring(frame[1]) == (fin ? "1" : "0") * "000" * bitstring(op)[end-3:end]
    # payload length bit
    @test frame[2] & 0b0111_1111 == length(statuscode)
    # ismasked bit
    hasmsk = (frame[2] & 0b1000_0000) >>> 7 != 0
    @test hasmsk  == clientwriting
    # the peer of the writer is
    peerws = dummyws(clientwriting)
    write(peerws.socket, frame)
    close(peerws.socket)
    wsf = read_frame(peerws)
    @test is_control_frame(wsf)
    @test wsf.opcode == WebSockets.OPCODE_CLOSE
    @test wsf.payload_len == length(statuscode)
    scode = Int(reinterpret(UInt16, reverse(wsf.data[1:2]))[1])
    reason = string(scode) * ":" *
            get(codeDesc, scode, "") *
            " " * String(wsf.data[3:end])
    @test reason == string(statusnumber) * ":" *
            get(codeDesc, statusnumber, "") *
            " " * freereason
end

close(io)


text:
# Autogenerated wrapper script for Qt5QuickControls2_jll for x86_64-apple-darwin
export libqt5quickcontrols2, libqt5quicktemplates2, quickcontrols2

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5QuickControls2")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "@rpath/QtQuickControls2")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "@rpath/QtQuickTemplates2")
JLLWrappers.@declare_file_product(quickcontrols2)
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "lib/QtQuickControls2.framework/Versions/5/QtQuickControls2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "lib/QtQuickTemplates2.framework/Versions/5/QtQuickTemplates2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        quickcontrols2,
        "qml/QtQuick/Controls.2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()