text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediapackage
using AWS.Compat
using AWS.UUIDs

"""
    configure_logs(id)
    configure_logs(id, params::Dict{String,<:Any})

Changes the Channel's properities to configure log subscription

# Arguments
- `id`: The ID of the channel to log subscription.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"egressAccessLogs"`:
- `"ingressAccessLogs"`:
"""
configure_logs(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/configure_logs"; aws_config=aws_config)
configure_logs(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/configure_logs", params; aws_config=aws_config)

"""
    create_channel(id)
    create_channel(id, params::Dict{String,<:Any})

Creates a new Channel.

# Arguments
- `id`: The ID of the Channel. The ID must be unique within the region and it
cannot be
  changed after a Channel is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: A short text description of the Channel.
- `"tags"`:
"""
create_channel(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/channels", Dict{String, Any}("id"=>id); aws_config=aws_config)
create_channel(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/channels", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), params)); aws_config=aws_config)

"""
    create_harvest_job(end_time, id, origin_endpoint_id, s3_destination, start_time)
    create_harvest_job(end_time, id, origin_endpoint_id, s3_destination, start_time, params::Dict{String,<:Any})

Creates a new HarvestJob record.

# Arguments
- `end_time`: The end of the time-window which will be harvested
- `id`: The ID of the HarvestJob. The ID must be unique within the region
and it cannot be
  changed after the HarvestJob is submitted
- `origin_endpoint_id`: The ID of the OriginEndpoint that the HarvestJob will harvest
  from.
This cannot be changed after the HarvestJob is submitted.
- `s3_destination`:
- `start_time`: The start of the time-window which will be harvested

"""
create_harvest_job(endTime, id, originEndpointId, s3Destination, startTime; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/harvest_jobs", Dict{String, Any}("endTime"=>endTime, "id"=>id, "originEndpointId"=>originEndpointId, "s3Destination"=>s3Destination, "startTime"=>startTime); aws_config=aws_config)
create_harvest_job(endTime, id, originEndpointId, s3Destination, startTime, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/harvest_jobs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "id"=>id, "originEndpointId"=>originEndpointId, "s3Destination"=>s3Destination, "startTime"=>startTime), params)); aws_config=aws_config)

"""
    create_origin_endpoint(channel_id, id)
    create_origin_endpoint(channel_id, id, params::Dict{String,<:Any})

Creates a new OriginEndpoint record.

# Arguments
- `channel_id`: The ID of the Channel that the OriginEndpoint will be associated with.
This
  cannot be changed after the OriginEndpoint is created.
- `id`: The ID of the OriginEndpoint.  The ID must be unique within the region
and it
  cannot be changed after the OriginEndpoint is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorization"`:
- `"cmafPackage"`:
- `"dashPackage"`:
- `"description"`: A short text description of the OriginEndpoint.
- `"hlsPackage"`:
- `"manifestName"`: A short string that will be used as the filename of the OriginEndpoint
  URL (defaults to \"index\").
- `"mssPackage"`:
- `"origination"`: Control whether origination of video is allowed for this OriginEndpoint.
  If set to ALLOW, the OriginEndpoint
may by requested, pursuant to any other form of access
  control. If set to DENY, the OriginEndpoint may not be
requested. This can be helpful for
  Live to VOD harvesting, or for temporarily disabling origination
- `"startoverWindowSeconds"`: Maximum duration (seconds) of content to retain for startover
  playback.
If not specified, startover playback will be disabled for the OriginEndpoint.
- `"tags"`:
- `"timeDelaySeconds"`: Amount of delay (seconds) to enforce on the playback of live
  content.
If not specified, there will be no time delay in effect for the OriginEndpoint.
- `"whitelist"`: A list of source IP CIDR blocks that will be allowed to access the
  OriginEndpoint.
"""
create_origin_endpoint(channelId, id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/origin_endpoints", Dict{String, Any}("channelId"=>channelId, "id"=>id); aws_config=aws_config)
create_origin_endpoint(channelId, id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/origin_endpoints", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelId"=>channelId, "id"=>id), params)); aws_config=aws_config)

"""
    delete_channel(id)
    delete_channel(id, params::Dict{String,<:Any})

Deletes an existing Channel.

# Arguments
- `id`: The ID of the Channel to delete.

"""
delete_channel(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/channels/$(id)"; aws_config=aws_config)
delete_channel(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/channels/$(id)", params; aws_config=aws_config)

"""
    delete_origin_endpoint(id)
    delete_origin_endpoint(id, params::Dict{String,<:Any})

Deletes an existing OriginEndpoint.

# Arguments
- `id`: The ID of the OriginEndpoint to delete.

"""
delete_origin_endpoint(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/origin_endpoints/$(id)"; aws_config=aws_config)
delete_origin_endpoint(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/origin_endpoints/$(id)", params; aws_config=aws_config)

"""
    describe_channel(id)
    describe_channel(id, params::Dict{String,<:Any})

Gets details about a Channel.

# Arguments
- `id`: The ID of a Channel.

"""
describe_channel(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/channels/$(id)"; aws_config=aws_config)
describe_channel(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/channels/$(id)", params; aws_config=aws_config)

"""
    describe_harvest_job(id)
    describe_harvest_job(id, params::Dict{String,<:Any})

Gets details about an existing HarvestJob.

# Arguments
- `id`: The ID of the HarvestJob.

"""
describe_harvest_job(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs/$(id)"; aws_config=aws_config)
describe_harvest_job(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs/$(id)", params; aws_config=aws_config)

"""
    describe_origin_endpoint(id)
    describe_origin_endpoint(id, params::Dict{String,<:Any})

Gets details about an existing OriginEndpoint.

# Arguments
- `id`: The ID of the OriginEndpoint.

"""
describe_origin_endpoint(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints/$(id)"; aws_config=aws_config)
describe_origin_endpoint(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints/$(id)", params; aws_config=aws_config)

"""
    list_channels()
    list_channels(params::Dict{String,<:Any})

Returns a collection of Channels.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Upper bound on number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
"""
list_channels(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/channels"; aws_config=aws_config)
list_channels(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/channels", params; aws_config=aws_config)

"""
    list_harvest_jobs()
    list_harvest_jobs(params::Dict{String,<:Any})

Returns a collection of HarvestJob records.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"includeChannelId"`: When specified, the request will return only HarvestJobs associated
  with the given Channel ID.
- `"includeStatus"`: When specified, the request will return only HarvestJobs in the given
  status.
- `"maxResults"`: The upper bound on the number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
"""
list_harvest_jobs(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs"; aws_config=aws_config)
list_harvest_jobs(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/harvest_jobs", params; aws_config=aws_config)

"""
    list_origin_endpoints()
    list_origin_endpoints(params::Dict{String,<:Any})

Returns a collection of OriginEndpoint records.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"channelId"`: When specified, the request will return only OriginEndpoints associated
  with the given Channel ID.
- `"maxResults"`: The upper bound on the number of records to return.
- `"nextToken"`: A token used to resume pagination from the end of a previous request.
"""
list_origin_endpoints(; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints"; aws_config=aws_config)
list_origin_endpoints(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/origin_endpoints", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource-arn)
    list_tags_for_resource(resource-arn, params::Dict{String,<:Any})

# Arguments
- `resource-arn`:

"""
list_tags_for_resource(resource_arn; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("GET", "/tags/$(resource-arn)", params; aws_config=aws_config)

"""
    rotate_channel_credentials(id)
    rotate_channel_credentials(id, params::Dict{String,<:Any})

Changes the Channel's first IngestEndpoint's username and password. WARNING - This API is
deprecated. Please use RotateIngestEndpointCredentials instead

# Arguments
- `id`: The ID of the channel to update.

"""
rotate_channel_credentials(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/credentials"; aws_config=aws_config)
rotate_channel_credentials(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/credentials", params; aws_config=aws_config)

"""
    rotate_ingest_endpoint_credentials(id, ingest_endpoint_id)
    rotate_ingest_endpoint_credentials(id, ingest_endpoint_id, params::Dict{String,<:Any})

Rotate the IngestEndpoint's username and password, as specified by the IngestEndpoint's id.

# Arguments
- `id`: The ID of the channel the IngestEndpoint is on.
- `ingest_endpoint_id`: The id of the IngestEndpoint whose credentials should be rotated

"""
rotate_ingest_endpoint_credentials(id, ingest_endpoint_id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/ingest_endpoints/$(ingest_endpoint_id)/credentials"; aws_config=aws_config)
rotate_ingest_endpoint_credentials(id, ingest_endpoint_id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)/ingest_endpoints/$(ingest_endpoint_id)/credentials", params; aws_config=aws_config)

"""
    tag_resource(resource-arn, tags)
    tag_resource(resource-arn, tags, params::Dict{String,<:Any})

# Arguments
- `resource-arn`:
- `tags`:

"""
tag_resource(resource_arn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource-arn, tag_keys)
    untag_resource(resource-arn, tag_keys, params::Dict{String,<:Any})

# Arguments
- `resource-arn`:
- `tag_keys`: The key(s) of tag to be deleted

"""
untag_resource(resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_channel(id)
    update_channel(id, params::Dict{String,<:Any})

Updates an existing Channel.

# Arguments
- `id`: The ID of the Channel to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: A short text description of the Channel.
"""
update_channel(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)"; aws_config=aws_config)
update_channel(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/channels/$(id)", params; aws_config=aws_config)

"""
    update_origin_endpoint(id)
    update_origin_endpoint(id, params::Dict{String,<:Any})

Updates an existing OriginEndpoint.

# Arguments
- `id`: The ID of the OriginEndpoint to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorization"`:
- `"cmafPackage"`:
- `"dashPackage"`:
- `"description"`: A short text description of the OriginEndpoint.
- `"hlsPackage"`:
- `"manifestName"`: A short string that will be appended to the end of the Endpoint URL.
- `"mssPackage"`:
- `"origination"`: Control whether origination of video is allowed for this OriginEndpoint.
  If set to ALLOW, the OriginEndpoint
may by requested, pursuant to any other form of access
  control. If set to DENY, the OriginEndpoint may not be
requested. This can be helpful for
  Live to VOD harvesting, or for temporarily disabling origination
- `"startoverWindowSeconds"`: Maximum duration (in seconds) of content to retain for
  startover playback.
If not specified, startover playback will be disabled for the
  OriginEndpoint.
- `"timeDelaySeconds"`: Amount of delay (in seconds) to enforce on the playback of live
  content.
If not specified, there will be no time delay in effect for the OriginEndpoint.
- `"whitelist"`: A list of source IP CIDR blocks that will be allowed to access the
  OriginEndpoint.
"""
update_origin_endpoint(id; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/origin_endpoints/$(id)"; aws_config=aws_config)
update_origin_endpoint(id, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediapackage("PUT", "/origin_endpoints/$(id)", params; aws_config=aws_config)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --compile=min -O0 --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using Raw2CSV

if length(ARGS) !== 2
	println(
		"usage: ./solved_mps.jl <solver_name>" *
		" <path_to_folder_with_LP_and_MIP_subdirs>"
	)
	println(
		"Note: the script will consider every file inside the LP and MIP" *
		" subdirectories."
	)
	exit()
end

# In the logfile the CPLEX and Gurobi parameters have their original
# names, so to normalize both, we create extractors with a normalized
# column name that is different from the name in the file.
solver_parameters = if ARGS[1] == "Gurobi"
	[
		("solver_time_limit",
			key_equals_extractor("TimeLimit", NoDefault{Float64}())),
		# for some reason Threads is a floating point in the files
		("solver_threads",
			key_equals_extractor("Threads", NoDefault{Float64}())),
	]
elseif ARGS[1] == "CPLEX"
	[
		("solver_time_limit",
			key_equals_extractor("CPXPARAM_TimeLimit", NoDefault{Float64}())),
		# for some reason Threads is a floating point in the files
		("solver_threads",
			key_equals_extractor("CPXPARAM_Threads", NoDefault{Float64}())),
	]
else
	error("The first argument of the script must be either CPLEX or Gurobi.")
end

function not_found(s)
	return "key '$s' not found in the log"
end

csv = gather_csv_from_folders(
	[joinpath(ARGS[2], "LP"), joinpath(ARGS[2], "MIP")],
	vcat([
		key_equals("formulation", NoDefault{String}()),
		key_equals("instance_name", NoDefault{String}()),
		key_equals("rotation", NoDefault{Int}()),
		key_equals("filename", NoDefault{String}()),
		key_equals("filepath", NoDefault{String}()),
		key_equals("solver", NoDefault{String}()),
		key_equals("relax", NoDefault{Int}()),
	], solver_parameters, [
		key_equals("termination_status", not_found("termination_status")),
		key_equals("raw_status", not_found("raw_status")),
		key_equals("primal_status", not_found("primal_status")),
		key_equals("dual_status", not_found("dual_status")),
		key_equals("simplex_iterations", -1.0),
		key_equals("barrier_iterations", -1.0),
		key_equals("number_of_variables", -1),
		key_equals("number_of_constraints", -1),
		key_equals("objective_bound", NaN),
		key_equals("objective_value", NaN),
		key_equals("dual_objective_value", NaN),
		key_equals("solve_time", NaN),
		key_equals("node_count", NaN), # uses float because it may be very large
		key_equals("relative_gap", NaN),
		# Unfortunately, many MPS were solved without exception catching and
		# without writing this field, so `true` implies we are 100% sure
		# an exception happened, but false means nothing.
		key_equals("exception_detected", false),
	])
)
print(csv)

ext: jl


text:
module Sample

using GaPLAC
using GaPLAC.AbstractGPs
using GaPLAC.CairoMakie

function run(args)
    @info "running 'sample'"
    @info args
    gpspec = GaPLAC.gp_formula(args["formula"])
    @debug "GP formula" GaPLAC.formula(gpspec)

    eq, vars = GaPLAC._apply_vars(GaPLAC.formula(gpspec))
    kernels = GaPLAC._walk_kernel(eq)
    length(vars) == length(kernels) || error("Something went wrong with equation parsing, number of variables should == number of kernels")

    @debug "GP equation" eq
    @debug "Model variables" vars
    gp = GP(eq)

    @debug "GP" gp

    atdict = Dict{Symbol, Any}()
    ats = Meta.parse(args["at"])
    @debug "Passed ranges" ats

    exprs = ats.head == :toplevel ? ats.args : [ats]

    for expr in exprs
        expr.head == Symbol("=") || error("Only assignments allowed in `--at` argument")
        var = expr.args[1]
        val = GaPLAC.eval(expr.args[2])
        atdict[var] = val
    end

    for (i, var) in enumerate(vars)
        if !in(var, keys(atdict))
            atdict[var] = GaPLAC._default_range(kernels[i])
        end
    end

    @debug "Inferred ranges" atdict
    @debug "Number of combinations" (*(length.(values(atdict))...))

    df = GaPLAC._make_test_df((atdict[v] for v in vars)...; vars)
    X = RowVecs(Matrix(df))
    df[!, GaPLAC.response(gpspec)] = rand(gp(X, 0.1))

    GaPLAC._df_output(df, args)

    if !isnothing(args["plot"])
        if length(vars) > 1
            @warn """
                Auto-generated sample plots not available for multivariable models,
                use `--output` to make a table and plot manually
                """
        else
            @info "Plotting output"
            x = df[!, first(vars)]
            y = df[!, GaPLAC.response(gpspec)]
            fx = AbstractGPs.FiniteGP(gp, x, 0.1)
            pgp = posterior(fx, y)
            xmin, xmax = extrema(x) .+ (-1, 1)
            xtest = range(xmin, xmax, length = 100)
            ym, yvar = mean_and_var(pgp, xtest)

            fig, ax, l = scatter(x, y, color=:purple, label="samples")
            # srt = sortperm(x)
            # lines!(x[srt],y[srt], color=:purple)

            lines!(xtest, ym, color=:dodgerblue, label="mean posterior")
            band!(xtest, ym .- yvar, ym .+ yvar, color=(:dodgerblue, 0.3), label="var posterior")

            ax.xlabel = string(first(vars))
            ax.ylabel = string(GaPLAC.response(gpspec))
            ax.title = "Sample from posterior, x from $(round(xmin, digits=2)) to $(round(xmax, digits=2))"
            axislegend()
            save(args["plot"], fig)
        end
    end

    return df, gp
end

end # module

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using CSyntax
using StaticArrays

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))
include(joinpath(@__DIR__, "camera.jl"))

# init window
width, height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# set camera
camera = PerspectiveCamera()
setposition!(camera, [0.0, 0.0, 2.0])

# vertex and normal
points = GLfloat[ 0.0, 0.5, 0.0,
                  0.5,-0.5, 0.0,
                 -0.5,-0.5, 0.0]

normals = GLfloat[0.0, 0.0, 1.0,
                  0.0, 0.0, 1.0,
                  0.0, 0.0, 1.0]

# create VBO
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

normal_vbo = GLuint(0)
@c glGenBuffers(1, &normal_vbo)
glBindBuffer(GL_ARRAY_BUFFER, normal_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(normals), normals, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, normal_vbo)
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "phong.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "phong.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)

model_loc = glGetUniformLocation(shader_prog, "model_mat")
view_loc = glGetUniformLocation(shader_prog, "view_mat")
proj_loc = glGetUniformLocation(shader_prog, "projection_mat")
glUseProgram(shader_prog)
model_mat = Matrix{GLfloat}(I, 4, 4)
glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_mat)
glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))

# enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

let
updatefps = FPSCounter()
# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindBuffer(GL_ARRAY_BUFFER, vao)
    model_mat[1,4] = sin(time())
    glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_mat)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # move camera
    updatecamera!(window, camera)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_path: OpenGL 4 Tutorials/08_phong/phong.jl

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_path: OpenGL 4 Tutorials/08_phong/phong.jl

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_path: OpenGL 4 Tutorials/08_phong/phong.jl

max_forks_repo_name: Gnimuc/Videre


text:
module IncomeTaxCalculations

#
# This module calculates income tax for an individual (and also, at the end a married couple).
# It's based mainly on the descriptions and examples in the excellent Melville's Taxation:
# Melville, Alan. 2019. Melville’s Taxation: Finance Act 2019. 25th ed. London. (and 2020 edition).
#

using Dates
using Dates: Date, now, TimeType, Year
using Parameters: @with_kw

using ScottishTaxBenefitModel

using .Definitions

using .ModelHousehold:
    Person

using .STBParameters:
    IncomeTaxSys

using .GeneralTaxComponents:
    RateBands,
    TaxResult,
    calctaxdue,
    delete_thresholds_up_to
    # , *

using .Utils:
    get_if_set

using .STBIncomes

using .Results:
    BenefitUnitResult,
    IndividualResult,
    ITResult

export
    apply_allowance,
    calc_income_tax!,
    calculate_company_car_charge,
    old_enough_for_mca

## FIXME just use the dict..
function guess_car_percentage_2020_21( sys :: IncomeTaxSys, company_car_fuel_type :: Fuel_Type )
    return sys.company_car_charge_by_CO2_emissions[company_car_fuel_type]
end

function calculate_company_car_charge(
    pers   :: Person,
    sys    :: IncomeTaxSys,
    calculator :: Function = guess_car_percentage_2020_21 ) :: Real
    value = max(0.0, pers.company_car_value-
        pers.company_car_contribution )
    if pers.fuel_supplied > 0.0
        value += sys.fuel_imputation
    end
    prop = calculator( sys, pers.company_car_fuel_type )
    value * prop
end

"""
Very rough approximation to MCA age - ignores all months since we don't have that in a typical dataset
TODO maybe overload this with age as a Date?
"""
function old_enough_for_mca(
    sys            :: IncomeTaxSys,
    age            :: Integer,
    model_run_date :: TimeType = now() ) :: Bool
    (model_run_date - Year(age)) < sys.mca_date
end

function calculate_allowance( pers::Person, sys :: IncomeTaxSys ) :: Real
    allowance = sys.personal_allowance
    if pers.registered_blind
        allowance += sys.blind_persons_allowance
    end
    allowance
end

function apply_allowance( allowance::Real, income::Real )::Tuple where RT<:Real
    r = max( 0.0, income - allowance )
    allowance = max(0.0, allowance-income)
    allowance,r
end

"""
  from Melville, ch13.
  Changes are in itres: add to pension fields and extend bands
  notes:
  Melville talks of "earned income" - using non-savings income
  FIXME: check: does the personal_allowance_income_limit get bumped up?
"""
function calculate_pension_taxation!(
    pres   :: IndividualResult,
    sys    :: IncomeTaxSys,
    pers   ::Person,
    total_income::Real,
    earned_income:: Real )

    pres.it.savings_thresholds = copy( sys.savings_thresholds )
    pres.it.dividend_thresholds = copy( sys.dividend_thresholds )
    pres.it.non_savings_thresholds = copy( sys.non_savings_thresholds )
    # fixme check avs here
    # FIXME not right?
    avc = pres.income[AVCS]
    pen = pres.income[PENSION_CONTRIBUTIONS_EMPLOYEE] +
        pres.income[PENSION_CONTRIBUTIONS_EMPLOYER]
    eligible_contribs = avc + pen
    if eligible_contribs <= 0.0
        return
    end

    max_relief = sys.pension_contrib_annual_allowance
    if total_income < sys.pension_contrib_basic_amount
        max_relief = sys.pension_contrib_basic_amount
    end
    if total_income > sys.pension_contrib_threshold_income
        excess = total_income - sys.pension_contrib_threshold_income
        max_relief = max( sys.pension_contrib_annual_minimum,
            sys.pension_contrib_annual_allowance - excess*sys.pension_contrib_withdrawal_rate )
    end
    # println( "total_income $total_income max_relief $max_relief eligible_contribs $eligible_contribs" )
    eligible_contribs = min( eligible_contribs, max_relief );
    pres.it.pension_eligible_for_relief = eligible_contribs
    basic_rate = sys.non_savings_rates[ sys.non_savings_basic_rate ]
    # println("total_income=$total_income max_relief=$max_relief eligible_contribs=$eligible_contribs basic_rate=$basic_rate sys.pension_contrib_withdrawal_rate=$(sys.pension_contrib_withdrawal_rate)")
    gross_contribs = eligible_contribs/(1-basic_rate)
    pres.it.pension_relief_at_source = gross_contribs - eligible_contribs
    pres.it.non_savings_thresholds .+= gross_contribs
    pres.it.savings_thresholds .+= gross_contribs
    pres.it.dividend_thresholds .+= gross_contribs
end

"""

Complete(??) income tax calculation, based on the Scottish/UK 2019 system.
Mostly taken from Melville (2019) chs 2-4.

FIXME this is too long and needs broken up.

problems:

1. we do this in strict non-savings, savings, dividends order; see 2(11) for examples where it's now advantageous to use a different order
2.

"""
function calc_income_tax!(
    pres   :: IndividualResult,
    pers   :: Person,
    sys    :: IncomeTaxSys,
    spouse_transfer :: Real = 0.0 )
    if sys.abolished
        return;
    end
    total_income = isum( pres.income, sys.all_taxable )
    non_savings_income = isum( pres.income, sys.non_savings_income )

    savings_income = isum( pres.income, sys.savings_income )
    dividends_income = isum( pres.income, sys.dividend_income )
    allowance = calculate_allowance( pers, sys )
    # allowance reductions goes here

    adjusted_net_income = total_income

    calculate_pension_taxation!( pres, sys, pers, total_income, non_savings_income )

    # adjusted_net_income -= pres.it.pension_eligible_for_relief

    adjusted_net_income += calculate_company_car_charge(pers, sys)
    # ...

    non_savings_tax = TaxResult(0.0, 0)
    savings_tax = TaxResult(0.0, 0)
    dividend_tax = TaxResult(0.0, 0)

    if adjusted_net_income > sys.personal_allowance_income_limit
        allowance =
            max(0.0,
                allowance -
                    sys.personal_allowance_withdrawal_rate*(
                        adjusted_net_income - sys.personal_allowance_income_limit ))
    end
    taxable_income = adjusted_net_income-allowance
    # note: we copy from the expanded versions from pension_contributions
    savings_thresholds = deepcopy( pres.it.savings_thresholds )
    savings_rates = deepcopy( sys.savings_rates )
    # FIXME model all this with parameters
    toprate = size( savings_thresholds )[1]
    non_savings_taxable = 0.0
    savings_taxable = 0.0
    dividends_taxable = 0.0

    if taxable_income > 0
        allowance,non_savings_taxable = apply_allowance( allowance, non_savings_income )
        non_savings_tax = calctaxdue(
            taxable=non_savings_taxable,
            rates=sys.non_savings_rates,
            thresholds=pres.it.non_savings_thresholds )

        # horrific savings calculation see Melville Ch2 "Savings Income" & examples 2-3
        # FIXME Move to separate function
        # delete the starting bands up to non_savings taxabke icome
        savings_rates, savings_thresholds = delete_thresholds_up_to(
            rates=savings_rates,
            thresholds=savings_thresholds,
            upto=non_savings_taxable );
        if sys.personal_savings_allowance > 0
            psa = sys.personal_savings_allowance
            # println( "taxable income $taxable_income sys.savings_thresholds[2] $(sys.savings_thresholds[2])")
            if taxable_income > sys.savings_thresholds[toprate]
                psa = 0.0
            elseif taxable_income > sys.savings_thresholds[2] # above the basic rate FIXME parameterise '2' here
                psa *= 0.5 # FIXME parameterise 0.5 here
            end
            if psa > 0.0 ## if we haven't deleted the zero band already, just widen it
                if savings_rates[1] == 0.0
                    savings_thresholds[1] += psa;
                else ## otherwise, insert a  new one.
                    savings_thresholds = vcat([psa], savings_thresholds )
                    savings_rates = vcat([0.0], savings_rates )
                end
            end
            pres.it.personal_savings_allowance = psa
        end # we have a personal_savings_allowance
        pres.it.savings_rates = savings_rates
        pres.it.savings_thresholds= savings_thresholds
        allowance,savings_taxable = apply_allowance( allowance, savings_income )
        savings_tax = calctaxdue(
            taxable=savings_taxable,
            rates=savings_rates,
            thresholds=savings_thresholds )

        # Dividends
        # see around example 8-9 ch2
        allowance,dividends_taxable =
            apply_allowance( allowance, dividends_income )
        dividend_rates=deepcopy(sys.dividend_rates)
        dividend_thresholds=deepcopy(pres.it.dividend_thresholds )
        # always preserve any bottom zero rate
        add_back_zero_band = false
        zero_band = 0.0
        used_thresholds = non_savings_taxable+savings_taxable
        copy_start = 1
        # handle the zero rate
        if dividend_rates[1] == 0.0
            add_back_zero_band = true
            zero_band = dividend_thresholds[1]
            used_thresholds += min( zero_band, dividends_taxable )
            copy_start = 2
        end
        dividend_rates, dividend_thresholds =
            delete_thresholds_up_to(
                rates=dividend_rates[copy_start:end],
                thresholds=dividend_thresholds[copy_start:end],
                upto=used_thresholds );
        if add_back_zero_band
            dividend_rates = vcat( [0.0], dividend_rates )
            dividend_thresholds .+= zero_band # push all up
            dividend_thresholds = vcat( zero_band, dividend_thresholds )
        end
        pres.it.dividend_rates = dividend_rates
        pres.it.dividend_thresholds= dividend_thresholds

        dividend_tax = calctaxdue(
            taxable=dividends_taxable,
            rates=dividend_rates,
            thresholds=dividend_thresholds )
    else # some allowance left
        allowance = -taxable_income # e.g. allowance - taxable_income
    end

    #
    # tax reducers
    #
    @assert non_savings_tax.due >= 0
    @assert savings_tax.due >= 0
    @assert dividend_tax.due >= 0
    total_tax = non_savings_tax.due+savings_tax.due+dividend_tax.due
    nst = non_savings_tax.due
    st = savings_tax.due
    dt = dividend_tax.due
    if spouse_transfer > 0
        # note that if scotland sets a high basic rate the sp_reduction here
        # can wipe out all savings taxes which should be remitted to uk gov
        sp_reduction =
            sys.non_savings_rates[sys.non_savings_basic_rate]*spouse_transfer
        total_tax = max( 0.0, total_tax - sp_reduction )
        # assign spouse reduction to non-savings - needed for Scottish income tax
        nst -= sp_reduction
        # println( "nst $nst sp_reduction=$sp_reduction total_tax=$total_tax")
        if nst < 0.0
            # assign any leftover to savings tax
            st += nst # subtract, really
            nst = 0.0
        end
        if st < 0.0
            # assign any leftover to savings tax
            dt = max(0.0, dt+st) # subtract, really
            st = 0.0
        end
    end

    # total_tax -= pres.it.pension_relief_at_source
    pres.income[INCOME_TAX] = total_tax

    pres.it.taxable_income = taxable_income
    pres.it.allowance = allowance
    pres.it.total_income = total_income
    pres.it.adjusted_net_income = adjusted_net_income

    # memo items, but we need nst for Scottish Income Tax
    pres.it.non_savings_tax = nst # non_savings_tax.due, less transferrable allowances
    pres.it.non_savings_income = non_savings_income
    pres.it.non_savings_band = non_savings_tax.end_band
    pres.it.non_savings_taxable = non_savings_taxable

    pres.it.savings_tax = st
    pres.it.savings_band = savings_tax.end_band
    pres.it.savings_income = savings_income
    pres.it.savings_taxable = savings_taxable

    pres.it.dividends_tax = dt
    pres.it.dividend_band = dividend_tax.end_band
    pres.it.dividends_income = dividends_income
    pres.it.dividends_taxable = dividends_taxable

    pres.it.unused_allowance = allowance
end

function allowed_to_transfer_allowance(
    sys  :: IncomeTaxSys;
    from :: ITResult,
    to   :: ITResult ) :: Bool

    can_transfer :: Bool = true
    if ! (from.unused_allowance > 0.0 &&
        to.unused_allowance <= 0.0)
         # nothing to transfer - this is actually wrong since
         # you can opt to transfer some allowance even if you
         # can technically use it.
        can_transfer = false
    elseif to.savings_band > sys.savings_basic_rate ||
        to.non_savings_band > sys.non_savings_basic_rate ||
        to.dividend_band > sys.dividend_basic_rate
        can_transfer = false
    end
    ## TODO disallow if mca claimed
    return can_transfer
end # can_transfer

function calculate_mca(
    pers :: Person,
    tax :: ITResult,
    sys :: IncomeTaxSys)::Real
    ## FIXME parameterise this
    mca = sys.married_couples_allowance
    if tax.adjusted_net_income > sys.mca_income_maximum
        mca = max( sys.mca_minimum, mca -
           (tax.adjusted_net_income-sys.mca_income_maximum)*sys.mca_withdrawal_rate)
    end
    mca * sys.mca_credit_rate
end

"""
 FIXME maybe send an actual benefit unit in here? With kids..
"""
function calc_income_tax!(
    bres   :: BenefitUnitResult,
    head   :: Person,
    spouse :: Union{Nothing, Person},
    sys    :: IncomeTaxSys )
    hdres = bres.pers[head.pid]
    calc_income_tax!( hdres, head, sys )
    # FIXME the transferable stuff here
    # is not quite right as you can elect to transfer more than
    # the surplus allowance in some cases.
    # also - add in restrictions on transferring to
    # higher rate payers.
    if spouse !== nothing
        spres = bres.pers[spouse.pid]
        calc_income_tax!( spres, spouse, sys )
        # This is not quite right - you can't claim the
        # MCA AND transfer an allowance. We're assuming
        # always MCA first (I think it's always more valuable?)
        if old_enough_for_mca( sys, head.age ) || old_enough_for_mca( sys, spouse.age )
            # shoud usually just go to the head but.. some stuff about partner
            # with greater income if married after 2005 and you can elect to do this if
            # married before, so:
            if hdres.it.adjusted_net_income > spres.it.adjusted_net_income
                hdres.it.mca = calculate_mca( head, hdres.it, sys )
                hdres.income[INCOME_TAX] = max( 0.0, hdres.income[INCOME_TAX]- hdres.it.mca )
                # FIXME not quite right - trickle down to savings/divs
                if hdres.it.mca > 0
                    hdres.it.non_savings_tax = max( 0.0, hdres.it.non_savings_tax - hdres.it.mca )
                end
            else
                spres.it.mca = calculate_mca( spouse, spres.it, sys )
                spres.income[INCOME_TAX] = max( 0.0, spres.income[INCOME_TAX] - spres.it.mca )
                if spres.it.mca > 0
                    spres.it.non_savings_tax = max( 0.0, spres.it.non_savings_tax - spres.it.mca )
                end
            end
        end
        if spres.it.mca == 0.0 == hdres.it.mca
            if allowed_to_transfer_allowance( sys, from=spres.it, to=hdres.it )
                transferable_allow = min( spres.it.unused_allowance, sys.marriage_allowance )
                calc_income_tax!( hdres, head, sys, transferable_allow )
                hdres.it.transferred_allowance = transferable_allow
            elseif allowed_to_transfer_allowance( sys, from=hdres.it, to=spres.it )
                transferable_allow = min( hdres.it.unused_allowance, sys.marriage_allowance )
                calc_income_tax!( spres, spouse, sys, transferable_allow )
                spres.it.transferred_allowance = transferable_allow
            end
        end
        if spres.income[INCOME_TAX] < 0.0
            println( spouse.income )
            println( spres.it )
            println( spres.income )
            @assert false "income tax can't be negative; spouse; pid=$(spouse.pid)"
        end
        if spres.it.non_savings_tax > spres.income[INCOME_TAX]
            println( spouse.income )
            println( spres.it )
            println( spres.income )
            @assert false "non savings part can't be greater than total income tax; spouse ; pid=$(spouse.pid)"
        end
    end
    if hdres.income[INCOME_TAX] < 0.0
        println( head.income )
        println( hdres.it )
        println( hdres.income )

        @assert false "income tax can't be negative; head; pid=$(head.pid)"
    end
    if hdres.it.non_savings_tax > hdres.income[INCOME_TAX]
        println( head.income )
        println( hdres.it )
        println( hdres.income )
        @assert false "non savings part can't be greater than total income tax; head ; pid=($head.pid)"
    end
end # calc_income_tax

end # module

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using GeometricFlux
using GeometricFlux: neighbors, laplacian_matrix
using LightGraphs
using LightGraphs: neighbors, laplacian_matrix
using SimpleWeightedGraphs
using Test

tests = [
    "layers/msgpass",
    "layers/conv",
    "layers/pool",
    "models",
    "linalg",
    "scatter",
    "graph/simplegraphs",
    "graph/weightedgraphs",
    "graph/utils",
]

@testset "GeometricFlux" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
"Convert PDDL name-with-hyphens to CamelCase type name."
function pddl_to_type_name(name)
    words = split(lowercase(string(name)), '-', keepempty=false)
    return join(uppercasefirst.(words))
end

"Returns a list of pairs from a collection, sorted by key."
sortedpairs(collection) = sort(collect(pairs(collection)), by=first)

"Returns the sorted list of keys for a collection."
sortedkeys(collection) = sort(collect(keys(collection)))


text:
#-- EventTime ---------------

## Type constructors

"""
    EventTime{T}

Immutable object containing the real-valued time to an event as well as an indicator of
whether the time corresponds to an observed event (`true`) or right censoring (`false`).
"""
struct EventTime{T<:Real}
    time::T
    status::Bool
end

EventTime(time::T) where {T<:Real} = EventTime{T}(time, true)

## Overloaded Base functions

Base.eltype(::EventTime{T}) where {T} = T
Base.show(io::IO, ev::EventTime) = print(io, ev.time, ifelse(ev.status, "", "+"))

Base.convert(T::Type{<:Real}, ev::EventTime) = convert(T, ev.time)
Base.convert(T::Type{EventTime}, x::Real) = EventTime(x)

function Base.isless(t1::EventTime, t2::EventTime)
    if t1.time == t2.time
        # When two EventTimes have the same observed time, we compare the event
        # status. Observed events compare less than censored events, since if the
        # censored event were to occur then it would happen at or after the given
        # time (by definition).
        return isevent(t1) && iscensored(t2)
    else
        return isless(t1.time, t2.time)
    end
end

## New functions

isevent(ev::EventTime) = ev.status
iscensored(ev::EventTime) = !ev.status


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter.

# Fields
$(TYPEDFIELDS)
"""
struct CovarianceMatrix <: AbstractConstraint
    function CovarianceMatrix()
        return new()
    end
end

############################################################################################
function construct_transform(constraint::CovarianceMatrix, val)
    transform = Bijectors.PDBijector()
    return transform, Bijectors.inverse(transform)
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::CovarianceMatrix,
    val::Matrix{R},
) where {R<:Real}
    ArgCheck.@argcheck LinearAlgebra.issymmetric(val)
    return true
end

idx_upper = tag(zeros(3,3), false, true)
sum(idx_upper)
############################################################################################
#!TODO: Works with flatten/unflatten - but constraint/unconstraint seems to deduce wrong type for ReverseDiff from Bijector - works fine with ForwardDiff/Zygote
#!NOTE: Bijectors map to lower triangular matrix while most AD libraries evaluate upper triangular matrices.
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenStrict,
    constraint::C,
    x::Matrix{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{CovarianceMatrix, Distributions.InverseWishart, Bijectors.PDBijector}
}
    #!NOTE: PDBijector seems to unconstrain to a Lower Diagonal Matrix
    idx_upper = tag(x, false, true)
    len = length(x)
    len_unflat = sum(idx_upper)
    function flatten_CovMatrix(x::AbstractMatrix{R}) where {R<:Real}
        ArgCheck.@argcheck length(x) == len
        return Vector{T}(flatten_Symmetric(x, idx_upper))
    end
    buffer_unflat = zeros(R, size(x))
    function CovMatrix_from_vec(v::Union{<:Real,AbstractVector{<:Real}})
        ArgCheck.@argcheck length(v) == len_unflat
        return Symmetric_from_flatten!(buffer_unflat, v, idx_upper)
    end
    return flatten_CovMatrix, CovMatrix_from_vec
end
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenFlexible,
    constraint::C,
    x::Matrix{R},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    R<:Real,
    C<:Union{CovarianceMatrix, Distributions.InverseWishart, Bijectors.PDBijector}
}
    #!NOTE: PDBijector seems to unconstrain to a Lower Diagonal Matrix
    idx_upper = tag(x, false, true)
    len = length(x)
    len_unflat = sum(idx_upper)
    function flatten_CovMatrix_AD(x::AbstractMatrix{R}) where {R<:Real}
        ArgCheck.@argcheck length(x) == len
        return Vector{R}(flatten_Symmetric(x, idx_upper))
    end
    dims = size(x)
    function CovMatrix_from_vecAD(v::Union{<:Real,AbstractVector{<:Real}})
        ArgCheck.@argcheck length(v) == len_unflat
        buffer_unflat = zeros(eltype(v), dims)
        return Symmetric_from_flatten!(buffer_unflat, v, idx_upper)
    end
    return flatten_CovMatrix_AD, CovMatrix_from_vecAD
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    CovarianceMatrix,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
#try
using Plots
#catch
#   import Pkg; Pkg.add("Plots")
#   Pkg.add("UnicodePlots")
#   using Plots
#end

# FIXME: using StatsBase
#using Plots
#using UnicodePlots
#unicodeplots()
function raster(p)
    fire = p.records[:fire]
    x, y = SNNFloat[], SNNFloat[]
    for t = eachindex(fire)
        for n in findall(fire[t])
            push!(x, t)
            push!(y, n)
        end
    end
    x, y
end

function raster(P::Array)
    y0 = SNNInt[0]
    X = SNNFloat[]; Y = SNNFloat[]
    for p in P
        x, y = raster(p)
        append!(X, x)
        append!(Y, y .+ sum(y0))
        push!(y0, p.N)
    end
    plt = scatter(X, Y, m = (1, :black), leg = :none,
                  xaxis=("t", (0, Inf)), yaxis = ("neuron",))
    y0 = y0[2:end-1]
    !isempty(y0) && hline!(plt, cumsum(y0), linecolor = :red)
    return plt
end

function vecplot(p, sym)
    v = getrecord(p, sym)
    y = hcat(v...)'
    x = 1:length(v)
    plot(x, y, leg = :none,
    xaxis=("t", extrema(x)),
    yaxis=(string(sym), extrema(y)))
end

function vecplot(P::Array, sym)
    plts = [vecplot(p, sym) for p in P]
    N = length(plts)
    plot(plts..., size = (600, 400N), layout = (N, 1))
end

function windowsize(p)
    A = sum.(p.records[:fire]) / length(p.N)
    W = round(SNNInt, 0.5p.N / mean(A)) # filter window, unit=1
end

function density(p, sym)
    X = getrecord(p, sym)
    t = 1:length(X)
    xmin, xmax = extrema(vcat(X...))
    edge = linspace(xmin, xmax, 50)
    c = center(edge)
    ρ = [fit(Histogram, x, edge).weights |> float for x in X] |> x -> hcat(x...)
    ρ = smooth(ρ, windowsize(p), 2)
    ρ ./= sum(ρ, 1)
    p = @gif for t = 1:length(X)
        bar(c, ρ[:, t], leg = false, xlabel = string(sym), yaxis = ("p", extrema(ρ)))
    end
    is_windows() && run(`powershell start $(p.filename)`)
    is_unix() && run(`xdg-open $(p.filename)`)
    p
end

function rateplot(p, sym)
    r = getrecord(p, sym)
    R = hcat(r...)
end

function rateplot(P::Array, sym)
    R = vcat([rateplot(p, sym) for p in P]...)
    y0 = [p.N for p in P][2:end-1]
    plt = heatmap(flipdim(R, 1), leg = :none)
    !isempty(y0) && hline!(plt, cumsum(y0), line = (:black, 1))
    plt
end

function activity(p)
    A = sum.(p.records[:fire]) / length(p.N)
    W = windowsize(p)
    A = smooth(A, W)
end

function activity(P::Array)
    A = activity.(P)
    t = 1:length(P[1].records[:fire])
    plot(t, A, leg=:none, xaxis=("t",), yaxis=("A", (0, Inf)))
end

function if_curve(model, current; neuron = 1, dt = 0.1ms, duration = 1second)
    E = model(neuron)
    monitor(E, [:fire])
    f = SNNFloat[]
    for I = current
        clear_records(E)
        E.I = [I]
        SNN.sim!([E], []; dt = dt, duration = duration)
        push!(f, activity(E))
    end
    plot(current, f)
end

# export density
# function density(p, sym)
#   X = getrecord(p, sym)
#   t = dt*(1:length(X))
#   xmin, xmax = extrema(vcat(X...))
#   edge = linspace(xmin, xmax, 100)
#   c = center(edge)
#   ρ = [fit(Histogram, x, edge).weights |> reverse |> float for x in X] |> x->hcat(x...)
#   ρ = smooth(ρ, windowsize(p), 2)
#   ρ ./= sum(ρ, 1)
#   surface(t, c, ρ, ylabel="p")
# end
# function density(P::Array, sym)
#   plts = [density(p, sym) for p in P]
#   plot(plts..., layout=(length(plts),1))
# end

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
include("../src/EDKit.jl")
using .EDKit
using LinearAlgebra


text:
module PipeLineTests
using Test

@testset "Builder.lt_page" begin
    using Documenter.Builder: lt_page
    # Checks to make sure that only exactly one of a<b, a==b and a>b is true for given a & b
    iscorrectisless(a,b) = sum([lt_page(a, b), a == b, lt_page(b, a)]) == 1
    # Test equal strings:
    for a in ["index.md", "foo/index.md", "foo.md", "bar/foo.md", "foo/bar/baz/qux", "", "α", "α/β", "α/index.md"]
        @test !lt_page(a, a)
        @test iscorrectisless(a, a)
    end
    # Test less thans
    for (a, b) in [
            ("a", "b"), ("α", "β"), ("b", "α"),
            # index.md takes precedence
            ("index.md", "a"),
            ("index.md", "index.mm"),
            ("index.md", "foo.md"),
            ("index.md", "bar/foo.md"),
            # Also in subdirectories:
            ("foo/index.md", "foo/a"),
            ("foo/index.md", "foo/index.mm"),
            ("foo/index.md", "foo/foo.md"),
            ("foo/index.md", "foo/bar/foo.md"),
            # But not over stuff that is outside of the subdirectory
            ("a", "foo/index.md"),
            ("α", "α/index.md"),
            ("foo/index.md", "g"),
            ("bar/index.md", "foo/index.md"),
            ("bar/qux/index.md", "foo/index.md"),
        ]
        @test lt_page(a, b)
        @test iscorrectisless(a, b)
    end

    @test sort(["foo", "bar"], lt=lt_page) == ["bar", "foo"]
    @test sort(["foo", "foo/bar"], lt=lt_page) == ["foo", "foo/bar"]
    @test sort(["foo", "f/bar"], lt=lt_page) == ["f/bar", "foo"]
    @test sort(["foo", "index.md"], lt=lt_page) == ["index.md", "foo"]
    @test sort(["foo.md", "foo/index.md", "index.md", "foo/foo.md"], lt=lt_page) == ["index.md", "foo.md", "foo/index.md", "foo/foo.md"]
    @test sort(["foo.md", "ϕωω/index.md", "index.md", "foo/foo.md"], lt=lt_page) == ["index.md", "foo.md", "foo/foo.md", "ϕωω/index.md"]
end

# Docstring signature syntax highlighting tests.
module HighlightSig
    using Test
    import Markdown
    import Documenter.Expanders: highlightsig!

    @testset "highlightsig!" begin
        s = """
                foo(bar::Baz)
            ---
                foo(bar::Baz)
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test isempty(original.content[1].language)
        @test md.content[1].language == "julia"
        @test original.content[end].language == md.content[end].language

        s = """
            ```lang
             foo(bar::Baz)
            ```
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test original == md
    end
end

end


text:
using SumOfSquares
using DynamicPolynomials

function sos_polynomial(n)
    @polyvar x
    monos = monomials(x, 0:2n)
    set = SumOfSquares.SOSPolynomialSet(
        SumOfSquares.FullSpace(), monos,
        SumOfSquares.Certificate.Remainder(
            SumOfSquares.SOSCone(),
            SumOfSquares.MonomialBasis,
            tuple()
        )
    )
    model = MOIU.UniversalFallback(MOIU.Model{Float64}())
    x = MOI.add_variables(model, n)
    scalar = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0)
    func = MOI.Utilities.vectorize([scalar for i in eachindex(monos)])
    BT = MOI.Bridges.Constraint.concrete_bridge_type(
        SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64},
        typeof(func),
        typeof(set))
    @time MOI.Bridges.Constraint.bridge_constraint(BT, model, func, set)
    return
end
#sos_polynomial(200)
sos_polynomial(2)


text:
export pairwise_marginal, pairwise_MI

using LinearAlgebra: diagind, diag
using CUDA: CUDA, CuMatrix, CuVector, CuArray

xlogx(x) =
    iszero(x) ? zero(x) : x * log(x)

xlogy(x, y) =
    iszero(x) && !isnan(y) ? zero(x) : x * log(y)

#############################
# Mutual Information on Binary Data
#############################

function pairwise_marginal(data::Union{Matrix{Bool}, CuMatrix{Bool}, BitMatrix};
        weights::Union{Vector, CuVector, Nothing}=nothing,
        pseudocount=0.0,
        Float=Float64)
    N = isnothing(weights) ? size(data, 1) : sum(weights)
    D = size(data, 2)
    base = N + pseudocount
    if data isa CuMatrix{Bool} && !isnothing(weights)
        weights = CuVector(weights)
    end

    not_data = Float.(.!data)
    data = Float.(data)
    pxy = similar(data, D, D, 4)

    if isnothing(weights)
        pxy[:, :, 4] = transpose(data) * data
        pxy[:, :, 3] = transpose(data) * not_data
        pxy[:, :, 2] = transpose(not_data) * data
    else
        pxy[:, :, 4] = transpose(data) * (data .* weights)
        pxy[:, :, 3] = transpose(data) * (not_data .* weights)
        pxy[:, :, 2] = transpose(not_data) * (data .* weights)
    end
    pxy[:, :, 1] = N .- pxy[:, :, 4] .- pxy[:, :, 3] .- pxy[:, :, 2]

    joint_count = pseudocount / 4
    pxy .+= joint_count

    # diagonal is the marginal of p(x)
    diag = diagind(pxy[:, :, 1])
    @view(pxy[:,:, 1])[diag] .+= joint_count
    @view(pxy[:,:, 2])[diag] .-= joint_count
    @view(pxy[:,:, 3])[diag] .-= joint_count
    @view(pxy[:,:, 4])[diag] .+= joint_count

    if data isa CuMatrix
        CUDA.unsafe_free!(data)
        CUDA.unsafe_free!(not_data)
    end
    pxy ./= base
end

function pairwise_MI(data::Union{Matrix{Bool}, CuMatrix{Bool}, BitMatrix};
        weights=nothing,
        pseudocount=0.0,
        Float=Float64)
    pxy = pairwise_marginal(data; weights, pseudocount, Float)
    p0 = diag(pxy[:, :, 1])
    p1 = diag(pxy[:, :, 4])
    D = size(data, 2)
    pxpy = Float.(similar(data, D, D, 4))
    pxpy[:,:,1] = p0 * p0'
    pxpy[:,:,2] = p0 * p1'
    pxpy[:,:,3] = p1 * p0'
    pxpy[:,:,4] = p1 * p1'
    pxy_log_pxy = @. xlogx(pxy)
    pxy_log_pxpy = @. xlogy(pxy, pxpy)
    MI = dropdims(sum(pxy_log_pxy - pxy_log_pxpy, dims=3), dims=3)
    if data isa CuMatrix{Bool}
        CUDA.unsafe_free!(pxy)
        CUDA.unsafe_free!(pxpy)
        CUDA.unsafe_free!(pxy_log_pxy)
        CUDA.unsafe_free!(pxy_log_pxpy)
    end
    MI
end

#############################
# Mutual Information on Categorical Data
#############################

function pairwise_marginal(data::Matrix;
        weights::Union{Vector, Nothing}=nothing,
        num_cats=maximum(data),
        pseudocount=0.0,
        Float=Float64)
    @assert minimum(data) > 0 "Categorical data labels are assumed to be indexed starting 1"

    num_samples = size(data, 1)
    num_vars = size(data, 2)

    if weights === nothing
        weights = ones(Int32, num_samples)
    else
        pseudocount = pseudocount * sum(weights) / num_samples
    end

    num_vars = size(data, 2)
    pair_margs = Array{Float}(undef, num_vars, num_vars, num_cats, num_cats)
    Z = sum(weights) + pseudocount
    single_smooth = pseudocount / num_cats
    pair_smooth = single_smooth / num_cats

    for i = 1:num_vars
        # note: multithreading needs to be on inner loop for thread-safe copying across diagonal
        Threads.@threads for j = 1:num_vars
            if i<=j
                if i!=j
                    @inbounds pair_margs[i,j,:,:] .= pair_smooth
                else
                    @inbounds pair_margs[i,j,:,:] .= zero(eltype(pair_margs))
                    for l = 1:num_cats
                        @inbounds pair_margs[i,j,l,l] = single_smooth
                    end
                end
                @simd for k = 1:size(data,1) # @avx here gives incorrect results
                    @inbounds pair_margs[i,j,data[k,i],data[k,j]] += weights[k]
                end
                @inbounds pair_margs[i,j,:,:] ./= Z
            else
                @inbounds pair_margs[i,j,:,:] .= (@view pair_margs[j,i,:,:])'
            end
            nothing
        end
    end

    return pair_margs
end

function pairwise_marginal(data::CuMatrix;
        weights::Union{CuVector, Nothing}=nothing,
        num_cats=maximum(data),
        pseudocount=0.0,
        Float=Float64)

    @assert minimum(data) > 0 "Categorical data labels are assumed to be indexed starting 1"

    num_samples = size(data, 1)
    num_vars = size(data, 2)

    if weights === nothing
        weights = ones(Int32, num_samples)
    else
        pseudocount = pseudocount * sum(weights) / num_samples
    end

    num_vars = size(data,2)
    pair_margs = CuArray{Float}(undef, num_vars, num_vars, num_cats, num_cats)
    Z = sum(weights) + pseudocount
    single_smooth = pseudocount / num_cats
    pair_smooth = single_smooth / num_cats

    data_device = CUDA.cudaconvert(data)
    weights_device = CUDA.cudaconvert(weights)
    pair_margs_device = CUDA.cudaconvert(pair_margs)

    var_indices = CuArray(1:num_vars)
    CUDA.@sync begin
        broadcast(var_indices, var_indices') do i,j
            if i <= j
                if i!=j
                    @inbounds pair_margs_device[i,j,:,:] .= pair_smooth
                else
                    @inbounds pair_margs_device[i,j,:,:] .= zero(Float)
                    for l = 1:num_cats
                        @inbounds pair_margs_device[i,j,l,l] = single_smooth
                    end
                end
                for k = 1:size(data_device,1)
                    pair_margs_device[i,j,data_device[k,i],data_device[k,j]] += weights_device[k]
                end
            end
            nothing
        end
        pair_margs ./= Z
        broadcast(var_indices, var_indices') do i,j
            if i > j
                for l = 1:num_cats, m = 1:num_cats
                    @inbounds pair_margs_device[i,j,l,m] = pair_margs_device[j,i,m,l]
                end
            end
            nothing
        end
    end
    return pair_margs
end

function pairwise_MI(data::Matrix;
            weights::Union{Vector, Nothing} = nothing,
            num_cats=maximum(data),
            pseudocount=0.0,
            Float=Float64)
    num_samples = size(data, 1)
    num_vars = size(data, 2)

    if weights === nothing
        weights = ones(Int32, num_samples)
    else
        pseudocount = pseudocount * sum(weights) / num_samples
    end

    # `joint_cont[i, j, k, w]' is the total weight of samples whose i- and j-th variable are k and w, respectively
    joint_cont = pairwise_marginal(data; weights, num_cats, pseudocount)

    # `marginal_cont[i, j]' is the total weight of sample whose i-th variable is j
    marginal_cont = zeros(Float, num_vars, num_cats)
    for i = 1:num_vars
        for j = 1:num_cats
            @inbounds marginal_cont[i,j] = joint_cont[i,i,j,j]
        end
    end

    # Compute mutual information
    MI = zeros(Float, num_vars, num_vars)
    for var1_idx = 1 : num_vars
        for var2_idx = var1_idx : num_vars
            @inbounds MI[var1_idx, var2_idx] = sum(
                xlogx.(joint_cont[var1_idx, var2_idx, :, :]) .-
                xlogy.(joint_cont[var1_idx, var2_idx, :, :], marginal_cont[var1_idx, :] .* marginal_cont[var2_idx, :]'))
        end
    end

    for var1_idx = 2 : num_vars
        for var2_idx = 1 : var1_idx - 1
            @inbounds MI[var1_idx, var2_idx] = MI[var2_idx, var1_idx]
        end
    end

    MI
end

# TODO categorical pairwise MI on GPU


text:
using LDUFacts: sqrtapp

@testset "lduupdate(A+u*u')" begin
    u = ones(5)
    v = [1.0, 2, 0, -1, 1]
    A = [-2.0 2 2 0 2; 2 1 1 0 1; 2 1 2 -1 1; 0 0 -1 -2 0; 2 1 1 0 -1]
    la = ldu(A, FullPivot())
    f = 0.5
    ldu_update!(la, v, f)
    @test norm(la.P' * (A + v * v' * f) * la.P - la.L * la.D * la.U) <= 1e-13
end

@testset "lduupdate(A+u*v'+v*u')" begin
    u = ones(5)
    v = [1.0, 2, 0, -1, 1]
    A = [-2.0 2 2 0 2; 2 1 1 0 1; 2 1 2 -1 1; 0 0 -1 -2 0; 2 1 1 0 -1]
    la = ldu(A, FullPivot())
    ldu_update!(la, u, v, 0.1)
    @test norm(la.P' * (A + (u * v' + v * u') * 0.1) * la.P - la.L * la.D * la.U) <= 1e-13
    # note: test fails when f = -1 - check!!!
    la = ldu(A, FullPivot())
    f = -1
    @test_throws SingularException ldu_update!(la, u, v, f)
    #@test_broken norm(la.P' * (A + (u * v' + v * u') * f) * la.P - la.L * la.D * la.U) <= 1e-13
end

@testset "sqrtapp" begin
    @test abs(sqrtapp(0.5) / sqrt(0.5) - 1) < 0.05
    @test abs(sqrtapp(3//2) / sqrt(1.5) - 1) < 0.05
end


text:
module TestExamples3DPart1

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "3d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "3D-Part1" begin

# Run basic tests
@testset "Examples 3D" begin
  # Linear scalar advection
  include("test_examples_3d_advection.jl")

  # Hyperbolic diffusion
  include("test_examples_3d_hypdiff.jl")

  # Compressible Euler
  include("test_examples_3d_euler.jl")

  # Compressible Euler with self-gravity
  include("test_examples_3d_eulergravity.jl")
end

@testset "Displaying components 3D" begin
  @test_nowarn include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"))

  # test both short and long printing formats
  @test_nowarn show(mesh); println()
  @test_nowarn println(mesh)
  @test_nowarn display(mesh)

  @test_nowarn show(equations); println()
  @test_nowarn println(equations)
  @test_nowarn display(equations)

  @test_nowarn show(solver); println()
  @test_nowarn println(solver)
  @test_nowarn display(solver)

  @test_nowarn show(solver.basis); println()
  @test_nowarn println(solver.basis)
  @test_nowarn display(solver.basis)

  @test_nowarn show(solver.mortar); println()
  @test_nowarn println(solver.mortar)
  @test_nowarn display(solver.mortar)

  @test_nowarn show(semi); println()
  @test_nowarn println(semi)
  @test_nowarn display(semi)

  @test_nowarn show(summary_callback); println()
  @test_nowarn println(summary_callback)
  @test_nowarn display(summary_callback)

  @test_nowarn show(amr_controller); println()
  @test_nowarn println(amr_controller)
  @test_nowarn display(amr_controller)

  @test_nowarn show(amr_callback); println()
  @test_nowarn println(amr_callback)
  @test_nowarn display(amr_callback)

  @test_nowarn show(stepsize_callback); println()
  @test_nowarn println(stepsize_callback)
  @test_nowarn display(stepsize_callback)

  @test_nowarn show(save_solution); println()
  @test_nowarn println(save_solution)
  @test_nowarn display(save_solution)

  @test_nowarn show(analysis_callback); println()
  @test_nowarn println(analysis_callback)
  @test_nowarn display(analysis_callback)

  @test_nowarn show(alive_callback); println()
  @test_nowarn println(alive_callback)
  @test_nowarn display(alive_callback)

  @test_nowarn println(callbacks)
end

@testset "Additional tests in 3D" begin
  @testset "compressible Euler" begin
    eqn = CompressibleEulerEquations3D(1.4)

    @test isapprox(energy_total([1.0, 2.0, 3.0, 4.0, 20.0], eqn), 20.0)
    @test isapprox(energy_kinetic([1.0, 2.0, 3.0, 4.0, 20], eqn), 14.5)
    @test isapprox(energy_internal([1.0, 2.0, 3.0, 4.0, 20], eqn), 5.5)
  end

  @testset "hyperbolic diffusion" begin
    @test_nowarn HyperbolicDiffusionEquations3D(nu=1.0)
    eqn = HyperbolicDiffusionEquations3D(nu=1.0)

    @test isapprox(initial_condition_sedov_self_gravity(collect(1:3), 4.5, eqn), zeros(4))
    @test isapprox(boundary_condition_sedov_self_gravity(collect(1:4), 1, 1, collect(11:13), 2.3, flux_central, eqn), [-1.0, -19.739208802178712, 0.0, 0.0])
    @test isapprox(boundary_condition_sedov_self_gravity(collect(1:4), 2, 2, collect(11:13), 4.5, flux_central, eqn), [-1.5, 0.0, -19.739208802178712, 0.0])
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # 3D-Part1

end #module


text:
using Revise
using DynamicalSystems
using Basins
using Plots
using BenchmarkTools

function M(x)
    y=0.
     if x ≤ 0; y=9*x/(4-5*x); end
     if 0 < x ≤ 4/9; y=9/4*x; end
     if 4/9 < x ≤ 5/9; y=81/4*(x-x^2)-4; end
     if  x > 5/9; y=9/4*(1-x); end
     return y
 end

# Basin bifurcation in quasiperiodically forced systems Ulrike Feudel, Annette Witt, Ying-Cheng Lai, and Celso Grebogi PRE 28, 1998
function chaotic_map(dz,z, p, n)

    xn = z[1]; yn = z[2]
    λ=p[1];
    dz[1]=M(xn)
    dz[2]=λ*yn+sin(2*π*xn)
    return
end

# dummy function to keep the initializator happy
function chaotic_map_J(J,z0, p, n)
   return
end
ds = DiscreteDynamicalSystem(chaotic_map,[1., 0.], [1.1] , chaotic_map_J)
integ  = integrator(ds)

function escape_function(y)

    if y > 10000
        return 1
    elseif y < -10000
        return 2
    end
    return 0
end

function  get_color(integ,x,y)
    reinit!(integ,[x,y])
    while escape_function(integ.u[2]) == 0
        step!(integ)
        integ.t > 2000 && break
    end
    return escape_function(integ.u[2]);
end

xg=range(0.,1.,length=600)
yg=range(-2.,6.,length=600)
integ.p[1] = 1.1
basin_t = zeros(length(xg),length(yg))

for (k,x) in enumerate(xg), (n,y) in enumerate(yg)
    basin_t[k,n]=get_color(integ,x,y);
end

@show bd = box_counting_dim(xg, yg, basin_t)
#
# @show b,d = Basins.uncertainty_dimension_sample_nr(xg,yg,basin_t)
#
# using LsqFit
# # get exponent
# @. model(x, p) = p[1]*x+p[2]
# fit = curve_fit(model, vec(log10.(b)), vec(log10.(d)), [2., 2.])
# @show D = coef(fit)
# @show bd
# @show 2-D[1]

#@time bsn=Basins.basin_general_ds(xg, θ, integ; dt=1, Ncheck=10)
#plot(xg,yg,basin_t',seriestype=:heatmap)

# Estimation using the method of the paper.
# The results do not match exactly. But the paper talks about
# approximate results
N=[]
eps_r =  10 .^range(-8,0,length=10)

for eps in eps_r
    @show eps
    N2=0
    N1=0

    for k in 1:8192*20

        k1 = rand(1:length(xg))
        k2 = rand(1:length(yg))

        c1 = get_color(integ, xg[k1],yg[k2])
        c2 = get_color(integ, xg[k1]+eps,yg[k2])
        c3 = get_color(integ, xg[k1]-eps, yg[k2])

        if length(unique([c1,Int(c2),Int(c3)]))>1
            N1 += 1
        end
        N2 += 1
    end
    @show N1/N2
    push!(N,N1/N2)
end

using LsqFit
# get exponent
@. model(x, p) = p[1]*x+p[2]
fit = curve_fit(model, vec(log10.(eps_r)), vec(log10.(N)), [2., 2.])
D = coef(fit)
@show 2-D[1]


text:
using ApproxFun, Base.Test

ef = Fun(exp,Interval())

@test ef == -(-ef)
@test ef == (ef-1) + 1

ef = Fun(exp)

@test ef == -(-ef)
@test ef == (ef-1) + 1

cf = Fun(cos)

ecf = Fun(x->cos(x).*exp(x))
eocf = Fun(x->cos(x)./exp(x))

@test_approx_eq ef(.5) exp(.5)
@test_approx_eq ecf(.123456) cos(.123456).*exp(.123456)

r=2.*rand(100) .- 1

@test maximum(abs(ef(r)-exp(r)))<100eps()
@test maximum(abs(ecf(r)-cos(r).*exp(r)))<100eps()

@test norm((ecf-cf.*ef).coefficients)<100eps()

@test maximum(abs((eocf-cf./ef).coefficients))<1000eps()

@test norm(((ef/3).*(3/ef)-1).coefficients)<1000eps()

## Diff and cumsum

@test norm((ef - ef').coefficients)<10E-11

@test norm((ef - cumsum(ef)').coefficients) < 20eps()
@test norm((cf - cumsum(cf)').coefficients) < 20eps()

@test_approx_eq sum(ef)  2.3504023872876028

@test_approx_eq norm(ef)  1.90443178083307

##Check other domains

ef = Fun(exp,[1,2])
cf = Fun(cos,[1,2])

ecf = Fun(x->cos(x).*exp(x),[1,2])
eocf = Fun(x->cos(x)./exp(x),[1,2])

r=rand(100) .+ 1
x=1.5

@test_approx_eq ef(x) exp(x)

@test maximum(abs(ef(r)-exp(r)))<100eps()
@test maximum(abs(ecf(r)-cos(r).*exp(r)))<100eps()

@test norm((ecf-cf.*ef).coefficients)<100eps()

@test maximum(abs((eocf-cf./ef).coefficients))<1000eps()

@test norm(((ef/3).*(3/ef)-1).coefficients)<1000eps()

## Diff and cumsum

@test norm((ef - ef').coefficients)<10E-11

@test norm((ef - cumsum(ef)').coefficients) < 10eps()
@test norm((cf - cumsum(cf)').coefficients) < 10eps()

@test_approx_eq sum(ef) 4.670774270471604

@test_approx_eq norm(ef) 4.858451087240335

##Roots

f=Fun(x->sin(10(x-.1)))
@test norm(f(roots(f)))< 1000eps()

## ALiasing

f=Fun(x->cos(50acos(x)))
@test norm(f.coefficients-eye(ncoefficients(f))[:,51])<100eps()

## Int values

@test_approx_eq Fun(x->2,10)(.1) 2
@test_approx_eq Fun(x->2)(.1) 2

@test_approx_eq Fun(Float64[],Chebyshev)([0.,1.]) [0.,0.]
@test_approx_eq Fun([],Chebyshev)(0.) 0.
@test_approx_eq Fun(x->[1.,0.])(0.) [1.,0.]


text:
# https://learnxinyminutes.com/docs/julia/
# Single line comments start with a hash (pound) symbol.
#= Multiline comments can be written
   by putting '#=' before the text  and '=#'
   after the text. They can also be nested.
=#

####################################################
## 1. Primitive Datatypes and Operators
####################################################

# Everything in Julia is an expression.

# There are several basic types of numbers.
3 # => 3 (Int64)
3.2 # => 3.2 (Float64)
2 + 1im # => 2 + 1im (Complex{Int64})
2//3 # => 2//3 (Rational{Int64})


text:
"""
    Transformation(source_crs, target_crs; always_xy=false, direction = PJ_FWD, area=C_NULL, ctx=C_NULL)

Create a Transformation that is a pipeline between two known coordinate reference systems.
Transformation implements the
[CoordinateTransformations.jl](https://github.com/JuliaGeometry/CoordinateTransformations.jl)
API.

To do the transformation on coordinates, call an instance of this struct like a function.
See below for an example. These functions accept either 2 to 4 numbers, or 1 collection of 2
to 4 numbers, with the first two required numbers representing x and y coordinates. Number 3
is the elevation (defaults to 0), and number 4 is the time in years (defaults to Inf).

`source_crs` and `target_crs` can be:
- a "AUTHORITY:CODE", like EPSG:25832. When using that syntax for a source CRS, the created
  pipeline will expect that the coordinates respect the axis order and axis unit of the
  official definition (so for example, for EPSG:4326, with latitude first and longitude
  next, in degrees). Similarly, when using that syntax for a target CRS, output values will
  be emitted according to the official definition of this CRS. This behavior can be
  overruled by passing `always_xy=true`.
- a PROJ string, like "+proj=longlat +datum=WGS84". When using that syntax, the axis order
  and unit for geographic CRS will be longitude, latitude, and the unit degrees.
- the name of a CRS as found in the PROJ database, e.g "WGS84", "NAD27", etc.
- more generally any string accepted by `proj_create` representing a CRS
- besides an `AbstractString`, it can also accept a `Ptr{PJ}`, pointing to a CRS that was
  already created with `proj_create`

`always_xy` can optionally fix the axis orderding to x,y or lon,lat order. By default it is
`false`, meaning the order is defined by the authority in charge of a given coordinate
reference system, as explained in [this PROJ FAQ
entry](https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent).

`direction` can be one of (`PJ_FWD`, `PJ_IDENT`, `PJ_INV`), which correspond to forward
(source to target), identity (do nothing) and inverse (target to source) transformations.

`area` sets the "area of use" for the Transformation. When it is supplied, the more accurate
transformation between two given systems can be chosen. When no area of use is specific and
several coordinate operations are possible depending on the area of use, this function will
internally store those candidate coordinate operations in the return PJ object. Each
subsequent coordinate transformation will then select the appropriate coordinate operation
by comparing the input coordinates with the area of use of the candidate coordinate
operations. The `area` pointer needs to be created using `proj_area_create`, filled using
`proj_area_set_bbox`, and destroyed using `proj_area_destroy`.

`ctx` determines the threading context. By default it is set to the global context. For
thread safety, use separate contexts created with `proj_context_create` or
`proj_context_clone`, and destroyed with `proj_context_destroy`.

# Examples
```julia
julia> trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy=true)
Transformation
    source: WGS 84 (with axis order normalized for visualization)
    target: Amersfoort / RD New

julia> trans(5.39, 52.16)  # this is in lon,lat order, since we set always_xy to true
(155191.3538124342, 463537.1362732911)
```
"""
mutable struct Transformation <: CoordinateTransformations.Transformation
    pj::Ptr{PJ}
    direction::PJ_DIRECTION
    function Transformation(pj::Ptr{PJ}, direction::PJ_DIRECTION = PJ_FWD)
        trans = new(pj, direction)
        finalizer(trans) do trans
            trans.pj = proj_destroy(trans.pj)
        end
        return trans
    end

end

function Transformation(
    source_crs::AbstractString,
    target_crs::AbstractString;
    always_xy::Bool = false,
    direction::PJ_DIRECTION = PJ_FWD,
    area::Ptr{PJ_AREA} = C_NULL,
    ctx::Ptr{PJ_CONTEXT} = C_NULL,
)
    pj = proj_create_crs_to_crs(source_crs, target_crs, area, ctx)
    pj = always_xy ? normalize_axis_order!(pj; ctx) : pj
    return Transformation(pj, direction)
end

function Transformation(
    source_crs::Ptr{PJ},
    target_crs::Ptr{PJ};
    always_xy::Bool = false,
    direction::PJ_DIRECTION = PJ_FWD,
    area::Ptr{PJ_AREA} = C_NULL,
    ctx::Ptr{PJ_CONTEXT} = C_NULL,
)
    pj = proj_create_crs_to_crs_from_pj(source_crs, target_crs, area, ctx)
    pj = always_xy ? normalize_axis_order!(pj; ctx) : pj
    return Transformation(pj, direction)
end

function Base.show(io::IO, trans::Transformation)
    source_crs = proj_get_source_crs(trans.pj)
    target_crs = proj_get_target_crs(trans.pj)
    source_info = proj_pj_info(source_crs)
    target_info = proj_pj_info(target_crs)
    source_description = unsafe_string(source_info.description)
    target_description = unsafe_string(target_info.description)

    dir = trans.direction
    direction_str = dir == PJ_FWD ? "forward" : dir == PJ_IDENT ? "identity" : "inverse"

    print(
        io,
        """Transformation
            source: $source_description
            target: $target_description
            direction: $direction_str
        """,
    )
end

"""
    normalize_axis_order!(pj::Ptr{PJ}; ctx = C_NULL)

Call proj_normalize_for_visualization on an object, and return the new object after freeing
the input object.
"""
function normalize_axis_order!(pj::Ptr{PJ}; ctx = C_NULL)
    pj_for_gis = proj_normalize_for_visualization(pj, ctx)
    proj_destroy(pj)
    return pj_for_gis
end

function Base.inv(
    trans::Transformation;
    always_xy::Bool = false,
    area::Ptr{PJ_AREA} = C_NULL,
    ctx::Ptr{PJ_CONTEXT} = C_NULL,
)
    source_crs = proj_get_source_crs(trans.pj)
    target_crs = proj_get_target_crs(trans.pj)

    return Transformation(
        source_crs,
        target_crs;
        direction = inv(trans.direction),
        area,
        ctx,
        always_xy,
    )
end

function (trans::Transformation)(x, y)::NTuple{2,Float64}
    p = proj_trans(trans.pj, trans.direction, (x, y))
    return p.x, p.y
end

function (trans::Transformation)(x, y, z)::NTuple{3,Float64}
    p = proj_trans(trans.pj, trans.direction, (x, y, z))
    return p.x, p.y, p.z
end

function (trans::Transformation)(x, y, z, t)::NTuple{4,Float64}
    p = proj_trans(trans.pj, trans.direction, (x, y, z, t))
    return p.x, p.y, p.z, p.t
end

function (trans::Transformation)(coord::Coord)::NTuple{4,Float64}
    p = proj_trans(trans.pj, trans.direction, coord)
    return p.x, p.y, p.z, p.t
end

function (trans::Transformation)(coord)::NTuple234
    n = length(coord)
    p = proj_trans(trans.pj, trans.direction, coord)
    return if n == 2
        p.x, p.y
    elseif n == 3
        p.x, p.y, p.z
    else
        p.x, p.y, p.z, p.t
    end
end

function CoordinateTransformations.compose(
    trans1::Transformation,
    trans2::Transformation;
    always_xy::Bool = false,
    area::Ptr{PJ_AREA} = C_NULL,
    ctx::Ptr{PJ_CONTEXT} = C_NULL,
)
    # create a new Transformation from trans1 source to trans2 target
    # can also be typed as trans1 ∘ trans2, typed with \circ
    # a → b ∘ c → d doesn't make much sense if b != c, though we don't enforce it

    source_crs =
        trans1.direction == PJ_FWD ? proj_get_source_crs(trans1.pj) :
        proj_get_target_crs(trans1.pj)
    target_crs =
        trans2.direction == PJ_FWD ? proj_get_target_crs(trans2.pj) :
        proj_get_source_crs(trans2.pj)
    return Transformation(source_crs, target_crs; area, ctx, always_xy)
end

"""
    enable_network!(active::Bool = true, ctx::Ptr{PJ_CONTEXT} = C_NULL)::Bool

Enable PROJ network access, if `active` is true, disable it if it is false. Optionally pass
a context to set it for that context, instead of the global one.

Returns true if network access is possible.
"""
function enable_network!(active::Bool = true, ctx::Ptr{PJ_CONTEXT} = C_NULL)
    enabled = proj_context_set_enable_network(Cint(active), ctx)
    return Bool(enabled)
end

"""
    network_enabled(ctx::Ptr{PJ_CONTEXT} = C_NULL)::Bool

Returns true if PROJ network access is enabled, false otherwise. Optionally pass a context
to check for that context, instead of the global one.
"""
function network_enabled(ctx::Ptr{PJ_CONTEXT} = C_NULL)
    enabled = proj_context_is_network_enabled(ctx)
    return Bool(enabled)
end

function Base.inv(direction::PJ_DIRECTION)
    return if direction == PJ_FWD
        PJ_INV
    elseif direction == PJ_IDENT
        PJ_IDENT
    else
        PJ_FWD
    end
end


text:
"""
    mutable struct Compressibility <: AbstractConstraint

    Compressibility(inv_zeta = 0.0)

A local incompressibility constraint on each grid point in a system.

Parameter `inv_zeta` controls how incompressibility is treated:
* ``ζ^{-1} = 0.0``: System is treated as incompressible
* ``ζ^{-1} > 0.0``: System is treated as compressible with bulk modulus ``ζ``
* ``ζ^{-1} < 0.0``: Incompressibility is not enforced
"""
mutable struct Compressibility <: AbstractConstraint
    inv_zeta :: Float64
    strict   :: Bool
    field    :: FieldGrid{Float64}
    rhosum   :: FieldGrid{Float64}
    system   :: Option{AbstractSystem}

    function Compressibility(inv_zeta::Real = 0.0)
        strict = inv_zeta ≈ 0.0 ? true : false
        return new(inv_zeta, strict, zeros(0,0,0), zeros(0,0,0), nothing)
    end
end

#==============================================================================#

Base.show(io::IO, c::Compressibility) = @printf(io, "Compressibility(inv_zeta = %.3f)", c.inv_zeta)

function setup!(c::Compressibility, sys::AbstractSystem)
    c.system = sys
    c.field = zeros(Float64, sys.dims)
    c.rhosum = zeros(Float64, sys.dims)
    return nothing
end

function update!(c::Compressibility)
    @assert !isnothing(c.system)
    sys = c.system

    # Update density sum
    c.rhosum .= 0.0
    for (mid, rho) in sys.density
        @. c.rhosum += rho
    end

    c.field .= 0.0
    if c.strict
        # Incmpressible system
        @. c.field = c.rhosum - 1.0
        for (mid, omega) in sys.fields
            @. c.field += omega - sys.potentials[mid]
        end
        c.field ./= nmonomers(sys)
    elseif c.inv_zeta > 0.0
        # Compressible system
        @. c.field = (1.0/c.inv_zeta) * (c.rhosum - 1.0)
    end

    return nothing
end


text:
"""
    DataArray

A type wrapping a predefined array of measurement data.
"""
immutable DataArray{T<:Number,D<:Number,A<:AbstractArray} <: AbstractData
    index::Vector
    values::A
    DataArray(i::AbstractVector{T},d::AbstractArray{D}) = new(collect(i),d)
end

"""
    data(:array,index::AbstractVector,values::AbstractArray)

Create a data object wrapping an array of fixed, predefined data values.

# Arguments
* `index::AbstractVector`: an index for each row of the values array (e.g., a time index)
* `values::AbstractArray`: each row is a data point and each column a variable
---
"""
function data(::Type{Val{:array}},i::AbstractVector,d::AbstractArray)
    @assert length(i) == size(d,1)
    DataArray{eltype(i),eltype(d),typeof(d)}(i,d)
end

numvalues(d::DataArray) = size(d.values,1)
numvars(d::DataArray) = size(d.values,2)
eltype(d::DataArray) = eltype(d.values)
generate!(d::DataArray) = d
dataindex(d::DataArray) = d.index
datavalues(d::DataArray) = d.values
datatypename(d::DataArray) = "DataArray"


text:
#NOTE: We do the error handling of nnp_initialize while loading NNPACK
function nnp_initialize()
    ccall((:nnp_initialize, libnnpack), nnp_status, (),)
end

function nnp_deinitialize()
    @nnpack_check ccall((:nnp_deinitialize, libnnpack), nnp_status, (),)
end

function pthreadpool_create(n = 0)
    ccall((:pthreadpool_create, libnnpack), Ptr{Cvoid}, (Csize_t,), n)
end

function nnp_relu_output(batch_size, channels, input, output, negative_slope, threadpool)
    @nnpack_check ccall((:nnp_relu_output, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, pthreadpool_t), batch_size, channels, input, output, negative_slope, threadpool)
end

function nnp_relu_output(x::Array{Float32,N}, y::Array{Float32,N}; negative_slope::AbstractFloat = 0.0, threadpool = C_NULL) where {N}
    # Investigate why the channel and batch dims need to specified like this
    nnp_relu_output(prod(size(x)[N-1:N]), prod(size(x)[1:N-2]), x, y, negative_slope, threadpool)
    y
end

function nnp_relu_input_gradient(batch_size, channels, grad_output, input, grad_input, negative_slope, threadpool)
    @nnpack_check ccall((:nnp_relu_input_gradient, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, pthreadpool_t), batch_size, channels, grad_output, input, grad_input, negative_slope, threadpool)
end

function nnp_relu_input_gradient(x::Array{Float32,N}, dy::Array{Float32,N}, dx::Array{Float32,N}; negative_slope::AbstractFloat = 0.0, threadpool = C_NULL) where {N}
    # Investigate why the channel and batch dims need to specified like this
    nnp_relu_input_gradient(Csize_t(prod(size(x)[N-1:N])), prod(size(x)[1:N-2]), dy, x, dx, negative_slope, threadpool)
    dx
end

function nnp_softmax_output(batch_size, channels, input, output, threadpool)
    @nnpack_check ccall((:nnp_softmax_output, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), batch_size, channels, input, output, threadpool)
end

function nnp_softmax_output(x::VecOrMat{Float32}, y::VecOrMat{Float32}; threadpool = C_NULL)
    nnp_softmax_output(ndims(x) == 2 ? size(x, 2) : 1, size(x, 1), x, y, threadpool)
    y
end

#FIXME: Output of fully connected not consistent with `kernel * input`
#NOTE: This most likely due to nnpack being row major. Investigate this.

function nnp_fully_connected_output(batch_size, input_channels, output_channels, input, kernel, output, threadpool, profile)
    @nnpack_check ccall((:nnp_fully_connected_output, libnnpack), nnp_status, (Csize_t, Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t, Ptr{Cvoid}), batch_size, input_channels, output_channels, input, kernel, output, threadpool, C_NULL)
end

function nnp_fully_connected_output(x::Array{Float32,2}, w::Array{Float32,2}, y::Array{Float32,2}; profile = nothing, threadpool = C_NULL)
    profile = profile == nothing ? nnp_profile() : profile
    nnp_fully_connected_output(size(x, 2), size(x, 1), size(w, 1), x, w, y, threadpool, profile)
    y
end

function nnp_fully_connected_inference_f16f32(input_channels, output_channels, input, kernel, output, threadpool)
    @nnpack_check ccall((:nnp_fully_connected_inference_f16f32, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cvoid}, Ptr{Cfloat}, pthreadpool_t), input_channels, output_channels, input, kernel, output, threadpool)
end

nnp_fully_connected_inference_f16f32(x::Array{Float32, 1}, w::Array{Float16,2}, y::Array{Float32, 1}; threadpool = C_NULL) =
    nnp_fully_connected_inference(reshape(x, size(x), 1), w, reshape(y, size(y), 1), threadpool = threadpool)

function nnp_fully_connected_inference_f16f32(x::Array{Float32, 2}, w::Array{Float16,2}, y::Array{Float32, 2}; threadpool = C_NULL)
    nnp_fully_connected_inference(size(x, 1), size(y, 1), x, w, y, threadpool)
    y
end

function nnp_fully_connected_inference(input_channels, output_channels, input, kernel, output, threadpool)
    @nnpack_check ccall((:nnp_fully_connected_inference, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), input_channels, output_channels, input, kernel, output, threadpool)
end

nnp_fully_connected_inference(x::Array{Float32, 1}, w::Array{Float32,2}; threadpool = C_NULL) =
    nnp_fully_connected_inference(reshape(x, size(x), 1), w, threadpool = threadpool)

function nnp_fully_connected_inference(x::Array{Float32, 2}, w::Array{Float32, 2}, y::Array{Float32, 2}; threadpool = C_NULL)
    nnp_fully_connected_inference(size(x, 1), size(y, 1), x, w, y, threadpool)
    y
end

function nnp_max_pooling_output(batch_size, channels, input_size, input_padding, pooling_size, pooling_stride, input, output, threadpool)
    @nnpack_check ccall((:nnp_max_pooling_output, libnnpack), nnp_status, (Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), batch_size, channels, input_size, input_padding, pooling_size, pooling_stride, input, output, threadpool)
end

function nnp_max_pooling_output(y::Array{Float32,4}, x::Array{Float32,4}, kernel::Tuple; padding = 0, stride = 1, threadpool = C_NULL)
    input_size = nnp_size(Csize_t.((size(x, 1), size(x, 2)))...)
    pooling_size = nnp_size(Csize_t.(kernel)...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    pooling_stride = nnp_size(Csize_t.(stride)...)
    nnp_max_pooling_output(size(x, 4), size(x, 3), input_size, input_padding, pooling_size, pooling_stride, x, y, threadpool)
    y
end

#TODO: Add wrapper for convolution inference

function nnp_convolution_input_gradient(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, grad_output, kernel, grad_input, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @nnpack_check ccall((:nnp_convolution_input_gradient, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, grad_output, kernel, grad_input, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_input_gradient(dx::Array{Float32,4}, dy::Array{Float32,4}, w::Array{Float32,4}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = C_NULL, profile = nothing)
    input_size = nnp_size(Csize_t.((size(dx,1), size(dx,2)))...)
    kernel_size = nnp_size(Csize_t.((size(w,1),size(w,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_input_gradient(UInt32(algo), size(dx,4), size(dx,3), size(w,4), input_size, input_padding, kernel_size, dy, w, dx, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    dx
end

function nnp_convolution_kernel_gradient(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, grad_output, grad_kernel, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @nnpack_check ccall((:nnp_convolution_kernel_gradient, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, grad_output, grad_kernel, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_kernel_gradient(dw::Array{Float32,4}, x::Array{Float32,4}, dy::Array{Float32,4}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = C_NULL, profile = nothing)
    input_size = nnp_size(Csize_t.((size(x,1), size(x,2)))...)
    kernel_size = nnp_size(Csize_t.((size(dw,1),size(dw,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_kernel_gradient(UInt32(algo), size(x,4), size(x,3), size(dw,4), input_size, input_padding, kernel_size, x, dy, dw, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    dw
end

function nnp_convolution_output(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, kernel, bias, output, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @nnpack_check ccall((:nnp_convolution_output, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, kernel, bias, output, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_output(y::Array{Float32,4}, x::Array{Float32,4}, w::Array{Float32,4}, b::Array{Float32,1}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = C_NULL, profile = nothing)
    input_size = nnp_size(Csize_t.((size(x,1), size(x,2)))...)
    kernel_size = nnp_size(Csize_t.((size(w,1),size(w,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_output(UInt32(algo), size(x,4), size(x,3), size(w,4), input_size, input_padding, kernel_size, x, w, b, y, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    y
end


text:
#=
Detection methods. Each method needs
1) Type with fields necessary for algorithm
2) function "detect" to operate on sort with type field defined above
3) any other necessary functions for detection algorithm

A method may also need a "detectprepare" function to use in its first iteration if the detection method uses a sliding window that depends on previous iterations (see power for an example)
=#

export DetectAbs, DetectNeg

function detectprepare(d::Detect,sort::Sorting,v::V) where V
end

#=
Power

Rutishauser et al 2006
=#
#=
type DetectPower <: Detect
    a::Int64
    b::Int64
    c::Int64
end

DetectPower()=DetectPower(0,0,0)

function detect{V}(d::DetectPower,sort::Sorting, i::Int64,v::V)

    @inbounds sort.d.a += v[i,sort.id] - sort.d.c
    @inbounds sort.d.b += v[i,sort.id]^2 - sort.d.c^2

    if i>=power_win
        @inbounds sort.d.c=v[i-power_win0,sort.id]
    else
        @inbounds sort.d.c=sort.sigend[i+sort.s.s_end-power_win0]
    end

    # equivalent to p = sqrt(1/n * sum( (f(t-i) - f_bar(t))^2))
    sort.thres < sqrt((sort.d.b - (sort.d.a^2/power_win))/power_win)
end

function detectprepare{V}(d::DetectPower,sort::Sorting,v::V)

    sort.d.a=0
    sort.d.b=0

    for i=1:power_win
        sort.d.a += v[i,sort.id]
        sort.d.b += v[i,sort.id]^2
    end

    sort.d.c=v[1,sort.id]

    for i=power_win1:sort.s.s_end
        sort.d.a += v[i,sort.id] - sort.d.c
        sort.d.b += v[i,sort.id]^2 - sort.d.c^2
        sort.d.c = v[i-power_win+1,sort.id]
    end

    nothing
end
=#
#=
Absolute Value

Quiroga et al 2004
=#

mutable struct DetectAbs <: Detect
end

function detect(d::DetectAbs,sort::Sorting,i::Int64,v)
    @inbounds return sort.thres < abs(v[i,sort.id])
end

#=
Negative Value
=#

mutable struct DetectNeg <: Detect
end

function detect(d::DetectNeg,sort::Sorting,i::Int64,v)
    @inbounds return sort.thres > v[i,sort.id]
end

#=
Nonlinear Energy Operator

Choi et al 2006
=#
#=
type DetectNEO <: Detect
end

function detect(d::DetectNEO, sort::Sorting,i::Int64, v::AbstractArray{Int64,2})

    if i==length(sort.rawSignal)
        #Will do spike detection next iteration due to edging
        psi=0
    elseif i>1
        psi=sort.rawSignal[i]^2 - sort.rawSignal[i+1] * sort.rawSignal[i-1]
    else

        #perform calculation for end of last step and this one, and return the larger
        psi1=sort.sigend[end]^2 - sort.rawSignal[i] * sort.sigend[end-1]
        psi2=sort.rawSignal[i]^2 - sort.rawSignal[i+1] * sort.sigend[end]

        if psi1>psi2
            return psi1
        else
            return psi2
        end
    end

    psi
end
=#
#=
Continuous Wavelet Transform

Nenadic et al 2005, Benitez et al 2008
=#

#=
Discrete Wavelet Transform

Kim and Kim 2003
=#

#=
Nonparametric Detection

Song et al 2006
=#

#=
Smoothed Nonlinear Energy Operator

Azami et al 2014
=#
#=
type DetectSNEO <: Detect
end

function detect(d::DetectSNEO, sort::Sorting,i::Int64, v::AbstractArray{Int64,2})

    if i==length(sort.rawSignal)
        #Will do spike detection next iteration due to edging
        psi=0
    elseif i>20

    elseif i>1
        psi=sort.rawSignal[i]^2 - sort.rawSignal[i+1] * sort.rawSignal[i-1]
    else

        #perform calculation for end of last step and this one, and return the larger
        psi1=sort.sigend[end]^2 - sort.rawSignal[i] * sort.sigend[end-1]
        psi2=sort.rawSignal[i]^2 - sort.rawSignal[i+1] * sort.sigend[end]

        if psi1>psi2
            return psi1
        else
            return psi2
        end
    end

    psi
end
=#

#=
Multiscale Correlation of Wavelet Coefficients

Yang et al 2011
=#

#=
type DetectMCWC <: Detect
    Tx::Array{Float64,1}
    rs::Array{Float64,1}
end

function DetectMCWC()
    DetectMCWC(zeros(Float64,11),zeros(Float64,10))
end

function detect(d::DetectMCWC, sort::Sorting, i::Int64, v::AbstractArray{Int64,2})

    p=0.0

    #calculate wavelet coefficients
    if i<=bigJ
        sort.d.Tx[:]=coiflets_scaled*[sort.sigend[(end-bigJ+i+1):end];sort.rawSignal[1:i]]
    else
        sort.d.Tx[:]=coiflets_scaled*sort.rawSignal[(i-bigJ+1):i]
    end

    for j=1:11
        sort.d.Tx[j] *= onesquarea[j]
    end

    #Correlation of wavelet coefficients
    for j=1:10
        sort.d.rs[j] = sort.d.Tx[j]*sort.d.Tx[j+1]
    end

    #Since I'm using a rolling window, I do this a little differently
    #maybe pre calculate the power to normalize?
    for j=1:10
        if sort.d.Tx[j]^2 < abs(sort.d.rs[j])
            p=2.0
            break
        end
    end

    for j=1:11
        sort.d.Tx[j]=0.0
    end

    p

end

=#
#=

=#


text:
# verify that our README.md examples are in working order

doc = Spiderman.http_get("http://en.wikipedia.org/")
news_headlines = text(select(css"#mp-itn b a", doc))

# most of the value in this test is the above code running...
# but it's still nice to know the query still works.
@test length(news_headlines) > 1

doc = Spiderman.http_get("https://news.ycombinator.com/")

scores = select(css".score", doc)
hrefs = parent(select(css".deadmark", doc))

items = []
for i=1:length(scores)
    score = parse(Int, replace(text(scores[i]), " points", ""))
    story_title = text(select(css"a", hrefs[i]))[1]
    story_href = getattr(select(css"a", hrefs[i])[1], "href")
    comment_href = getattr(select(css"a", parent(scores[i]))[end], "href")

    push!(items, (score, story_title, story_href, comment_href))
end

@test length(items) > 2


text:
@testset "Long Sequence" begin
    str = random_dna(50000)
    @test ReferenceSequence(str) == DNASequence(str)
    str = ("N"^300 * random_dna(1000))^5 * "N"^300
    @test ReferenceSequence(str) == DNASequence(str)
end