text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ResourcesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks whether a resource exists.
Param: resourceGroupName::String (required)
Param: resourceProviderNamespace::String (required)
Param: parentResourcePath::String (required)
Param: resourceType::String (required)
Param: resourceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function resourcesCheckExistence(_api::ResourcesApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "resourcesCheckExistence", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "resourcesCheckExistence", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "HEAD", Nothing, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/{resourceProviderNamespace}/{parentResourcePath}/{resourceType}/{resourceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "parentResourcePath", parentResourcePath)  # type String
    Swagger.set_param(_ctx.path, "resourceType", resourceType)  # type String
    Swagger.set_param(_ctx.path, "resourceName", resourceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Checks by ID whether a resource exists.
Param: resourceId::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function resourcesCheckExistenceById(_api::ResourcesApi, resourceId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "HEAD", Nothing, "/{resourceId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceId", resourceId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Creates a resource.
Param: resourceGroupName::String (required)
Param: resourceProviderNamespace::String (required)
Param: parentResourcePath::String (required)
Param: resourceType::String (required)
Param: resourceName::String (required)
Param: api_version::String (required)
Param: parameters::GenericResource (required)
Param: subscriptionId::String (required)
Return: GenericResource
"""
function resourcesCreateOrUpdate(_api::ResourcesApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, parameters, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "resourcesCreateOrUpdate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "resourcesCreateOrUpdate", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "PUT", GenericResource, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/{resourceProviderNamespace}/{parentResourcePath}/{resourceType}/{resourceName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "parentResourcePath", parentResourcePath)  # type String
    Swagger.set_param(_ctx.path, "resourceType", resourceType)  # type String
    Swagger.set_param(_ctx.path, "resourceName", resourceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Create a resource by ID.
Param: resourceId::String (required)
Param: api_version::String (required)
Param: parameters::GenericResource (required)
Return: GenericResource
"""
function resourcesCreateOrUpdateById(_api::ResourcesApi, resourceId::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", GenericResource, "/{resourceId}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceId", resourceId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes a resource.
Param: resourceGroupName::String (required)
Param: resourceProviderNamespace::String (required)
Param: parentResourcePath::String (required)
Param: resourceType::String (required)
Param: resourceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function resourcesDelete(_api::ResourcesApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "resourcesDelete", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "resourcesDelete", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/{resourceProviderNamespace}/{parentResourcePath}/{resourceType}/{resourceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "parentResourcePath", parentResourcePath)  # type String
    Swagger.set_param(_ctx.path, "resourceType", resourceType)  # type String
    Swagger.set_param(_ctx.path, "resourceName", resourceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes a resource by ID.
Param: resourceId::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function resourcesDeleteById(_api::ResourcesApi, resourceId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/{resourceId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceId", resourceId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a resource.
Param: resourceGroupName::String (required)
Param: resourceProviderNamespace::String (required)
Param: parentResourcePath::String (required)
Param: resourceType::String (required)
Param: resourceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: GenericResource
"""
function resourcesGet(_api::ResourcesApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "resourcesGet", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "resourcesGet", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", GenericResource, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/{resourceProviderNamespace}/{parentResourcePath}/{resourceType}/{resourceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "parentResourcePath", parentResourcePath)  # type String
    Swagger.set_param(_ctx.path, "resourceType", resourceType)  # type String
    Swagger.set_param(_ctx.path, "resourceName", resourceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a resource by ID.
Param: resourceId::String (required)
Param: api_version::String (required)
Return: GenericResource
"""
function resourcesGetById(_api::ResourcesApi, resourceId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", GenericResource, "/{resourceId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceId", resourceId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Get all the resources in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: expand::String
Param: top::Int32
Return: ResourceListResult
"""
function resourcesList(_api::ResourcesApi, api_version::String, subscriptionId::String; filter=nothing, expand=nothing, top=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ResourceListResult, "/subscriptions/{subscriptionId}/resources", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""
Moves resources from one resource group to another resource group.
The resources to move must be in the same source resource group. The target resource group may be in a different subscription. When moving resources, both the source group and the target group are locked for the duration of the operation. Write and delete operations are blocked on the groups until the move completes.
Param: sourceResourceGroupName::String (required)
Param: parameters::ResourcesMoveInfo (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function resourcesMoveResources(_api::ResourcesApi, sourceResourceGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("sourceResourceGroupName", "resourcesMoveResources", :maxLength, sourceResourceGroupName, 90)
    Swagger.validate_param("sourceResourceGroupName", "resourcesMoveResources", :minLength, sourceResourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{sourceResourceGroupName}/moveResources", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "sourceResourceGroupName", sourceResourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""
Validates whether resources can be moved from one resource group to another resource group.
This operation checks whether the specified resources can be moved to the target. The resources to move must be in the same source resource group. The target resource group may be in a different subscription. If validation succeeds, it returns HTTP response code 204 (no content). If validation fails, it returns HTTP response code 409 (Conflict) with an error message. Retrieve the URL in the Location header value to check the result of the long-running operation.
Param: sourceResourceGroupName::String (required)
Param: parameters::ResourcesMoveInfo (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function resourcesValidateMoveResources(_api::ResourcesApi, sourceResourceGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("sourceResourceGroupName", "resourcesValidateMoveResources", :maxLength, sourceResourceGroupName, 90)
    Swagger.validate_param("sourceResourceGroupName", "resourcesValidateMoveResources", :minLength, sourceResourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{sourceResourceGroupName}/validateMoveResources", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "sourceResourceGroupName", sourceResourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export resourcesCheckExistence, resourcesCheckExistenceById, resourcesCreateOrUpdate, resourcesCreateOrUpdateById, resourcesDelete, resourcesDeleteById, resourcesGet, resourcesGetById, resourcesList, resourcesMoveResources, resourcesValidateMoveResources


text:
# License for this file: MIT (expat)
# Copyright 2016-2018, DLR Institute of System Dynamics and Control

"""
    module test_withoutMacros_withoutVariables1

Test function ModiaMath.simulate(..) with standard models (no events, no sparse Jacobian),
and without Variables and without Macros
"""
module test_withoutMacros_withoutVariables1

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
using ModiaMath.Test

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "PT1.jl"))
import .PT1

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "PendulumODE.jl"))
import .PendulumODE

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "PendulumDAE.jl"))
import .PendulumDAE

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "FreeBodyRotation.jl"))
import .FreeBodyRotation

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "StateSelection.jl"))
import .StateSelection

@testset "ModiaMath: withoutMacros_withoutVariables/*.jl without events" begin
    # Simulate PT1 block "T*der(x) + x = 0" in different variants
    m = PT1.Model(T=2.0, x0=1.5)
    interval = 0.1
    stopTime = 10.0
    t = collect(0.0:interval:stopTime)
    nt = length(t)

    # Analytic solution
    function fExact(m::PT1.Model, t::Vector{Float64})
        nt = length(t)
        xExact     = zeros(nt)
        der_xExact = zeros(nt)
        for i = 1:nt
            xExact[i,1]     = exp(-t[i] / m.T) * m.x0
            der_xExact[i,1] = -xExact[i] / m.T
        end
        return (xExact, der_xExact)
    end

    # Test whether simulation results in the exact solution
    @testset "Simulate PT1 with default tolerance" begin
        tolerance = 1e-4
        result    = ModiaMath.simulate!(m, stopTime=stopTime, interval=interval, tolerance=tolerance)
        (xExact, der_xExact) = fExact(m, t)

        @test isapprox(result["x[1]"], xExact    ; atol=10 * tolerance )
        @test isapprox(result["der(x[1])"], der_xExact; atol=100 * tolerance )
        @test m.simulationState.statistics.sparseSolver == false

    end

    @testset "Simulate PT1 with tolerance = 1e-8 and log=true" begin
        tolerance = 1.0e-8
        result = ModiaMath.simulate!(m, stopTime=stopTime, log=true, interval=interval, tolerance=tolerance)
        (xExact, der_xExact) = fExact(m, t)

        @test isapprox(result["x[1]"], xExact    ; atol=10 * tolerance )
        @test isapprox(result["der(x[1])"], der_xExact; atol=100 * tolerance )
        @test m.simulationState.statistics.sparseSolver == false
    end

    @testset "Simulate PendulumODE" begin
        tolerance = 1e-8
        m3 = PendulumODE.Model(L=0.5, m=1.0, d=0.1, phi0_deg=60.0)   # L,m,d,phi_deg
        result = ModiaMath.simulate!(m3, stopTime=stopTime, interval=interval, tolerance=tolerance)
        phi = result["phi"]
        w   = result["w"]

        @test isapprox(phi[end], 0.120289; atol=1e-3 )
        @test isapprox(w[end], 0.273965; atol=1e-2 )
    end

    @testset "Simulate PendulumDAE" begin
        stopTime  = 2.0
        tolerance = 1e-8
        m4 = PendulumDAE.Model()
        result = ModiaMath.simulate!(m4, stopTime=stopTime, interval=interval, tolerance=tolerance)
        x  = result["x"]
        y  = result["y"]
        vx = result["vx"]
        vy = result["vy"]
        dervx = result["der(vx)"]
        dervy = result["der(vy)"]

        @test isapprox(x[end],  0.542797; atol=1e-5 )
        @test isapprox(y[end], -0.839864; atol=1e-5 )
        @test isapprox(vx[end],  1.80223 ; atol=1e-5 )
        @test isapprox(vy[end],  1.16476 ; atol=1e-5 )
        @test isapprox(dervx[end], -6.97155 ; atol=1e-5 )
        @test isapprox(dervy[end],  0.977022; atol=1e-5 )
    end

    @testset "Simulate PendulumDAE with x_fixed=true" begin
        stopTime  = 2.0
        tolerance = 1e-8
        m4 = PendulumDAE.Model(x_fixed=true)
        result = ModiaMath.simulate!(m4, stopTime=stopTime, interval=interval, tolerance=tolerance)
        x  = result["x"]
        y  = result["y"]
        vx = result["vx"]
        vy = result["vy"]
        dervx = result["der(vx)"]
        dervy = result["der(vy)"]

        @test isapprox(x[end],  0.431501; atol=1e-5 )
        @test isapprox(y[end], -0.902112; atol=1e-5 )
        @test isapprox(vx[end],  1.2889  ; atol=1e-5 )
        @test isapprox(vy[end],  0.616512; atol=1e-5 )
        @test isapprox(dervx[end], -4.69952 ; atol=1e-5 )
        @test isapprox(dervy[end],  0.01497 ; atol=1e-5 )
    end

    @testset "Simulate FreeBodyRotation" begin
        stopTime  = 5.0
        tolerance = 1e-8
        m5 = FreeBodyRotation.Model()
        result = ModiaMath.simulate!(m5, stopTime=stopTime, interval=interval, tolerance=tolerance)

        #=
        @test isapprox(x[end]    ,  0.542797; atol=1e-3 )
        @test isapprox(y[end]    , -0.839864; atol=1e-3 )
        @test isapprox(vx[end]   ,  1.80223 ; atol=1e-3 )
        @test isapprox(vy[end]   ,  1.16476 ; atol=1e-3 )
        @test isapprox(dervx[end], -6.97155 ; atol=1e-3 )
        @test isapprox(dervy[end],  0.97702 ; atol=1e-3 )
        =#
    end

    @testset "Simulate State Selection" begin
        m6 = StateSelection.Model()
        result = ModiaMath.simulate!(m6, stopTime=1.0)
		f3 = result["f[3]"][end]

        @test isapprox(f3[end], 9.81; atol=1e-3 )
    end
end

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_PT1.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_PendulumODE.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_PendulumDAE.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_FreeBodyRotation.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_StateSelection.jl"))

end

ext: jl


text:
abstract type Functional end

quaddata(fn::Functional, refs, cells) = quadpoints(refs, cells, [5])
quadrule(fn::Functional, refs, p, cell, qd) = qd[1,p]

"""
    assemble(fn, tfs)

Assemble the vector of test coefficients corresponding to functional
`fn` and test functions `tfs`.
"""
function assemble(field::Functional, tfs)

    b = zeros(ComplexF64, numfunctions(tfs))
    store(v,m) = (b[m] += v)
    assemble!(field, tfs, store)
    return b
end

function assemble!(field::Functional, tfs::DirectProductSpace, store)
    I = Int[0]
    for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
    for (i,s) in enumerate(tfs.factors)
        store1(v,m) = store(v, m + I[i])
        assemble!(field, s, store1)
    end
end

function assemble!(field::Functional, tfs::Space, store)

    tels, tad = assemblydata(tfs)

    trefs = refspace(tfs)
    qd = quaddata(field, trefs, tels)

    for (t, tcell) in enumerate(tels)

        # compute the testing with the reference elements
        qr = quadrule(field, trefs, t, tcell, qd)
        blocal = celltestvalues(trefs, tcell, field, qr)

        for i in 1 : numfunctions(trefs)
            for (m,a) in tad[t,i]
                store(a*blocal[i], m)
            end
        end

    end

end

function assemble!(field::Functional, tfs::subdBasis, store)

    tels, tad = assemblydata(tfs)

    trefs = refspace(tfs)
    qd = quaddata(field, trefs, tels)

    for (t, tcell) in enumerate(tels)

        # compute the testing with the reference elements
        qr = quadrule(field, trefs, t, tcell, qd)
        blocal = celltestvalues(trefs, tcell, field, qr)

        for i in 1 : length(tad[t])
            for (m,a) in tad[t][i]
                store(a*blocal[i], m)
            end
        end

    end

end

function celltestvalues(tshs, tcell, field, qr)

    num_tshs = numfunctions(tshs)
    interactions = zeros(ComplexF64, num_tshs)

    num_oqp = length(qr)

    for p in 1 : num_oqp
        mp = qr[p].point

        dx =qr[p].weight

        fval = field(mp)
        tvals = qr[p].value

        for m in 1 : num_tshs
            tval = tvals[m]

            igd = integrand(field, tval, fval)
            interactions[m] += igd * dx
        end
    end

    return interactions
end

function celltestvalues(tshs::subReferenceSpace, tcell, field, qr)

    num_oqp = length(qr)
    num_tshs = length(qr[1].value[1])
    interactions = (ComplexF64, num_tshs)
    for p in 1 : num_oqp
        mp = qr[p].point

        dx =qr[p].weight

        fval = field(mp)
        tvals = qr[p].value

        interactions = zeros(ComplexF64, num_tshs)
        for m in 1 : num_tshs
            tval = tvals[1][m]

            igd = integrand(field, tval, fval)
            interactions[m] += igd * dx
        end
    end

    return interactions
end


text:
###############################################################################
#
# Preview data
#
###############################################################################
"""
    preview_data(
                xs::Array,
                ys::Array;
                title = randstring(10),
                figsize::Tuple{Number,Number} = (4,3),
                xlab::String = "X label",
                ylab::String = "Y label",
                marker::String = "",
                linestyle::String = "-",
                label_fontsize::Int = 12
    )

Preview data, given
- `xs` Array of x
- `ys` Array of y
- `title` Axis titile
- `figsize` Canvas size
- `xlab` X axis label
- `ylab` Y axis label
- `marker` Marker style
- `linestyle` Line style
- `label_fontsize` X and Y axes label font size
"""
function preview_data(
            xs::Array,
            ys::Array;
            title = randstring(10),
            figsize::Tuple{Number,Number} = (4,3),
            xlab::String = "X label",
            ylab::String = "Y label",
            marker::String = "",
            linestyle::String = "-",
            label_fontsize::Int = 12
)
    # create canvas and axis
    fig,axs = create_canvas(title; figsize=figsize);
    ax0     = axs[1];

    # plot the data
    ax0.plot(xs, ys; marker=marker, linestyle=linestyle);
    set_xylabels!(axs, xlab, ylab, fontsize=label_fontsize);
    fig.set_tight_layout(true);

    return fig
end

###############################################################################
#
# Preview NC and TIFF files
#
###############################################################################
"""
    preview_dataset!(ax::PyObject, filename::String, label)
    preview_dataset!(
                ax::PyObject,
                filename::String,
                label::String,
                format::FormatNC
    )

Preview dataset, given
- `filename` Dataset file to preview
- `label` Label of the data to preview, variable name in NC files, band name in
    TIFF files
- `format` [`AbstractFormat`](@ref) type file format
"""
function preview_dataset!(ax::PyObject, filename::String, label)
    if filename[end-2:end]==".nc" && typeof(label)==String
        preview_dataset!(ax, filename, label, FormatNC());
    else
        @warn "File format and label type is not supported!"
    end

    return nothing
end

function preview_dataset!(
            ax::PyObject,
            filename::String,
            label::String,
            format::FormatNC
)
    data = read_nc(filename, label);
    dims = size(data);

    # plot figure
    if length(dims)==2
        @info "Dataset is 2D, plotting it using imshow...";
        ax.imshow(data', origin="lower");
    elseif length(dims)==3
        @info "Dataset is 3D, plotting the 1st layer using imshow...";
        ax.imshow(data[:,:,1]', origin="lower");
    end

    return nothing
end


text:
using Random, Distributions, HypothesisTests, Plots; pyplot()

muRange = 51:0.02:55
n = 20
N = 10^4
mu0 = 53.0
powerT, powerU = [], []

for muActual in muRange

    dist = Normal(muActual, 1.2)
    rejectT, rejectU = 0, 0
    Random.seed!(1)

    for _ in 1:N
        data = rand(dist,n)
        xBar, stdDev = mean(data), std(data)

        tStatT = (xBar - mu0)/(stdDev/sqrt(n))
        pValT  = 2*ccdf(TDist(n-1), abs(tStatT))

        xPositive = sum(data .> mu0)
        uStat = max(xPositive, n-xPositive)
        pValSign = 2*cdf(Binomial(n,0.5), uStat)

        rejectT += pValT < 0.05
        rejectU += pValSign < 0.05
    end

    push!(powerT, rejectT/N)
    push!(powerU, rejectU/N)

end

plot(muRange, powerT, c=:blue, label="t test")
plot!(muRange, powerU, c=:red, label="Sign test",
	xlims=(51,55), ylims=(0,1),
	xlabel="Different values of muActual",
	ylabel="Proportion of times H0 rejected", legend=:bottomleft)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
abstract type MwlsObject end

"""
# User provided member variables
- `inputs::Array{Float64, 2}`: a 2d array of input points where each point is in a single column,
- `outputs::Array{Float64, N}`: a 2d array or a vector of output points where each output is in a single column,
- `EPS::Float64`: ε of the method (cell edge length and the default distance for range search),
- `weightfunc::Function`: weighting function θ of the method.

# Automatically created member variables
If created by the `mwls_naive` function, the member variables in this section are created automatically.

- `vars::Vector{PolyVar{true}}`: variables of the polynomial,
- `b::Vector{Monomial{true}}`: the basis of the polynomial,
- `matrix`: the result of `b * transpose(b)`.
"""
struct MwlsNaiveObject <: MwlsObject
  inputs::Array{Real, 2}
  outputs::Array{Real, N} where {N} # accepts both scalar and vector outputs
  EPS::Real
  weightfunc::Function
  vars::Vector{PolyVar{true}}
  b::Vector{Monomial{true}}
  matrix::Array{Term{true, Int}, 2}
end

function MwlsNaiveObject(inputs, outputs, EPS, weightfunc, vars, b)
  t = transpose(inputs)
  o = transpose(outputs)
  return MwlsNaiveObject(t, o, EPS, weightfunc, vars, b, b * transpose(b))
end

"""
Documentation is left out on purpose to discourage the use of MwlsNaiveObject.
K-d trees are always faster and cell linked lists are faster on datasets larger than 20.
`mwls_naive` has the same signature as `mwls_cll`.
"""
function mwls_naive(inputs::Array{T, N}, outputs::Array{U, M},
                    EPS::Real, weightfunc::Function;
                    maxDegree::Integer = 2) where {T <: Real, U <: Real, N, M}
  size(outputs, 1) != size(inputs, 1) &&
    error("The amount of inputs and outputs differs")

  vars::Vector{PolyVar{true}}, b :: Vector{Monomial{true}} =
    genmonomials(size(inputs, 2), maxDegree)

  return MwlsNaiveObject(inputs, outputs, EPS, weightfunc, vars, b)
end

"""
Documentation is left out on purpose to discourage the use of MwlsNaiveObject.
K-d trees are always faster and cell linked lists are faster on datasets larger than 20.
`mwls_naive` has the same signature as `mwls_cll`.
"""
function mwls_naive(input::Array{T, 2}, EPS::Real, weightfunc::Function;
                    outputDim::Integer = 1, maxDegree::Integer = 2) where {T <: Real}
  width = size(input, 2)
  inputEnd = width - outputDim
  outputStart = inputEnd + 1
  return mwls_naive(input[:, inputEnd == 1 ? 1 : 1:inputEnd],
                    input[:, outputStart == width ? width : outputStart:width],
                    EPS, weightfunc,
                    maxDegree = maxDegree)
end

mwlsNaive(a...;b...) = warn("`mwlsNaive` is deprecated, use `mwls_naive` instead.")

"""
# User provided member variables
- `inputs::Array{Float64, 2}`: a 2d array of input points where each point is in a single column,
- `outputs::Array{Float64, N}`: a 2d array or a vector of output points where each output is in a single column,
- `EPS::Float64`: ε of the method (cell edge length and the default distance for range search),
- `weightfunc::Function`: weighting function θ of the method.

# Automatically created member variables
If created by the `mwls_cll` function, the member variables in this section are created automatically.

- `vars::Vector{PolyVar{true}}`: variables of the polynomial,
- `b::Vector{Monomial{true}}`: the basis of the polynomial,
- `matrix`: the result of `b * transpose(b)`,
- `tree`: a k-d tree for range search.
"""
struct MwlsKdObject <: MwlsObject
  inputs::Array{Real, 2}
  outputs::Array{Real, N} where {N} # accepts both scalar and vector outputs
  EPS::Real
  weightfunc::Function
  vars::Vector{PolyVar{true}}
  b::Vector{Monomial{true}}
  matrix::Array{Term{true, Int}, 2}
  tree::KDTree
end

function MwlsKdObject(inputs, outputs, EPS, weightfunc, vars, b;
                      leafSize::Integer = 10)
  t = transpose(inputs)
  o = transpose(outputs)
  return MwlsKdObject(t, o, EPS, weightfunc,
                      vars, b, b * transpose(b),
                      KDTree(t, leafsize = leafSize))
end

"""
    mwls_kd(inputs::Array{T, N}, outputs::Array{U}, EPS::Real, weightfunc::Function) where {T <: Real, U <: Real, N}
    mwls_kd(inputs::Array{T, N}, outputs::Array{U}, EPS::Real, weightfunc::Function; leafsize::Int = 10, maxDegree::Int = 2) where {T <: Real, U <: Real, N}

Creates `MwlsKdObject` from sample input and sample output data, the cutoff distance ε and a weighting function θ.

# Arguments
- `inputs`: a 2d array or a vector of input points where each point is on a single row,
- `outputs`: a 2d array or a vector of output scalars where each output is on a single row,
- `EPS::Real`: ε of the method (default distance threshold for neighbor search),
- `weightfunc::Function`: weighting function of the method. It should be in form `(distance, EPS) -> Float64`.

# Keyword arguments
- `leafSize::Int`: the size of the leaves in the k-d-tree, 10 by default.
- `maxDegree::Int`: the maximal degree of polynomials used for approximation, 2 by default.
"""
function mwls_kd(inputs::Array{T, N}, outputs::Array{U},
                 EPS::Real, weightfunc::Function;
                 leafSize::Integer = 10, maxDegree::Integer = 2) where {T <: Real, U <: Real, N}
  size(outputs, 1) != size(inputs, 1) &&
    error("The amount of inputs and outputs differs.")

  vars::Vector{PolyVar{true}}, b::Vector{Monomial{true}} =
    genmonomials(size(inputs, 2), maxDegree)
  return MwlsKdObject(N == 1 ? hcat(inputs, zeros(size(inputs, 1))) : inputs,
                      outputs, EPS, weightfunc,
                      vars, b, leafSize = leafSize)
end

"""
    mwls_kd(input::Array{T, 2}, EPS::Real, weightfunc::Function) where {T <: Real}
    mwls_kd(input::Array{T, 2}, EPS::Real, weightfunc::Function; outputDim::Int = 1, leafSize::Int = 10, maxDegree::Int = 2) where {T <: Real}

In this `mwls_kd` function, the sample input and sample output data are passed in a single array.
It is assumed that each pair of input and output is on a single row.
Dimension of the output is specified with kwarg `outputDim`.
"""
function mwls_kd(input::Array{T, 2}, EPS::Real, weightfunc::Function;
                outputDim::Integer = 1, leafSize::Integer = 10, maxDegree::Int = 2) where {T <: Real}
  width = size(input, 2)
  inputEnd = width - outputDim
  outputStart = inputEnd + 1
  return mwls_kd(input[:, inputEnd == 1 ? 1 : 1:inputEnd],
                 input[:, outputStart == width ? width : outputStart:width],
                 EPS, weightfunc,
                 leafSize = leafSize, maxDegree = maxDegree)
end

mwlsKd(a...;b...) = warn("`mwlsKd` is deprecated, use `mwls_kd` instead.")

"""
# User provided member variables
- `inputs::Array{Float64, 2}`: a 2d array of input points where each point is in a single column,
- `outputs::Array{Float64, N}`: a 2d array or a vector of output points where each output is in a single column,
- `EPS::Float64`: ε of the method (cell edge length and the default distance for range search),
- `weightfunc::Function`: weighting function θ of the method.

# Automatically created member variables
If created by the `mwls_cll` function, the member variables in this section are created automatically.

- `vars::Vector{PolyVar{true}}`: variables of the polynomial,
- `b::Vector{Monomial{true}}`: the basis of the polynomial,
- `matrix`: the result of `b * transpose(b)`,
- `cll`: a cell linked list for range search.
"""
struct MwlsCllObject <: MwlsObject
  inputs::Array{Real, 2}
  outputs::Array{Real, N} where {N} # accepts both scalar and vector outputs
  EPS::Real
  weightfunc::Function
  vars::Vector{PolyVar{true}}
  b::Vector{Monomial{true}}
  matrix::Array{Term{true, Int}, 2}
  cll::CellLinkedList
end

function MwlsCllObject(inputs, outputs, EPS, weightfunc, vars, b)
  t = transpose(inputs)
  o = transpose(outputs)
  return MwlsCllObject(t, o, EPS, weightfunc,
                       vars, b, b * transpose(b),
                       CellLinkedList(t, EPS))
end

"""
    mwls_cll(inputs::Array{T, N}, outputs::Array{U}, EPS::Real, weightfunc::Function) where {T <: Real, U <: Real, N}
    mwls_cll(inputs::Array{T, N}, outputs::Array{U}, EPS::Real, weightfunc::Function; maxDegree::Int = 2) where {T <: Real, U <: Real, N}

Creates `MwlsCllObject` from sample input and sample output data, the cutoff distance ε and a weighting function θ.

# Arguments
- `inputs`: a 2d array of input points where each point is on a single row,
- `outputs`: a 2d array or a vector of output scalars where each output is on a single row,
- `EPS::Real`: ε of the method (cell edge length and the default distance for range search),
- `weightfunc::Function`: weighting function θ of the method. It should be in form `(distance between two vectors, EPS) -> Float64`.

# Keyword arguments
- `maxDegree::Int`: the maximal degree of polynomials used for approximation, 2 by default.
"""
function mwls_cll(inputs::Array{T, N}, outputs::Array{U},
                  EPS::Real, weightfunc::Function;
                  maxDegree::Integer = 2) where {T <: Real, U <: Real, N}
  size(outputs, 1) != size(inputs, 1) &&
    error("The amount of inputs and outputs differs.")

  vars::Vector{PolyVar{true}}, b::Vector{Monomial{true}} =
    genmonomials(size(inputs, 2), maxDegree)

  return MwlsCllObject(N == 1 ? hcat(inputs, zeros(size(inputs,1))) : inputs,
                       outputs, EPS, weightfunc, vars, b)
end

"""
    mwls_cll(input::Array{T, 2}, EPS::Real, weightfunc::Function) where {T <: Real}
    mwls_cll(input::Array{T, 2}, EPS::Real, weightfunc::Function; outputDim::Int = 1, maxDegree::Int = 2) where {T <: Real}

In this `mwls_cll` function, the sample input and sample output data are passed in a single array.
It is assumed that each pair of input and output is on a single row.
Dimension of the output is specified with kwarg `outputDim`.
"""
function mwls_cll(input::Array{T, 2}, EPS::Real, weightfunc::Function;
                  outputDim::Integer = 1, maxDegree::Integer = 2) where {T <: Real}
  width = size(input, 2)
  inputEnd = width - outputDim
  outputStart = inputEnd + 1
  return mwls_cll(input[:, inputEnd == 1 ? 1 : 1:inputEnd],
                  input[:, outputStart == width ? width : outputStart:width],
                  EPS, weightfunc, maxDegree = maxDegree)
end

mwlsCll(a...;b...) = warn("`mwlsCll` is deprecated, use `mwls_cll` instead.")


text:
using EngEconomics

function benefitCostRatio(EUAB, EUAC)
	if EUAB / EUAC >= 1.0
		println("Acceptable")
	end

	return EUAB / EUAC
end

P = 2000000
A1 = 100000
A2 = 820000
A3 = 400000
n = 20
i = 0.08

costs = P * capitalRecoveryFactor(i, n)
benefits = A2
disbenefits = A3
O_AND_M_Costs = A1

EUAB = benefits - disbenefits - O_AND_M_Costs
EUAC = costs

bcRatio = benefitCostRatio(EUAB, EUAC)


text:
function safeRound(num)
  if VERSION < v"0.4-"
    iround(num)
  else
    round(Integer,num)
  end
end

function safeFloor(num)
  if VERSION < v"0.4-"
    ifloor(num)
  else
    floor(Integer,num)
  end
end

function safeCeil(num)
  if VERSION < v"0.4-"
    iceil(num)
  else
    ceil(Integer,num)
  end
end

ceilNegLog10{F<:(@compat AbstractFloat)}(x::F) = safeCeil(-log10(x))
roundNegLog10{F<:(@compat AbstractFloat)}(x::F) = safeRound(-log10(x))
roundUpToTick{F<:(@compat AbstractFloat),R<:Real}(x::F,m::R) = x == 0. ? 0.: (x > 0 ? ceil(x, ceilNegLog10(m)) : -floor(-x, ceilNegLog10(m)))
roundDownToTick{F<:(@compat AbstractFloat),R<:Real}(x::F,m::R) = x == 0. ? 0.: (x > 0 ? floor(x, ceilNegLog10(m)) : -ceil(-x, ceilNegLog10(m)))
roundUpToSubTick{F<:(@compat AbstractFloat),R<:Real}(x::F,m::R) = x == 0. ? 0.: (x > 0 ? ceil(x, ceilNegLog10(m)+1) : -floor(-x, ceilNegLog10(m)+1))
roundDownToSubTick{F<:(@compat AbstractFloat),R<:Real}(x::F,m::R) = x == 0. ? 0.: (x > 0 ? floor(x, ceilNegLog10(m)+1) : -ceil(-x, ceilNegLog10(m)+1))
floatRoundLog10{F<:(@compat AbstractFloat),R<:Real}(x::F,m::R) = x == 0. ? 0.: (x > 0 ? round(x, ceilNegLog10(m)+1) : -round(-x, ceilNegLog10(m)+1))
floatRoundLog10{F<:(@compat AbstractFloat)}(x::F) = x > 0 ? floatRoundLog10(x,x): floatRoundLog10(x,-x)

function plottingRange{F<:(@compat AbstractFloat),R<:(@compat AbstractFloat)}(xmin::F, xmax::R)
  diffX = xmax - xmin
  xmax = roundUpToTick(xmax, diffX)
  xmin = roundDownToTick(xmin, diffX)
  xmin, xmax
end

function plottingRangeNarrow{F<:(@compat AbstractFloat),R<:(@compat AbstractFloat)}(xmin::F, xmax::R)
  diffX = xmax - xmin
  xmax = roundUpToSubTick(xmax, diffX)
  xmin = roundDownToSubTick(xmin, diffX)
  xmin, xmax
end

const borderMap = Dict{Symbol,Dict{Symbol,(@compat AbstractString)}}()
const borderSolid = Dict{Symbol,(@compat AbstractString)}()
borderSolid[:tl]="┌"
borderSolid[:tr]="┐"
borderSolid[:bl]="└"
borderSolid[:br]="┘"
borderSolid[:t]="─"
borderSolid[:l]="│"
borderSolid[:b]="─"
borderSolid[:r]="│"
const borderBold = Dict{Symbol,(@compat AbstractString)}()
borderBold[:tl]="┏"
borderBold[:tr]="┓"
borderBold[:bl]="┗"
borderBold[:br]="┛"
borderBold[:t]="━"
borderBold[:l]="┃"
borderBold[:b]="━"
borderBold[:r]="┃"
const borderNone = Dict{Symbol,(@compat AbstractString)}()
borderNone[:tl]=" "
borderNone[:tr]=" "
borderNone[:bl]=" "
borderNone[:br]=" "
borderNone[:t]=" "
borderNone[:l]=" "
borderNone[:b]=" "
borderNone[:r]=" "
const borderDashed = Dict{Symbol,(@compat AbstractString)}()
borderDashed[:tl]="┌"
borderDashed[:tr]="┐"
borderDashed[:bl]="└"
borderDashed[:br]="┘"
borderDashed[:t]="╌"
borderDashed[:l]="│"
borderDashed[:b]="╌"
borderDashed[:r]="│"
const borderDotted = Dict{Symbol,(@compat AbstractString)}()
borderDotted[:tl]="⡤"
borderDotted[:tr]="⢤"
borderDotted[:bl]="⠓"
borderDotted[:br]="⠚"
borderDotted[:t]="⠤"
borderDotted[:l]="⡇"
borderDotted[:b]="⠒"
borderDotted[:r]="⢸"
borderMap[:solid]=borderSolid
borderMap[:bold]=borderBold
borderMap[:none]=borderNone
borderMap[:dashed]=borderDashed
borderMap[:dotted]=borderDotted

const autoColors = [:blue, :red, :yellow, :magenta, :green, :cyan]

const colorEncode = Dict{Symbol,(@compat UInt8)}()
colorEncode[:white]=0b000
colorEncode[:blue]=0b001
colorEncode[:red]=0b010
colorEncode[:magenta]=0b011
colorEncode[:yellow]=0b100
colorEncode[:green]=0b101
colorEncode[:cyan]=0b110
const colorDecode = Dict{(@compat UInt8),Symbol}()
for k in keys(colorEncode)
  v = colorEncode[k]
  colorDecode[v]=k
end
colorDecode[0b111]=:white

# ▖▗▘▙▚▛▜▝▞▟
# ▁▂▃▄▅▆▇█
# ░▒▓█
# ⬛

function printColor(color::(@compat UInt8), io::IO, args...)
  #if isa(io, Base.TTY)
    col = colorDecode[color]
    str = string(args...)
    print_with_color(col, io, str)
  #else
  #  print(io, args...)
  #end
end


text:
module BasisSpecification

using ..GraphSignal, ..GraphPartition, LinearAlgebra

export BasisSpec #, bsfull, bs_haar, bs_level

# BasisSpec data structure and constructors
"""
    BS = BasisSpec(dvec, dvec_loc, c2f, description)

is a data structure for a BasisSpec object contaning the following fields:
* `levlist::Vector{Tuple{Int, Int}}`: the integer sequence that specifies a particular basis
* `c2f::Bool`: if true (default), this indicates that the basis comes from the coarse-to-fine dictionary; if false, this indicates that the basis comes from the fine-to-coarse dictionary
* `description::String`: a description of the specified basis

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito) |
Translated to julia and revised by Naoki Saito, Feb. 22, 2017
Modified by Yiqun Shao, May. 20, 2018
"""
mutable struct BasisSpec
    levlist::Vector{Tuple{Int, Int}}
    c2f::Bool
    description::String

    # An inner constructor here.
    function BasisSpec(levlist::Vector{Tuple{Int, Int}};
                       c2f::Bool = true,
                       description::String = "")
        new(levlist, c2f, description )
    end # of BasisSpec constructor

end # of type BasisSpec

end  # of module BasisSpecification


text:
# Testing partitionize
for n = 5:20
    for r = 5:10
        p = rand(1:n,r)
        d = countmap(p)
        v = sort(collect(values(d)),rev = true)
        @test HyperModularity.partitionize(p) == v
    end
end

v = vec([3 4 5 0 0])
@test HyperModularity.sortedremovezeros(v) == [3;4;5]

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
using Test, WAVI

@testset "Grids" begin
    @info "Testing Grids...."

    @testset "Testing Grid construction" begin
    @info "Testing Grid construction..."

    #no arguments
    nx = 80
    ny = 10
    grid = Grid(nx = nx, ny = ny)
    @test grid isa Grid

    #test that nothing default arguments picked up
    @test grid.h_mask == trues(nx,ny)
    @test grid.u_iszero == falses(nx+1,ny)
    @test grid.v_iszero == falses(nx,ny+1)
    @test grid.quadrature_weights ==  [0.5;ones(grid.nσ-2);0.5]/(grid.nσ-1)

    #test the string boundary conditions method
    grid = Grid(nx = nx, ny = ny, u_iszero = ["north"])
    @test all(grid.u_iszero[1,:] .== 1)
    grid = Grid(nx = nx, ny = ny, u_iszero = ["south"])
    @test all(grid.u_iszero[end,:] .== 1)
    grid = Grid(nx = nx, ny = ny, u_iszero = ["west"])
    @test all(grid.u_iszero[:,1] .== 1)
    grid = Grid(nx = nx, ny = ny, u_iszero = ["east"])
    @test all(grid.u_iszero[:,end] .== 1)
    grid = Grid(nx = nx, ny = ny, v_iszero = ["North"]) #with upper case
    @test all(grid.v_iszero[1,:] .== 1)
    grid = Grid(nx = nx, ny = ny, v_iszero = ["SouTh"])
    @test all(grid.v_iszero[end,:] .== 1)
    grid = Grid(nx = nx, ny = ny, v_iszero = ["WeSt"])
    @test all(grid.v_iszero[:,1] .== 1)
    grid = Grid(nx = nx, ny = ny, v_iszero = ["East"])
    @test all(grid.v_iszero[:,end] .== 1)
    end

    @testset "Testing Grid construction" begin
        @info "Testing Grid construction errors..."

        #h_mask size incorrect
        @test_throws DimensionMismatch Grid(nx = 10, ny = 10, h_mask = trues(5,5))
        @test_throws DimensionMismatch Grid(nx = 5, ny = 10, h_mask = trues(5,5))
        @test_throws DimensionMismatch Grid(nx = 10, ny = 5, h_mask = trues(5,5))
        nσ = 10
        @test_throws DimensionMismatch Grid(nσ = nσ, quadrature_weights =  [0.5;ones(nσ-1);0.5]/(nσ-1)) #quadrature weights one short

        #u_iszero size incorrect
        @test_throws DimensionMismatch Grid(nx = 10, ny = 10, u_iszero = ones(10,10))

        #v_iszero size incorrect
        @test_throws DimensionMismatch Grid(nx = 10, ny = 10, v_iszero = ones(10,10))

        #h_mask non-boolean
        @test_throws ArgumentError Grid(nx = 10, ny = 10, h_mask = 2.0*ones(10,10))

        #u_iszero non-boolean
        @test_throws ArgumentError Grid(nx = 10, ny = 10, u_iszero = 2.0*ones(11,10))

        #v_iszero non-boolean
        @test_throws ArgumentError Grid(nx = 10, ny = 10, v_iszero = 2.0*ones(10,11))

        #non-positive integer number of grid points
        @test_throws ArgumentError Grid(nx = 10.,)
        @test_throws ArgumentError Grid(ny = 10.,)
        @test_throws ArgumentError Grid(nσ = 10.,)
        @test_throws ArgumentError Grid(nx = -5)
        @test_throws ArgumentError Grid(ny = -5)
        @test_throws ArgumentError Grid(nσ = -5)

    end
end


text:
module SSH

    using DataStructures
    include("constants.jl")

    type Session
        transport::IO
        is_client::Bool
        mac_length::UInt8
        sequence_number::UInt32
        recv_sequence_number::UInt32
        block_size::UInt16
        encrypt!::Any
        decrypt!::Any
        hmac::Any
        channels::Vector{Any}
        allocated_channels::DataStructures.IntSet
        V_C::String
        V_S::String
        I_C::Vector{UInt8}
        I_S::Vector{UInt8}
        session_id::Vector{UInt8}
        Session(transport, is_client) = new(transport, is_client, 0, 0, 0, 8, x->x, x->x, (data,seqno)->UInt8[],
            Vector{Any}(), DataStructures.IntSet())
    end

    type Channel
        session::Session
        isopen::Bool
        input_buffer::IOBuffer
        data_available::Condition
        remote_number::UInt32
        local_number::UInt32
        window_size::UInt32
        max_packet_size::UInt32
        on_channel_request::Any
    end

    immutable PacketBuffer
        buf::IOBuffer
        function PacketBuffer(packet_type)
            buf = IOBuffer()
            write(buf, UInt8(packet_type))
            new(buf)
        end
    end
    Base.copy(buf::PacketBuffer) = copy(buf.buf)

    Base.write(buf::PacketBuffer, args...) = write(buf.buf, args...)

    our_ident = "SSH-2.0-SSH.jlv0.1 Do not use in production systems"
    function Base.connect(::Type{Session}, transport::IO; client = true)
        # Send the indentification string
        write(transport, our_ident, "\r\n")
        ident = readuntil(transport, "\r\n")[1:end-2]
        # Some servers may announce dual-version support by sending SSH-1.99-. They
        # are not supported.
        startswith(ident, "SSH-2.0-") || error("Invalid identification string `$ident`")
        sess = Session(transport, client)
        if client
            sess.V_C = our_ident; sess.V_S = ident
        else
            sess.V_C = ident; sess.V_S = our_ident
        end
        sess
    end

    block_size(session) = session.block_size
    function read_packet(session::Session)
        encrypted_begin = read(session.transport, block_size(session))
        buf = IOBuffer(session.decrypt!(encrypted_begin), true, true)
        packet_length = bswap(read(buf, UInt32))
        padding_length = read(buf, UInt8)
        remaining_enc_length = (packet_length + sizeof(UInt32)) - block_size(session)
        remaining = session.decrypt!(read(session.transport, remaining_enc_length))
        pos = position(buf); seekend(buf); write(buf, remaining); seek(buf, pos)
        payload_length = packet_length - padding_length - sizeof(UInt8)
        payload = read(buf, payload_length)
        skip(buf, padding_length)
        # mac is not encrypted
        mac = read(session.transport, session.mac_length)
        # TODO: Verify MAC here
        session.recv_sequence_number += 1
        IOBuffer(payload)
    end

    function Base.write(session::Session, buf::PacketBuffer)
        payload = takebuf_array(copy(buf.buf))
        block_size = session.block_size
        packet_length = sizeof(payload) + sizeof(UInt8)
        padding_size = mod(block_size-mod(packet_length+sizeof(UInt32), block_size), block_size)
        if padding_size < 4
            padding_size += block_size
        end
        packet_length += padding_size
        # Blackhole writes if the connection is closed, but still allow reads
        # to be processed
        if isopen(session.transport)
            buf = IOBuffer()
            write(buf, bswap(UInt32(packet_length)))
            write(buf, UInt8(padding_size))
            write(buf, payload)
            write(buf, rand(UInt8, padding_size))
            unencrypted_data = takebuf_array(buf)
            hmac = session.hmac(unencrypted_data, session.sequence_number)
            encrypted = session.encrypt!(unencrypted_data)
            sendbuf = IOBuffer()
            write(sendbuf, encrypted); write(sendbuf, hmac)
            # Make sure this is task-atomic
            session.sequence_number += 1
            write(session.transport, takebuf_array(sendbuf))
        end

        # mac
    end

    function require_packet(session::Session)
        local packet
        while true
            packet = read_packet(session)
            if Base.peek(packet) == SSH_MSG_IGNORE
                continue
            end
            break
        end
        packet
    end
    function require_packet(session, kind)
        packet = require_packet(session)
        got_kind = read(packet, UInt8)
        if got_kind != kind
            error("Expected package $kind got $got_kind")
        end
        packet
    end

    const kex_algorithms = ["diffie-hellman-group14-sha1"]
    const server_host_key_algorithms = ["ssh-rsa"]
    const encryption_algorithms = ["aes128-ctr"]
    const mac_algorithms = ["hmac-sha1"]
    const compression_algorithms = ["zlib","none"]

    function write_name_list(io, list)
        nl = join(list, ",")
        write(io, bswap(UInt32(sizeof(nl))))
        write(io, nl)
    end

    function send_kexinit!(session)
        packet = PacketBuffer(SSH_MSG_KEXINIT)
        write(packet, rand(UInt8, 16))
        write_name_list(packet, kex_algorithms)
        write_name_list(packet, server_host_key_algorithms)
        write_name_list(packet, encryption_algorithms)
        write_name_list(packet, encryption_algorithms)
        write_name_list(packet, mac_algorithms)
        write_name_list(packet, mac_algorithms)
        write_name_list(packet, compression_algorithms)
        write_name_list(packet, compression_algorithms)
        write_name_list(packet, [])
        write_name_list(packet, [])
        write(packet, UInt8(0)) # first_kex_packet_follows
        write(packet, UInt32(0))
        if session.is_client
            session.I_C = takebuf_array(copy(packet))
        else
            session.I_S = takebuf_array(copy(packet))
        end
        write(session, packet)
    end

    function negotiate_algorithm(packet, client_list; allow_none = false)
        list = read_name_list(packet)
        idx = findfirst(alg->alg in list, client_list)
        !allow_none && idx == 0 && error("Could not negotiate")
        idx == 0 ? "" : client_list[idx]
    end

    function read_name_list(packet)
        nbytes = bswap(read(packet, UInt32))
        data = read(packet,nbytes)
        Set(split(String(data),','))
    end
    function negotiate_algorithms!(session)
        send_kexinit!(session)
        packet = require_packet(session, SSH_MSG_KEXINIT)
        if session.is_client
            session.I_S = takebuf_array(copy(packet))
        else
            session.I_C = takebuf_array(copy(packet))
        end
        skip(packet, 16) # Cookie
        remote_kex_algorithms = read_name_list(packet)
        remote_host_key_algorithms = read_name_list(packet)
        kex_algorithm = findfirst(alg->alg in remote_kex_algorithms, kex_algorithms)
        cs_crypt_algorithm = negotiate_algorithm(packet, encryption_algorithms)
        sc_crypt_algorithm = negotiate_algorithm(packet, encryption_algorithms)
        cs_mac_algorithm = negotiate_algorithm(packet, mac_algorithms)
        sc_mac_algorithm = negotiate_algorithm(packet, mac_algorithms)
        cs_comp_algorithm = negotiate_algorithm(packet, compression_algorithms)
        sc_comp_algorithm = negotiate_algorithm(packet, compression_algorithms)
        cs_lang_algorithm = negotiate_algorithm(packet, []; allow_none = true)
        sc_lang_algorithm = negotiate_algorithm(packet, []; allow_none = true)
        first_kex_packet_follows = read(packet, UInt8)
        reserved = read(packet, UInt32)
        @assert first_kex_packet_follows == false && reserved == 0
    end

    function write_mpint(packet, mpint::BigInt)
        size = ndigits(mpint, 2)
        nbytes = div(size+8-1,8)
        data = Array(UInt8, nbytes)
        count = Ref{Csize_t}(0)
        ccall((:__gmpz_export,:libgmp), Ptr{Void},
                (Ptr{Void}, Ptr{Csize_t}, Cint, Csize_t, Cint, Csize_t, Ptr{BigInt}),
                data, count, 1, 1, 1, 0, &mpint)
        @assert count[] == nbytes
        # Need extra padding
        need_padding = false
        if size % 8 == 0
            nbytes += 1
            need_padding = true
        end
        write(packet, bswap(UInt32(nbytes)))
        if need_padding
            write(packet, UInt8(0))
        end
        write(packet, data)
    end

    function import_bigint{T}(data::Vector{T})
        b = BigInt()
        ccall((:__gmpz_import,:libgmp), Ptr{Void},
            (Ptr{BigInt}, Csize_t, Cint, Csize_t, Cint, Csize_t, Ptr{Void}),
            &b, length(data), 1, sizeof(data[]), 1, 0, data)
        b
    end

    read_string(packet) = read(packet, bswap(read(packet, UInt32)))
    write_string(packet, data) = (write(packet, bswap(UInt32(sizeof(data)))); write(packet, data))
    function read_mpint(packet)
        data = read_string(packet)
        import_bigint(data)
    end

    immutable DHGroup
        g::BigInt
        p::BigInt
        x_max::BigInt
    end
    DHGroup(g::BigInt,data::Vector) = DHGroup(g, import_bigint(data), BigInt(2)^1024)

    # RFC 3526 Group 14
    const group14 = DHGroup(BigInt(2),map(bswap,[
        0xFFFFFFFF, 0xFFFFFFFF, 0xC90FDAA2, 0x2168C234, 0xC4C6628B, 0x80DC1CD1,
        0x29024E08, 0x8A67CC74, 0x020BBEA6, 0x3B139B22, 0x514A0879, 0x8E3404DD,
        0xEF9519B3, 0xCD3A431B, 0x302B0A6D, 0xF25F1437, 0x4FE1356D, 0x6D51C245,
        0xE485B576, 0x625E7EC6, 0xF44C42E9, 0xA637ED6B, 0x0BFF5CB6, 0xF406B7ED,
        0xEE386BFB, 0x5A899FA5, 0xAE9F2411, 0x7C4B1FE6, 0x49286651, 0xECE45B3D,
        0xC2007CB8, 0xA163BF05, 0x98DA4836, 0x1C55D39A, 0x69163FA8, 0xFD24CF5F,
        0x83655D23, 0xDCA3AD96, 0x1C62F356, 0x208552BB, 0x9ED52907, 0x7096966D,
        0x670C354E, 0x4ABC9804, 0xF1746C08, 0xCA18217C, 0x32905E46, 0x2E36CE3B,
        0xE39E772C, 0x180E8603, 0x9B2783A2, 0xEC07A28F, 0xB5C55DF0, 0x6F4C52C9,
        0xDE2BCBF6, 0x95581718, 0x3995497C, 0xEA956AE5, 0x15D22618, 0x98FA0510,
        0x15728E5A, 0x8AACAA68, 0xFFFFFFFF, 0xFFFFFFFF ]))

    using MbedTLS

    function update_string!(hasher, string)
        update!(hasher, reinterpret(UInt8,[bswap(UInt32(sizeof(string)))]))
        update!(hasher, string)
    end

    function mpint_arr(mpint)
        buf = IOBuffer()
        write_mpint(buf, mpint)
        takebuf_array(buf)
    end

    function compute_kex_hash(session, K_S, e_data, f_data, K)
        hasher=MbedTLS.MD(MD_SHA1)
        write_string(hasher, session.V_C)
        write_string(hasher, session.V_S)
        write_string(hasher, session.I_C)
        write_string(hasher, session.I_S)
        write_string(hasher, K_S)
        write_string(hasher, e_data)
        write_string(hasher, f_data)
        write_mpint(hasher, K)
        MbedTLS.finish!(hasher)
    end

    function setup_crypt!(session, K, H)
        K_data = mpint_arr(K)
        # Key derivation (rfc4253 - 7.2)
        # HASH(K || H || X || session_id)
        function derive_key(X)
            hasher=MbedTLS.MD(MD_SHA1)
            write(hasher, K_data)
            write(hasher, H)
            write(hasher, UInt8[X])
            write(hasher, session.session_id)
            MbedTLS.finish!(hasher)
        end
        cipher_type = MbedTLS.CipherInfo(MbedTLS.CIPHER_AES_128_CTR)
        block_size = 16
        cs_IV = derive_key('A')[1:block_size]
        sc_IV = derive_key('B')[1:block_size]
        cs_enc = derive_key('C')[1:block_size]
        sc_enc = derive_key('D')[1:block_size]
        cs_int = derive_key('E')[1:20]
        sc_int = derive_key('F')[1:20]

        # If client, Validate signature here, I guess

        # Newkeys barrier
        write(session, PacketBuffer(SSH_MSG_NEWKEYS))
        require_packet(session, SSH_MSG_NEWKEYS)

        # All further traffic will be encrypted
        encryptor = MbedTLS.Cipher(cipher_type)
        MbedTLS.set_key!(encryptor, session.is_client ? cs_enc : sc_enc, MbedTLS.ENCRYPT)
        MbedTLS.set_iv!(encryptor, session.is_client ? cs_IV : sc_IV)
        decryptor = MbedTLS.Cipher(cipher_type)
        MbedTLS.set_key!(decryptor, session.is_client ? sc_enc : cs_enc, MbedTLS.DECRYPT)
        MbedTLS.set_iv!(decryptor, session.is_client ? sc_IV : cs_IV)
        session.encrypt! = function(data)
            MbedTLS.update!(encryptor, data, data)
            @assert MbedTLS.finish!(encryptor, data) == 0
            data
        end
        session.decrypt! = function(data)
            MbedTLS.update!(decryptor, data, data)
            @assert MbedTLS.finish!(decryptor, data) == 0
            data
        end
        session.hmac = function(data, seqno::UInt32)
            hmac = Array(UInt8, 20)
            hmacx = MbedTLS.MD(MD_SHA1, session.is_client ? cs_int : sc_int)
            write(hmacx, reinterpret(UInt8,[bswap(seqno)]))
            write(hmacx, data)
            MbedTLS.finish!(hmacx, hmac)
            hmac
        end
        session.mac_length = 20
        session.block_size = 16
    end

    function client_dh_kex!(session)
        packet = PacketBuffer(SSH_MSG_KEXDH_INIT)
        group = group14
        # Perform Diffie-Hellman key exchange
        x = rand(1:group.x_max)
        e = powermod(group.g, x, group.p)
        e_data = mpint_arr(e)
        write(packet, e_data)
        write(session, packet)
        response = require_packet(session, SSH_MSG_KEXDH_REPLY)
        K_S = read(response, bswap(read(response, UInt32)))
        f_data = read_string(response)
        f = import_bigint(f_data)
        signature = read(response, bswap(read(response, UInt32)))
        # Compute K
        K = powermod(f, x, group.p)
        K_data = mpint_arr(K)
        # Compute H
        session.session_id = H = compute_kex_hash(session, K_S, e_data[5:end], f_data, K)
        setup_crypt!(session, K, H)
    end

    function server_dh_kex!(session, hostkey_pub, hostkey_priv)
        group = group14
        # Get client's public value
        packet = require_packet(session, SSH_MSG_KEXDH_INIT)
        e_data = read_string(packet)
        e = import_bigint(e_data)

        # Prepare response
        packet = PacketBuffer(SSH_MSG_KEXDH_REPLY)

        # Get K_S from public key
        pubkeydata = read(open(hostkey_pub))
        K_S = base64decode(String(split(String(pubkeydata),' ')[2]))

        write_string(packet, K_S)

        # Compute y, f
        y = rand(0:group.x_max)
        f = powermod(group.g, y, group.p)
        f_data = mpint_arr(f)

        write(packet, f_data)

        # Compute K and H
        K = powermod(e, y, group.p)
        session.session_id = H = compute_kex_hash(session, K_S, e_data, f_data[5:end], K)

        # Sign H
        write_string(packet, generate_signature(pubkeydata, hostkey_priv, H))

        # Send back the packet
        write(session, packet)

        # Set up crypto
        setup_crypt!(session, K, H)
    end

    function disconnect(session)
        packet = PacketBuffer(SSH_MSG_DISCONNECT)
        write(packet, bswap(UInt32(SSH_DISCONNECT_BY_APPLICATION)))
        write_string(packet, "Disconnect requested")
        write_string(packet, "en")
        write(session, packet)
        close(session.transport)
    end

    # Authentication ("ssh-userauth" service) (RFC 4252)
    function enter_userauth!(session)
        packet = PacketBuffer(SSH_MSG_SERVICE_REQUEST)
        write_string(packet, "ssh-userauth")
        write(session, packet)
        @assert read_string(require_packet(session, SSH_MSG_SERVICE_ACCEPT)) == "ssh-userauth"
    end

    function wait_for_userauth(session; publickey = nothing, keyboard_interactive=nothing, success_message = true)
        @assert String(read_string(require_packet(session, SSH_MSG_SERVICE_REQUEST))) == "ssh-userauth"
        packet = PacketBuffer(SSH_MSG_SERVICE_ACCEPT)
        write_string(packet, "ssh-userauth"); write(session, packet)

        while true
            packet = require_packet(session, SSH_MSG_USERAUTH_REQUEST)
            username = String(read_string(packet))
            servicename = String(read_string(packet))
            methodname = String(read_string(packet))
            @assert servicename == "ssh-connection"
            if methodname == "none"
            elseif methodname == "publickey"
                has_sig = read(packet, UInt8) != 0
                algorithm = read_string(packet)
                blob = read_string(packet)
                if publickey !== nothing && publickey(username, algorithm, blob)
                    if !has_sig
                        # Inform the client that this public key is acceptable
                        packet = PacketBuffer(SSH_MSG_USERAUTH_PK_OK)
                        write_string(packet, algorithm)
                        write_string(packet, blob)
                        write(session, packet)
                        continue
                    else
                        signature = IOBuffer(read_string(packet))
                        sigmethod = read_string(signature)
                        @assert sigmethod == algorithm
                        sigblob = read_string(signature)
                        # Load the public key
                        blobbuf = IOBuffer(blob)
                        algname = read_string(blobbuf)
                        @assert algname == algorithm
                        e = read_mpint(blobbuf)
                        n = read_mpint(blobbuf)
                        pubkey = MbedTLS.pubkey_from_vals!(MbedTLS.RSA(
                            MbedTLS.MBEDTLS_RSA_PKCS_V15, MD_SHA1), e, n)
                        # Generate the data over which to verify the signature
                        sigbuf = IOBuffer()
                        write_string(sigbuf, session.session_id)
                        write(sigbuf, UInt8(SSH_MSG_USERAUTH_REQUEST))
                        write_string(sigbuf, username); write_string(sigbuf, servicename);
                        write_string(sigbuf, methodname); write(sigbuf, UInt8(1))
                        write_string(sigbuf, algorithm); write_string(sigbuf, blob)
                        # Verify signature (throws on failure)
                        MbedTLS.verify(pubkey, MD_SHA1,
                            MbedTLS.digest(MD_SHA1, takebuf_array(sigbuf)),
                            sigblob)
                        # Indicate authentication success
                        success_message && write(session, PacketBuffer(SSH_MSG_USERAUTH_SUCCESS))
                        break
                    end
                    # Fall through to USERAUTH_FAILURE
                end
                # Fall through to USERAUTH_FAILURE
            elseif methodname == "keyboard-interactive"
                language_tag = read_string(packet)
                submethods = read_string(packet)
                if keyboard_interactive(session, submethods)
                    success_message && write(session, PacketBuffer(SSH_MSG_USERAUTH_SUCCESS))
                    break
                end
                # Fall through to USERAUTH_FAILURE
            else
                error("Unsupported Method `$methodname`")
            end
            userauth_failure(session, publickey=publickey,keyboard_interactive=keyboard_interactive)
        end
    end
    function userauth_failure(session;publickey=nothing,keyboard_interactive=nothing,partial_success=false)
        packet = PacketBuffer(SSH_MSG_USERAUTH_FAILURE)
        avilable_methods = []
        publickey !== nothing && push!(avilable_methods, "publickey")
        keyboard_interactive !== nothing && push!(avilable_methods, "keyboard-interactive")
        write_name_list(packet, avilable_methods); write(packet, UInt8(partial_success))
        write(session, packet)
    end
    function userauth_partial_success(session;kwargs...)
        userauth_failure(session;kwargs...,partial_success=true)
        wait_for_userauth(session; kwargs..., success_message = false)
        return true
    end

    function keyboard_interactive_prompt(session, name, instruction, prompts = [], echo = true)
        packet = PacketBuffer(SSH_MSG_USERAUTH_INFO_REQUEST)
        write_string(packet, name)
        write_string(packet, instruction)
        write_string(packet, "")
        nprompts = isa(prompts, Array) ? length(prompts) : 1
        write(packet, bswap(Cint(nprompts)))
        isa(prompts, Array) && !isa(echo, Array) && (echo = fill!(Vector{Bool}(length(prompts)),echo))
        for (prompt, do_echo) in (isa(prompts, Array) ? zip(prompts,echo) : zip((prompts,),(echo,)))
            write_string(packet, prompt)
            write(packet, UInt8(do_echo))
        end
        write(session, packet)
        reply = require_packet(session, SSH_MSG_USERAUTH_INFO_RESPONSE)
        nreplys = bswap(read(reply, Cint))
        if nreplys != nprompts
            error("Reply number mismatch")
        end
        replys = map(1:nreplys) do _
            read_string(reply)
        end
        return !isa(prompts, Array) ? replys[] : replys
    end

    function clientauth_list(session, username, servicename = "ssh-connection")
        packet = PacketBuffer(SSH_MSG_USERAUTH_REQUEST)
        write_string(packet, username)
        write_string(packet, servicename)
        write_string(packet, "none")
        write(session, packet)
        fp = require_packet(session, SSH_MSG_USERAUTH_FAILURE)
        available_auth = read_name_list(fp)
        @show available_auth
    end

    function generate_signature(pubkeydata, privkey, data)
        pk = MbedTLS.parse_keyfile(privkey)
        sig = Array(UInt8, 1024)
        rng = Base.MersenneTwister()
        len = MbedTLS.sign!(pk, MD_SHA1, MbedTLS.digest(MD_SHA1, data), sig, rng)
        sig_buf = IOBuffer()
        write_string(sig_buf, pubkeydata[1:7])
        write_string(sig_buf, sig[1:len])
        takebuf_array(sig_buf)
    end

    function clientauth_pubkey(session, username, pubkey, privkey, servicename = "ssh-connection")
        packet = PacketBuffer(SSH_MSG_USERAUTH_REQUEST)

        # Open publickey
        pubkeydata = read(open(pubkey))

        # Put together request
        function write_request(io)
            write_string(io, username)
            write_string(io, servicename)
            write_string(io, "publickey")
            write(io, UInt8(1))

            write_string(io, pubkeydata[1:7])
            write_string(io, base64decode(String(split(String(pubkeydata),' ')[2])))
        end
        write_request(packet)

        # Generate rsa signature
        # Step 1: Assemble data to take the signature over
        buf = IOBuffer()
        write_string(buf, session.session_id)
        write(buf, UInt8(SSH_MSG_USERAUTH_REQUEST))
        write_request(buf)

        # Step 2: Open private key and compute signature
        write_string(packet, generate_signature(pubkeydata, privkey, takebuf_array(buf)))

        write(session, packet)
        packet = require_packet(session)
        kind = read(packet, UInt8)
        if kind == SSH_MSG_USERAUTH_SUCCESS
            return true
        elseif kind == SSH_MSG_USERAUTH_FAILURE
            return false
        else
            error(kind)
        end
    end

    # ssh-connection protocol
    Base.isopen(chan::Channel) = chan.isopen
    function Base.open(chan::Channel)
        chan.isopen = true
        packet = PacketBuffer(SSH_MSG_CHANNEL_OPEN_CONFIRMATION)
        write(packet, bswap(chan.remote_number))
        write(packet, bswap(chan.local_number))
        write(packet, bswap(chan.window_size))
        write(packet, bswap(chan.max_packet_size))
        write(chan.session, packet)
    end

    function Base.write(chan::Channel, data)
        packet = PacketBuffer(SSH_MSG_CHANNEL_DATA)
        write(packet, bswap(chan.remote_number))
        write_string(packet, data)
        write(chan.session, packet)
    end

    function Base.readavailable(chan::Channel)
        while nb_available(chan.input_buffer) == 0
            wait(chan.data_available)
        end
        readavailable(chan.input_buffer)
    end

    on_channel_request(f, chan) = chan.on_channel_request = f

    function allocate_channel_no(session)
        idx = DataStructures.nextnot(session.allocated_channels, 0)[2]
        (idx > length(session.channels)) && resize!(session.channels, idx)
        idx
    end

    function perform_ssh_connection(new_chan_cb, session::Session)
        while true
            packet = require_packet(session)
            kind = read(packet, UInt8)
            if kind == SSH_MSG_CHANNEL_OPEN
                chan_kind = read_string(packet)
                remote_no = bswap(read(packet, UInt32))
                initial_window = bswap(read(packet, UInt32))
                max_packet_size = bswap(read(packet, UInt32))
                local_no = allocate_channel_no(session)
                chan = Channel(session, false, PipeBuffer(), Condition(),
                    remote_no, local_no, initial_window, max_packet_size, nothing)
                session.channels[local_no] = chan
                new_chan_cb(String(chan_kind), chan)
                if !isopen(chan)
                    packet = PacketBuffer(SSH_MSG_CHANNEL_OPEN_FAILURE)
                    write(packet, bswap(UInt32(remote_no)))
                    write(packet, bswap(UInt32(SSH_OPEN_UNKNOWN_CHANNEL_TYPE)))
                    write_string(packet, "Unknown channel type")
                    write_string(packet, "en")
                end
            elseif kind == SSH_MSG_CHANNEL_REQUEST
                channel_no = bswap(read(packet, UInt32))
                chan = session.channels[channel_no]
                req_type = String(read_string(packet))
                if chan.on_channel_request !== nothing
                    if chan.on_channel_request(req_type, packet)
                        reply_packet = PacketBuffer(SSH_MSG_CHANNEL_SUCCESS)
                        write(reply_packet, chan.remote_number)
                        write(session, reply_packet)
                    end
                else
                    reply_packet = PacketBuffer(SSH_MSG_CHANNEL_FAILURE)
                    write(reply_packet, chan.remote_number)
                    write(session, reply_packet)
                end
            elseif kind == SSH_MSG_CHANNEL_DATA
                channel_no = bswap(read(packet, UInt32))
                chan = session.channels[channel_no]
                data = read_string(packet)
                write(chan.input_buffer, data)
                notify(chan.data_available)
            else
                @show kind
                packet = PacketBuffer(SSH_MSG_UNIMPLEMENTED)
                write(packet, bswap(UInt32(session.recv_sequence_number)))
                write(session, packet)
            end
        end
    end

    # Encoded terminal modes
    const NCCS = is_linux() ? 32 : 20
    const tcflag_t = is_linux() ? Cuint : Culong
    const speed_t = tcflag_t
    immutable termios
        c_iflag::tcflag_t
        c_oflag::tcflag_t
        c_cflag::tcflag_t
        c_lflag::tcflag_t
        @static if is_linux()
            c_line::UInt8
        end
        c_cc::NTuple{NCCS, UInt8}
        c_uispeed::speed_t
        c_ospeed::speed_t
    end

    const maps_idx = is_linux() ? 1 : 2

    op_char_map = Dict(
    #  SSH => (linux, Apple/BSD)
         1 => (  0,  8),  # VINTR
         2 => (  1,  9),  # VQUIT
         3 => (  2,  3),  # VERASE
         4 => (  3,  5),  # VKILL
         5 => (  4,  0),  # VEOF
         6 => ( 11,  1),  # VEOL
         7 => ( 16,  2),  # VEOL2
         8 => (  8, 12),  # VSTART
         9 => (  9, 13),  # VSTOP
        10 => ( 10, 10),  # VSUSP
        11 => ( -1, 11),  # VDUSP
        12 => ( 12,  6),  # VREPRINT
        13 => ( 14,  4),  # VWERASE
        14 => ( 15, 14),  # VLNEXT
        15 => ( -1, -1),  # VFLUSH
        16 => ( -1, -1),  # VSWTCH
        17 => ( -1, 18),  # VSTATUS
        18 => ( -1, 15),  # VDISCARD
    )

    iflag_map = Dict(
    #  SSH => ( Linux   , Apple/BSD )
        30 => (0o0000004, 0x00000004), # IGNPAR
        31 => (0o0000010, 0x00000008), # PARMRK
        32 => (0o0000020, 0x00000010), # INPCK
        33 => (0o0000040, 0x00000020), # ISTRIP
        34 => (0o0000100, 0x00000040), # INLCR
        35 => (0o0000200, 0x00000080), # IGNCR
        36 => (0o0000400, 0x00000100), # ICRNL
        37 => (0o0001000, 0x00000000), # IUCLC
        38 => (0o0002000, 0x00000200), # IXON
        39 => (0o0004000, 0x00000800), # IXANY
        40 => (0o0010000, 0x00000400), # IXOFF
        41 => (0o0020000, 0x00002000), # IMAXBEL
        42 => (0o0040000, 0x00004000), # IUTF8
    )

    lflag_map = Dict(
    #  SSH => ( Linux   , Apple/BSD )
        50 => (0o0000001, 0x00000080), # ISIG
        51 => (0o0000002, 0x00000100), # ICANON
        52 => (0o0000004, 0x00000000), # XCASE
        53 => (0o0000010, 0x00000008), # ECHO
        54 => (0o0000020, 0x00000002), # ECHOE
        55 => (0o0000040, 0x00000004), # ECHOK
        56 => (0o0000100, 0x00000010), # ECHONL
        57 => (0o0000200, 0x80000000), # NOFLSH
        58 => (0o0000400, 0x00400000), # TOSTOP
        59 => (0o0100000, 0x00000400), # IEXTEN
        60 => (0o0001000, 0x00000040), # ECHOCTL
        61 => (0o0004000, 0x00000001), # ECHOKE
        62 => (0o0040000, 0x20000000), # PENDIN
    )

    oflag_map = Dict(
    #  SSH => ( Linux   , Apple/BSD )
        70 => (0o0000001, 0x00000001), # OPOST
        71 => (0o0000002, 0x00000000), # OLCUC
        72 => (0o0000004, 0x00000002), # ONLCR
        73 => (0o0000010, 0x00000010), # OCRNL
        74 => (0o0000020, 0x00000020), # ONOCR
        75 => (0o0000040, 0x00000040), # ONLRET
    )

    cflag_map = Dict(
    #  SSH => ( Linux   , Apple/BSD )
        90 => (0o0000040, 0x00000200), # CS7
        91 => (0o0000060, 0x00000300), # CS8
        92 => (0o0000400, 0x00001000), # PARENB
        93 => (0o0001000, 0x00002000), # PARODD
    )

    function process_flags(flags, mask, operand)
        if operand != 0
            flags |= mask
        else
            flags &= ~mask
        end
        flags
    end
    function decode_modes(buf, old_termios)
        iflags = UInt32(old_termios.c_iflag)
        lflags = UInt32(old_termios.c_lflag)
        oflags = UInt32(old_termios.c_oflag)
        cflags = UInt32(old_termios.c_cflag)
        c_cc = zeros(UInt8, NCCS)
        for i = 1:NCCS
            c_cc[i] = old_termios.c_cc[i]
        end
        while !eof(buf)
            opcode = read(buf, UInt8)
            if 1 <= opcode <= 159
                operand = bswap(read(buf, UInt32))
                if haskey(op_char_map, opcode)
                    op_char_map[opcode][maps_idx] != -1 &&
                        (c_cc[op_char_map[opcode][maps_idx]+1] = operand == 255 ?
                         0 : operand)
                elseif haskey(iflag_map, opcode)
                    iflags |= process_flags(iflags, iflag_map[opcode][maps_idx], operand)
                elseif haskey(lflag_map, opcode)
                    lflags |= process_flags(lflags, lflag_map[opcode][maps_idx], operand)
                elseif haskey(oflag_map, opcode)
                    oflags |= process_flags(oflags, oflag_map[opcode][maps_idx], operand)
                elseif haskey(cflag_map, opcode)
                    cflags |= process_flags(cflags, cflag_map[opcode][maps_idx], operand)
                end
            else
                break
            end
        end
        return termios(iflags, oflags, cflags, lflags,
            (is_linux() ? (0,) : ())..., tuple(c_cc...), 0, 0)
    end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 2264b716-a37a-11eb-1393-7f19142f6d9c
begin
	using LinearAlgebra
	using BandedMatrices

	using Plots; Plots.gr()
	using LaTeXStrings
end

# ╔═╡ 44c0e18e-a42b-11eb-1810-3fe2e9e81d53
begin
	# Taken from https://github.com/JuliaMatrices/BandedMatrices.jl/issues/180
	function kron(a::BandedMatrix{T}, b::BandedMatrix{V}) where {T,V}
		R = BandedMatrix{Base.promote_op(*,T,V)}(undef, (size(a,1) * size(b,1), size(a,2) * size(b,2)), (bandwidth(a,1) * size(b,1) + bandwidth(b,1), bandwidth(a,2) * size(b,2) + bandwidth(b,2)))
		m = 0
		@inbounds for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
			aij = a[i,j]
			for k = 1:size(b,1)
				R[m += 1] = aij * b[k,l]
			end
		end
		R
	end

	kron(A::BandedMatrix, B::BandedMatrix, C::AbstractMatrix...) = kron_(kron_(A,B), C...)
	kron(A) = A
end

# ╔═╡ 3cec498e-a37b-11eb-0b97-e33d6bec2dea
begin
	# Simulation parameters
	const N_max = 5 	# Fock space truncation

	# Single site operators
	const bₙ = BandedMatrix(1 => (1:N_max) .^ (1/2))
	const Iₙ = BandedMatrix(0 => ones(Float64, N_max+1))

	@assert bₙ' * bₙ ≈ Iₙ .* (0:N_max)

	"""
		H(t, μ, b̄, L)

	Bose-Hubbard `L`-chain mean-field Hamiltonian
	"""
	function H(t, μ, b̄, L)
		b = [kron((Iₙ for _ in 1:(i-1))..., bₙ, (Iₙ for _ in (i+1):L)...) for i in 1:L]
		I = BandedMatrix(0 => ones(Float64, (N_max+1)^L))

		# on-site
		h = sum(-μ * b[i]' * b[i] + 0.5 * b[i]' * b[i] .* (b[i]' * b[i] .- I) for i in 1:L)

		# mean-field hopping
		h += -t * (b[1]' * b̄ + b̄' * b[1])

		if L > 1
			# mean-field hopping
			h += -t * (b[L]' * b̄ + b̄' * b[L])

			# chain hopping
			h += -t * sum(b[i+1]' * b[i] + b[i]' * b[i+1] for i in 1:(L-1))
		end

		return h
	end

	"""
		mean_field(t, μ, b̄, L)

	Bose-Hubbard mean-field at T=0 (ϕ being the ground-state)
		<ϕ|b|ϕ>
	"""
	function mean_field(t, μ, b̄, L)
		ϕ = begin
			_, vecs = eigen(Matrix(H(t, μ, b̄, L)))
			vecs[:,1] # eigen-vector of the lowest eigen-value
		end

		# If H is huge, prefer `using KrylovKit`
		# ϕ = eigsolve(H(t, μ, b̄, L), 1, :SR)[2][1]

		return ϕ' * kron(bₙ, (Iₙ for _ in 2:L)...) * ϕ
	end

	"""
		fixedpoint(f, x0; xtol)

	Finds the fixed point
		x* = f(x*)
	"""
	function fixedpoint(f, x; xtol=1e-4, β=0.1)
		x0, x = Inf, x
		while norm(x - x0, Inf) > xtol
			x0, x = x, β * x + (1 - β) * f(x)
		end
		return x
	end
end;

# ╔═╡ 80506ae8-a40e-11eb-117c-6b22718cc945
begin
	# Physical parameters (units of U)
	μs = 0:0.03:4
	ts = 0:0.006:0.2

	L = 1 	# chain size

	# Parameter scan
	data = [fixedpoint(b̄ -> mean_field(t, μ, b̄, L), 1.0) for μ in μs, t in ts]
end;

# ╔═╡ 3958ae48-a37d-11eb-26f7-b7511953016e
heatmap(ts, μs, data;
	xlabel=L"t/U",
	ylabel=L"\mu/U",
	title=L"\textrm{mean-field} <b>",
	xlims=extrema(ts))

# ╔═╡ Cell order:
# ╠═2264b716-a37a-11eb-1393-7f19142f6d9c
# ╟─44c0e18e-a42b-11eb-1810-3fe2e9e81d53
# ╠═3cec498e-a37b-11eb-0b97-e33d6bec2dea
# ╠═80506ae8-a40e-11eb-117c-6b22718cc945
# ╠═3958ae48-a37d-11eb-26f7-b7511953016e

max_stars_repo_name: FMeirinhos/unsorted-notebooks

max_issues_repo_name: FMeirinhos/unsorted-notebooks

max_forks_repo_name: FMeirinhos/unsorted-notebooks


text:
"""
    getbuiltins(ssa::Array)

Get an array of strings with any used builtin .wat functions as specified by the dict builtinfuncs.
"""
function getbuiltins(ssa::Array)
    builtins = Dict()
    builtins!(builtins, ssa)
    return values(builtins)
end

builtins!(builtins, item) = isbuiltin(item) ? builtins[item.name] = builtinfuncs[item.name] : nothing
builtins!(builtins, items::Array) = builtins!.((builtins,), items)

isbuiltin(item) = false
isbuiltin(item::GlobalRef) = haskey(builtinfuncs, item.name)

"""
builtinfuncs: a Dict with handwritten .wat of some julia builtins.
"""
builtinfuncs = Dict(
:(zero) => "wat",
:(zeros) => "wat",
:(randn) => "wat",
:(copy) => "wat",
:(min) => "wat",
:(max) => "wat",
:(size) => "wat",
:(getindex) => raw"""(func $getindex (param $v i32) (param $i i32) (result f32)
(f32.load (i32.add (local.get $v) (i32.shl (local.get $i) (i32.const 2)))))""",
:(setindex!) => raw"""(func $setindex! (param $v i32) (param $x f32) (param $i i32)
(f32.store (i32.add (local.get $v) (i32.shl (local.get $i) (i32.const 2))) (local.get $x)))""",
:(firstindex) => raw"""(func $firstindex (param $v i32) (result i32)
(i32.const 1))""",
:(lastindex) => raw"""(func $lastindex (param $v i32) (result i32)
(i32.trunc_f32_s (f32.load (local.get $v))))""",
:(length) => raw"""(func $length (param $v i32) (result i32)
(i32.trunc_f32_s (f32.load (local.get $v))))""",
:(iteratearray_init) => raw"""(func $iteratearray_init (param $v i32) (result f32 i32)
(f32.load (i32.add (local.get $v) (i32.const 4)))
(i32.const 1))""",
:(iteratearray) => raw"""(func $iteratearray (param $v i32) (param $i i32) (result f32 i32)
(f32.load (i32.add (local.get $v) (i32.shl (local.tee $i (i32.add (local.get $i) (i32.const 1))) (i32.const 2))))
(select (local.get $i) (i32.const 0) (i32.le_s (local.get $i) (i32.trunc_f32_s (f32.load $v)))))""",
:(iterateunitrange_init) => raw"""(func $iterateunitrange_init (param $n i32) (param $N i32) (result i32 i32)
(local.get $n)
(select (i32.const 1) (i32.const 0) (i32.le_s (local.get $n) (local.get $N))))""",
:(iterateunitrange) => raw"""(func $iterateunitrange (param $n i32) (param $N i32) (param $i i32) (result i32 i32)
(local.tee $i (i32.add (local.get $i) (i32.const 1)))
(i32.le_s (local.get $i) (local.get $N)))""",
:(iteratesteprange_init) => raw"""(func $iteratesteprange_init (param $n i32) (param $k i32) (param $N i32) (result i32 i32)
(local.get $n) (i32.const 1))""",
:(iteratesteprange) => raw"""(func $iteratesteprange (param $n i32) (param $k i32) (param $N i32) (param $i i32) (result i32 i32)
(local.tee $i (i32.add (local.get $i) (local.get $k)))
(select (i32.const 1) (i32.const 0) (i32.and (i32.ge_s (local.get $i) (local.get $n)) (i32.le_s (local.get $i) (local.get $N)))))""",
:(iteratecollection_init) => raw"""(func $iteratecollection_init (param $v i32) (result i32) (result f32)
TODO
""",
:(iteratecollection) => raw"""(func $iteratecollection (param $v i32) (param $i i32) (result i32) (result f32)
TODO
"""
)

# TODO
# https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array-1


text:
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Utility: makeminmax
import GreyDecision.Moora: moora

@testset "Moora with white numbers" begin
    tol = 0.0001
    w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
    Amat = [
        GreyNumber{Float64}(100.0, 100.0) GreyNumber{Float64}(92.0, 92.0) GreyNumber{Float64}(10.0, 10.0) GreyNumber{Float64}(2.0, 2.0) GreyNumber{Float64}(80.0, 80.0) GreyNumber{Float64}(70.0, 70.0) GreyNumber{Float64}(95.0, 95.0) GreyNumber{Float64}(80.0, 80.0)
        GreyNumber{Float64}(80.0, 80.0) GreyNumber{Float64}(70.0, 70.0) GreyNumber{Float64}(8.0, 8.0) GreyNumber{Float64}(4.0, 4.0) GreyNumber{Float64}(100.0, 100.0) GreyNumber{Float64}(80.0, 80.0) GreyNumber{Float64}(80.0, 80.0) GreyNumber{Float64}(90.0, 90.0)
        GreyNumber{Float64}(90.0, 90.0) GreyNumber{Float64}(85.0, 85.0) GreyNumber{Float64}(5.0, 5.0) GreyNumber{Float64}(0.0, 0.0) GreyNumber{Float64}(75.0, 75.0) GreyNumber{Float64}(95.0, 95.0) GreyNumber{Float64}(70.0, 70.0) GreyNumber{Float64}(70.0, 70.0)
        GreyNumber{Float64}(70.0, 70.0) GreyNumber{Float64}(88.0, 88.0) GreyNumber{Float64}(20.0, 20.0) GreyNumber{Float64}(18.0, 18.0) GreyNumber{Float64}(60.0, 60.0) GreyNumber{Float64}(90.0, 90.0) GreyNumber{Float64}(95.0, 95.0) GreyNumber{Float64}(85.0, 85.0)
    ]
    fns = makeminmax([maximum, maximum, maximum, maximum, maximum, maximum, maximum, maximum])

    result = moora(Amat, w, fns)
    scores = result[:scores]

    @test result[:bestindex] == 4

    @test isapprox(scores[1], GreyNumber(0.33159387), atol = tol)
    @test isapprox(scores[2], GreyNumber(0.29014464), atol = tol)
    @test isapprox(scores[3], GreyNumber(0.37304311), atol = tol)
    @test isapprox(scores[4], GreyNumber(0.01926526), atol = tol)
 end


text: 1 != "foo"

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia