text:
module Schemas

import Base.==
import Base.fullname
import Base.hash
import Base.show

using Avro.Common
using JSON

export Schema,
       NullSchema,
       BooleanSchema,
       StringSchema,
       BytesSchema,
       IntSchema,
       LongSchema,
       FloatSchema,
       DoubleSchema,
       RecordSchema,
       EnumSchema,
       FixedSchema,
       ArraySchema,
       MapSchema,
       UnionSchema,
       PRIMITIVE_TYPES,
       parse

"""
The set of builtin types.
"""
const PRIMITIVE_TYPES = [
    "null",
    "boolean",
    "int",
    "long",
    "float",
    "double",
    "bytes",
    "string"
    ]

"""
The set of named types.
"""
const NAMED_TYPES = [
    "record",
    "enum",
    "fixed"
    ]

"""
Valid types may appear in the `name` field of a JSON object.
"""
const VALID_TYPES = [
    PRIMITIVE_TYPES;
    NAMED_TYPES;
    ["array", "map", "union"]
    ]

function get_required(data::Dict, key::String, error_message::String)
    get(data, key) do
        throw(SchemaParseException(string(error_message, ": ", data)))
    end
end

"""
Exception during Schema parsing.
"""
struct SchemaParseException <: Exception
    message::String
end

"""
The parent of the Avro Schema hierarchy.
"""
abstract type Schema end;

abstract type PrimitiveSchema <: Schema end;

# Generate the primitive type schemas
for primitive_type in PRIMITIVE_TYPES
    classname = Symbol(capitalize(primitive_type), "Schema")
    @eval begin
        struct $(classname) <: PrimitiveSchema
        end
    end
    @eval const $(Symbol(uppercase(primitive_type))) = $(classname)()
end

function PrimitiveSchema(typename::String)
    if typename == "null"
        return NULL
    elseif typename == "boolean"
        return BOOLEAN
    elseif typename == "string"
        return STRING
    elseif typename == "bytes"
        return BYTES
    elseif typename == "int"
        return INT
    elseif typename == "long"
        return LONG
    elseif typename == "float"
        return FLOAT
    elseif typename == "double"
        return DOUBLE
    else
        throw(Exception("Invalid schema typename"))
    end
end

"""
A fully qualified name in Avro.
"""
struct FullName
    value::String
end

"""
Construct a fullname from a name and space.
"""
function FullName(name::String, space::String)
    if '.' in name
        FullName(name)
    else
        FullName(space == "" ? name : string(space, '.', name))
    end
end

"""
A context object for maintaining state during schema parsing. The the state
includes the current namespace and a dictionary of the parsed schemas.
"""
mutable struct ParseContext
    space::String
    schemas::Dict{FullName, Schema}
end

"""
A schema that can be uniquely identified with a fullname. A fullname is
comprised of a name and a namespace.
"""
abstract type NamedSchema <: Schema end;

"""
Constructs named schemas from a JSON object. The named schema object returned
depends on the type provided by the JSON object.

    "fixed" | "error" => FixedSchema
    "enum" => EnumSchema
    "record" => RecordSchema

If a named schema with the specified name has already been parsed, the
previously parsed schema object is returned from the parse context.
"""
function NamedSchema(json_data::Dict, context::ParseContext, typename::String)
    # Push the current namespace onto the stack
    parent_space = context.space

    # Parse the name, namespace
    name = get_required(json_data, "name", "No name in schema")
    space = parent_space
    if haskey(json_data, "namespace")
        space = json_data["namespace"]
        # Set the current namespace if present
        context.space = space
    end
    fullname = FullName(name, space)

    # Parse the optional aliases and doc
    doc = get(json_data, "doc", "")
    alias_names = get(json_data, "aliases", [])
    aliases = [FullName(name, context.space) for name in alias_names]

    # Return early if a schema with the given fullname already exists
    if haskey(context.schemas, fullname)
        return context.schemas[fullname]
    end

    # Parse the correct schema
    if typename in ["record", "error"]
        schema = RecordSchema(json_data, context, fullname, doc, aliases)
    elseif typename == "fixed"
        schema = FixedSchema(json_data, context, fullname, doc, aliases)
    elseif typename == "enum"
        schema = EnumSchema(json_data, context, fullname, doc, aliases)
    end

    # Pop the namespace on the stack
    context.space = parent_space

    schema
end

"""
An enum of possible orders.
"""
@enum Order Ascending Descending Ignore

"""
Construct an Order object from a name. This constructor is case-insensitive. If
the name is invalid, throws an exception.
"""
function Order(name::String)
    clean_name = uppercase(name)
    if clean_name == "ASCENDING"
        Ascending
    elseif clean_name == "DESCENDING"
        Descending
    elseif clean_name == "IGNORE"
        Ignore
    else
        throw(Exception("Invalid order name: $name"))
    end
end

"""
A field of a record schema.
"""
struct Field
    name::String
    position::Int
    schema::Schema

    # TODO
    # default::T

    doc::String
    order::Order
    aliases::Vector{FullName}
end

function Field(
        name::String,
        position::Int,
        schema::Schema;
        doc::String = "",
        order::Order = Ascending,
        aliases::Vector{FullName} = FullName[])
    Field(name, position, schema, doc, order, aliases)
end

"""
Construct a field of a record given a JSON object.
"""
function Field(field_data::Dict, context::ParseContext, position::Int)
    name = get_required(field_data, "name", "No field name")
    doc = get(field_data, "doc", "")
    schema = parse_schema(field_data["type"], context)

    # TODO
    # default = field_data["default"]

    order = Order(get(field_data, "order", "ascending"))
    alias_names = get(field_data, "aliases", [])
    aliases = [FullName(name, context.space) for name in alias_names]

    Field(name, position, schema, doc, order, aliases)
end

"""
A record schema.
"""
mutable struct RecordSchema <: NamedSchema
    fullname::FullName
    fields::Vector{Field}
    doc::String
    aliases::Vector{FullName}
end

function RecordSchema(
        fullname::FullName,
        fields::Vector{Field};
        doc::String = "",
        aliases::Vector{FullName} = FullName[])
    RecordSchema(fullname, fields, doc, aliases)
end

function RecordSchema(fullname::FullName, doc::String, aliases::Vector{FullName})
    RecordSchema(fullname, Field[], doc, aliases)
end

"""
Construct a record schema from a JSON object.
"""
function RecordSchema(
        json_data::Dict,
        context::ParseContext,
        fullname::FullName,
        doc::String,
        aliases::Vector{FullName})
    schema = RecordSchema(fullname, doc, aliases)
    context.schemas[fullname] = schema

    # Parse the fields
    json_fields = json_data["fields"]
    fields = Array{Field}(undef, length(json_fields))
    for (i, json_field) in enumerate(json_fields)
        fields[i] = Field(json_field, context, i - 1)
    end
    schema.fields = fields

    schema
end

"""
An enum schema.
"""
struct EnumSchema <: NamedSchema
    fullname::FullName
    symbols::Vector{String}
    doc::String # optional
    aliases::Vector{FullName} # optional
end

function EnumSchema(
        fullname::FullName,
        symbols::Vector{String};
        doc::String = "",
        aliases::Vector{FullName} = FullName[])
    EnumSchema(fullname, symbols, doc, aliases)
end

"""
Construct an enum schema from a JSON object.
"""
function EnumSchema(
        json_data::Dict,
        context::ParseContext,
        fullname::FullName,
        doc::String,
        aliases::Vector{FullName})
    symbols = Array{String}(json_data["symbols"])
    schema = EnumSchema(fullname, symbols, doc, aliases)
    context.schemas[fullname] = schema
    schema
end

"""
An array schema.
"""
struct ArraySchema{T <: Schema} <: Schema
    items::T
end

"""
Construct an array schema from a JSON object.
"""
function ArraySchema(json_data::Dict, context::ParseContext)
    items = parse_schema(json_data["items"], context)
    ArraySchema(items)
end

"""
A map schema.
"""
struct MapSchema{T <: Schema} <: Schema
    values::T
end

"""
Construct a map schema from a JSON object.
"""
function MapSchema(json_data::Dict, context::ParseContext)
    values = parse_schema(json_data["values"], context)
    MapSchema(values)
end

"""
A union schema.
"""
struct UnionSchema <: Schema
    schemas::Vector{Schema}
end

"""
Construct a union schema from a JSON array.
"""
function UnionSchema(json_data::Array, context::ParseContext)
    schemas = map(schema_data -> parse_schema(schema_data, context), json_data)
    UnionSchema(schemas)
end

"""
A fixed schema.
"""
struct FixedSchema <: NamedSchema
    fullname::FullName
    size::Int
    doc::String
    aliases::Vector{FullName}
end

function FixedSchema(
        fullname::FullName,
        size::Int;
        doc::String = "",
        aliases::Vector{FullName} = FullName[])
    FixedSchema(fullname, size, doc, aliases)
end

"""
Construct a fixed schema from a JSON object.
"""
function FixedSchema(
        json_data::Dict,
        context::ParseContext,
        fullname::FullName,
        doc::String,
        aliases::Vector{FullName})
    size = json_data["size"]
    schema = FixedSchema(fullname, size, doc, aliases)
    context.schemas[fullname] = schema
    schema
end

"""
Parse a schema object from a JSON object.
"""
function parse_schema(json_data::Dict, context::ParseContext)
    # Throws KeyError if not found
    typename = get_required(json_data, "type", "No type")
    if typename in PRIMITIVE_TYPES
        PrimitiveSchema(typename)
    elseif typename in NAMED_TYPES
        NamedSchema(json_data, context, typename)
    elseif typename in VALID_TYPES
        if typename == "array"
            ArraySchema(json_data, context)
        elseif typename == "map"
            MapSchema(json_data, context)
        end
    else
        throw(SchemaParseException("Schema not yet supported: $typename"))
    end
end

"""
Parse a schema object from a JSON object. This schema must be a Union according
to the Avro specification.
"""
function parse_schema(json_data::Array, context::ParseContext)
    UnionSchema(json_data, context)
end

"""
Parse a schema object from a parsed JSON string.
"""
function parse_schema(schema_data::String, context::ParseContext)
    if schema_data in PRIMITIVE_TYPES
        # Choose the primitive schema based on the type name
        PrimitiveSchema(schema_data)
    else
        fullname = FullName(schema_data, context.space)
        get(context.schemas, fullname) do
            fullname = FullName(schema_data)
            context.schemas[fullname]
        end
    end
end

"""
Parse a schema object from a JSON string.
"""
function parse(json_string::String)
    json_data = JSON.parse(json_string)

    # Three cases:
    # 1. JSON object
    # 2. JSON array (union)
    # 3. JSON string (primitive)
    context = ParseContext("", Dict())
    parse_schema(json_data, context)
end

"""
Equality definitions for fullnames.
"""
==(fullname1::FullName, fullname2::FullName) = fullname1.value == fullname2.value
hash(fullname::FullName) = hash(fullname.value)

"""
Equality definitions for schemas.
"""
==(::Schema, ::Schema) = false
==(::A, ::A) where A <: PrimitiveSchema = true
==(a::RecordSchema, b::RecordSchema) = a.fullname == b.fullname
==(a::EnumSchema, b::EnumSchema) = a.fullname == b.fullname && a.symbols == b.symbols
==(a::ArraySchema, b::ArraySchema) = a.items == b.items
==(a::MapSchema, b::MapSchema) = a.values == b.values
==(a::UnionSchema, b::UnionSchema) = a.schemas == b.schemas
==(a::FixedSchema, b::FixedSchema) = a.fullname == b.fullname && a.size == b.size

"""
Show definitions for schemas to show them as JSON.
"""
show(io::IO, name::FullName) = serialize(io::IO, name::FullName)
serialize(io::IO, name::FullName, seen=Set{FullName}()) = write(io, "\"$(name.value)\"")

for primitive_type in PRIMITIVE_TYPES
    primitive_json = "\"$primitive_type\""
    classname = Symbol(capitalize(primitive_type), "Schema")
    @eval begin
        serialize(io::IO, ::$classname, seen=Set{FullName}()) = print(io, $primitive_json)
        show(io::IO, x::$classname) = serialize(io::IO, x::$classname)
    end
end

show(io::IO, field::Field) = serialize(io, field)
function serialize(io::IO, field::Field, seen::Set{FullName}=Set{FullName}())
    write(io, "{\"name\":\"$(field.name)\",\"type\":")
    serialize(io, field.schema, seen)
    write(io, "}")
end

show(io::IO, schema::RecordSchema) = serialize(io, schema)
function serialize(io::IO, schema::RecordSchema, seen::Set{FullName}=Set{FullName}())
    if schema.fullname in seen
        serialize(io, schema.fullname, seen)
    else
        push!(seen, schema.fullname)
        write(io, "{\"name\":")
        serialize(io, schema.fullname, seen)
        write(io, ",\"type\":\"record\",\"fields\":[")
        serialize(io, schema.fields[1], seen)
        for field in schema.fields[2:end]
            write(io, ",")
            serialize(io, field, seen)
        end
        write(io, "]}")
    end
end

show(io::IO, schema::EnumSchema) = serialize(io::IO, schema::EnumSchema)
function serialize(io::IO, schema::EnumSchema, seen::Set{FullName}=Set{FullName}())
    if schema.fullname in seen
        serialize(io, schema.fullname, seen)
    else
        push!(seen, schema.fullname)
        write(io, "{\"name\":")
        serialize(io, schema.fullname, seen)
        write(io, ",\"type\":\"enum\",\"symbols\":[")
        write(io, join(["\"$symbol\"" for symbol in schema.symbols], ","))
        write(io, "]}")
    end
end

show(io::IO, schema::ArraySchema) = serialize(io::IO, schema::ArraySchema)
function serialize(io::IO, schema::ArraySchema, seen::Set{FullName}=Set{FullName}())
    write(io, "{\"type\":\"array\",\"items\":")
    serialize(io, schema.items, seen)
    write(io, "}")
end

show(io::IO, schema::MapSchema) = serialize(io::IO, schema::MapSchema)
function serialize(io::IO, schema::MapSchema, seen::Set{FullName}=Set{FullName}())
    write(io, "{\"type\":\"map\",\"values\":")
    serialize(io, schema.values, seen)
    write(io, "}")
end

show(io::IO, schema::UnionSchema) = serialize(io::IO, schema::UnionSchema)
function serialize(io::IO, schema::UnionSchema, seen::Set{FullName}=Set{FullName}())
    write(io, "[")
    serialize(io, schema.schemas[1], seen)
    for item in schema.schemas[2:end]
        write(io, ",")
        serialize(io, item, seen)
    end
    write(io, "]")
end

show(io::IO, schema::FixedSchema) = serialize(io::IO, schema::FixedSchema)
function serialize(io::IO, schema::FixedSchema, seen::Set{FullName}=Set{FullName}())
    if schema.fullname in seen
        serialize(io, schema.fullname, seen)
    else
        push!(seen, schema.fullname)
        write(io, "{\"name\":")
        serialize(io, schema.fullname)
        write(io, ",\"type\":\"fixed\",\"size\":$(schema.size)}")
    end
end

end


text:
using ACTRModels, BenchmarkTools
n = 10
parms = (sa=false, mmp=true, γ=1.3, δ=1.0)
chunks = [Chunk(;s1=j, s2=j) for j in 1:n for k in 1:n]
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Partial Activation n $n")
@btime compute_activation!($actr; s1=1, s2=1);

parms = (sa=true, mmp=true, γ=1.3, δ=1.0)
memory = Declarative(;memory=chunks)
imaginal = Imaginal(;buffer=Chunk(;s1=2,s2=1))
actr = ACTR(;declarative=memory, imaginal=imaginal, parms...)
println("Spreading Activation n $n")
@btime compute_activation!($actr; s1=1, s2=1);

parms = (noise=true,)
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Activation Noise n $n")
@btime compute_activation!($actr);

n = 100
parms = (sa=false, mmp=true, γ=1.3, δ=1.0)
chunks = [Chunk(;s1=j, s2=j) for j in 1:n for k in 1:n]
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Partial Activation n $n")
@btime compute_activation!($actr; s1=1, s2=1);

parms = (sa=true, mmp=true, γ=1.3, δ=1.0)
memory = Declarative(;memory=chunks)
imaginal = Imaginal(;buffer=Chunk(;s1=2,s2=1))
actr = ACTR(;declarative=memory, imaginal=imaginal, parms...)
println("Spreading Activation n $n")
@btime compute_activation!($actr; s1=1, s2=1);

parms = (noise=true,)
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Activation Noise n $n")
@btime compute_activation!($actr);

n = 10
parms = (sa=false, mmp=true, γ=1.3, δ=1.0)
chunks = [Chunk(;s1=j, s2=j) for j in 1:n for k in 1:n]
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Partial Activation n $n")
@btime retrieval_probs($actr; s1=1, s2=1);

parms = (sa=true, mmp=true, γ=1.3, δ=1.0)
memory = Declarative(;memory=chunks)
imaginal = Imaginal(;buffer=Chunk(;s1=2,s2=1))
actr = ACTR(;declarative=memory, imaginal=imaginal, parms...)
println("Spreading Activation n $n")
@btime retrieval_probs($actr; s1=1, s2=1);

n = 100
parms = (sa=false, mmp=true, γ=1.3, δ=1.0)
chunks = [Chunk(;s1=j, s2=j) for j in 1:n for k in 1:n]
memory = Declarative(;memory=chunks)
actr = ACTR(;declarative=memory, parms...)
println("Partial Activation n $n")
@btime retrieval_probs($actr; s1=1, s2=1);

parms = (sa=true, mmp=true, γ=1.3, δ=1.0)
memory = Declarative(;memory=chunks)
imaginal = Imaginal(;buffer=Chunk(;s1=2,s2=1))
actr = ACTR(;declarative=memory, imaginal=imaginal, parms...)
println("Spreading Activation n $n")
@btime retrieval_probs($actr; s1=1, s2=1);

ext: jl


text:
# Autogenerated wrapper script for Mineos_jll for x86_64-w64-mingw32-libgfortran3
export cucss2sac, eigcon, eigen2asc, endi, green, minos_bran, simpledit, syndat

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Mineos")
JLLWrappers.@declare_executable_product(cucss2sac)
JLLWrappers.@declare_executable_product(eigcon)
JLLWrappers.@declare_executable_product(eigen2asc)
JLLWrappers.@declare_executable_product(endi)
JLLWrappers.@declare_executable_product(green)
JLLWrappers.@declare_executable_product(minos_bran)
JLLWrappers.@declare_executable_product(simpledit)
JLLWrappers.@declare_executable_product(syndat)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        cucss2sac,
        "bin\\cucss2sac.exe",
    )

    JLLWrappers.@init_executable_product(
        eigcon,
        "bin\\eigcon.exe",
    )

    JLLWrappers.@init_executable_product(
        eigen2asc,
        "bin\\eigen2asc.exe",
    )

    JLLWrappers.@init_executable_product(
        endi,
        "bin\\endi.exe",
    )

    JLLWrappers.@init_executable_product(
        green,
        "bin\\green.exe",
    )

    JLLWrappers.@init_executable_product(
        minos_bran,
        "bin\\minos_bran.exe",
    )

    JLLWrappers.@init_executable_product(
        simpledit,
        "bin\\simpledit.exe",
    )

    JLLWrappers.@init_executable_product(
        syndat,
        "bin\\syndat.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Oceananigans.Fields: broadcasted_to_abstract_operation

using Base.Broadcast: Broadcasted
using Base: identity

const BroadcastedIdentity = Broadcasted{<:Any, <:Any, typeof(identity), <:Any}

broadcasted_to_abstract_operation(loc, grid, bc::BroadcastedIdentity) =
    interpolate_operation(loc, Tuple(broadcasted_to_abstract_operation(loc, grid, a) for a in bc.args)...)

broadcasted_to_abstract_operation(loc, grid, op::AbstractOperation) = at(loc, op)

function broadcasted_to_abstract_operation(loc, grid, bc::Broadcasted{<:Any, <:Any, <:Any, <:Any})
    abstract_op = bc.f(loc, Tuple(broadcasted_to_abstract_operation(loc, grid, a) for a in bc.args)...)
    return interpolate_operation(loc, abstract_op) # For "stubborn" BinaryOperations
end


text:
let n::Int64=0
	global increment_n() = (n+=1)
	global reset_n() = (n=0)
end

"""
	read_tree(
		nw_file::AbstractString;
		NodeDataType=DEFAULT_NODE_DATATYPE, force_new_labels=false
	)

Read Newick file `nw_file` and create a `Tree{NodeDataType}` object from it.
`NodeDataType` must be a subtype of `TreeNodeData`, and must have a *callable default outer constructor*: the call `NodeDataType()` must exist and return a valid instance of `NodeDataType`.
Use `force_new_labels=true` to force the renaming of all internal nodes.
"""
function read_tree(
	nw_file::AbstractString;
	NodeDataType=DEFAULT_NODE_DATATYPE, force_new_labels=false
)
	tree = node2tree(read_newick(nw_file; NodeDataType); force_new_labels)
	check_tree(tree)
	return tree
end

"""
	parse_newick_string(
		nw::AbstractString;
		NodeDataType=DEFAULT_NODE_DATATYPE, force_new_labels=false
	)

Parse newick string into a tree. See `read_tree` for more informations.
"""
function parse_newick_string(
	nw::AbstractString;
	NodeDataType=DEFAULT_NODE_DATATYPE, force_new_labels=false, strict=true,
)
	reset_n()
	root = TreeNode(NodeDataType())
	parse_newick!(nw, root, NodeDataType)
	root.isroot = true
	tree = node2tree(root; force_new_labels)
	check_tree(tree; strict)
	return tree
end

"""
	read_newick(nw_file::AbstractString)

Read Newick file `nw_file` and create a graph of `TreeNode` objects in the process.
  Return the root of said graph.
  `node2tree` or `read_tree` must be called to obtain a `Tree` object.
"""
function read_newick(nw_file::AbstractString; NodeDataType=DEFAULT_NODE_DATATYPE)
	@assert NodeDataType <: TreeNodeData
	f = open(nw_file)
	nw = readlines(f)
	close(f)
	if length(nw) > 1
		error("File $nw_file has more than one line.")
	elseif length(nw) == 0
		error("File $nw_file is empty")
	end
	nw = nw[1]
	if nw[end] != ';'
		error("File $nw_file does not end with ';'")
	end
	nw = nw[1:end-1]

	reset_n()
	root = parse_newick(nw; NodeDataType)
	return root
end

"""
	parse_newick(nw::AbstractString; NodeDataType=DEFAULT_NODE_DATATYPE)

Parse newick string into a `TreeNode`.
"""
function parse_newick(nw::AbstractString; NodeDataType=DEFAULT_NODE_DATATYPE)
	reset_n()
	root = TreeNode(NodeDataType())
	parse_newick!(nw, root, NodeDataType)
	root.isroot = true # Rooting the tree with outer-most node of the newick string
	root.tau = missing
	return root
end

"""
	parse_newick!(nw::AbstractString, root::TreeNode)

Parse the tree contained in Newick string `nw`, rooting it at `root`.
"""
function parse_newick!(nw::AbstractString, root::TreeNode, NodeDataType)

	# Setting isroot to false. Special case of the root is handled in main calling function
	root.isroot = false
	# Getting label of the node, after last parenthesis
	parts = map(x->String(x), split(nw, ")"))
	lab, tau = nw_parse_name(String(parts[end]))
	if lab == ""
		lab = "NODE_$(increment_n())"
	elseif length(lab) > 1 && lab[1:2] == "[&" # Deal with extended newick annotations
		lab = "NODE_$(increment_n())" * lab
	end

	root.label, root.tau = (lab,tau)

	if length(parts) == 1 # Is a leaf. Subtree is empty
		root.isleaf = true
	else # Has children
		root.isleaf = false
		if parts[1][1] != '('
			println(parts[1][1])
			error("Parenthesis mismatch. This may be caused by spaces in the newick string.")
		else
			parts[1] = parts[1][2:end] # Removing first bracket
		end
		children = join(parts[1:end-1], ")") # String containing children, now delimited with ','
		l_children = nw_parse_children(children) # List of children (array of strings)

		for sc in l_children
			nc = TreeNode(NodeDataType())
			parse_newick!(sc, nc, NodeDataType) # Will set everything right for subtree corresponding to nc
			nc.anc = root
			push!(root.child, nc)
		end
	end

end

"""
	nw_parse_children(s::AbstractString)

Idea from http://stackoverflow.com/a/26809037
Split a string of children in newick format to an array of strings.

## Example
`"A,(B,C),D"` --> `["A","(B,C)","D"]`
"""
function nw_parse_children(s::AbstractString)
	parcount = 0
	annotation = false # For reading extended newick grammar, with [&key=val] after label
	l_children = []
	current = ""
	cstart = 1
	cend = 1
	for (i,c) in enumerate("$(s),")
		if c == '['
			annotation = true
		elseif c == ']'
			annotation = false
		elseif c == ',' && parcount == 0 && annotation == false
			cend = i-1
			push!(l_children, s[cstart:cend])
			cstart = i+1
		else
			if c == '('
				parcount +=1
			elseif c == ')'
				parcount -=1
			end
		end
	end
	return l_children
end

"""
	nw_parse_name(s::AbstractString)

Parse Newick string of child into name and time to ancestor. Default value for missing time is `missing`.
"""
function nw_parse_name(s::AbstractString)
	temp = split(s, ":")
	if occursin(':', s) # Node has a time
		if length(temp) == 2 # Node also has a name, return both
			tau = (tau = tryparse(Float64,temp[2]); isnothing(tau) ? missing : tau) # Dealing with unparsable times
			return string(temp[1]), tau
		else # Return empty name
			tau = (tau = tryparse(Float64,temp[1]); isnothing(tau) ? missing : tau) # Dealing with unparsable times
			return "", tau
		end
	else # Node does not have a time, return string as name
		return s, missing
	end
end

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools

max_forks_repo_name: PierreBarrat/TreeTools


text:
import ..BalanceLaws: projection

# Zero-out vertical momentum tendencies based on compressibility
function projection(pv::Momentum, atmos::AtmosModel, td::TendencyDef, args, x)
    return projection(pv, compressibility_model(atmos), td, args, x)
end

# Zero-out vertical momentum fluxes for Anelastic1D:
function projection(
    pv::Momentum,
    ::Anelastic1D,
    ::TendencyDef{Flux{O}},
    args,
    x,
) where {O}
    return x .* SArray{Tuple{3, 3}}(1, 1, 1, 1, 1, 1, 0, 0, 0)
end

# Zero-out vertical momentum sources for Anelastic1D:
function projection(::Momentum, ::Anelastic1D, ::TendencyDef{Source}, args, x)
    return x .* SVector(1, 1, 0)
end


text:
using Sockets

const portnum = 12345

function canopen()
    try
        server = listen(portnum)
        println("This is the only instance.")
        sleep(20)
    catch y
        if findfirst("EADDRINUSE", string(y)) != nothing
            println("There is already an instance running.")
        end
    end
end

canopen()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JuLIP: Atoms, AbstractCalculator
export EnergyMixing, prepare_qmmm!

import JuLIP: energy, forces

# ------------------------------------------------------------------

function get_domains(at::Atoms)
   Iqm = get_data(at, "qmmm_e_Iqm")
   # Iqmbuf = get_data(at, "qmmm_e_Iqmbuf")
   Imm = setdiff(1:length(at), Iqm)
   return Iqm, Imm, nothing
end

struct EnergyMixing{T} <: AbstractCalculator
   atqm::Atoms{T}
   QM::AbstractCalculator
   MM::AbstractCalculator
end

function energy(calc::EnergyMixing{T}, at::Atoms{T}) where {T}
   Iqm, Imm, Iqmbuf = get_domains(at)
   # set_positions!(calc.atqm, positions(at)[Iqmbuf])
   Eqm = energy(calc.QM, at; domain = Iqm)
   Emm = energy(calc.MM, at; domain = Imm)
   return Eqm + Emm
end

function forces(calc::EnergyMixing{T}, at::Atoms{T}) where {T}
   Iqm, Imm, Iqmbuf = get_domains(at)
   # set_positions!(calc.atqm, positions(at)[Iqmbuf])
   Fqm = forces(calc.QM, at; domain = Iqm)
   Fmm = forces(calc.MM, at; domain = Imm)
   return Fqm + Fmm
end

# ------------------------------------------------------------------

struct ForceMixing{T}
   atqm::Atoms{T}
   QM::AbstractCalculator
   MM::AbstractCalculator
end

function forces(calc::ForceMixing{T}, at::Atoms{T}) where {T}
   Iqm, Imm, Iqmbuf = get_domains(at)
   Fqm = forces(calc.QM, at; domain = Iqm)
   Fmm = forces(calc.MM, at; domain = Imm)
   return Fqm + Fmm
end

# ------------------------------------------------------------------

prepare_qmmm!(at; kwargs...) = prepare_qmmm!(at, EnergyMixing; kwargs...)

function prepare_qmmm!(at::Atoms, ::Type{EnergyMixing};
                       Vqm = nothing, Vmm = nothing,
                       Iqm = nothing)
   # perform some sanity checks
   Iqm = collect(Iqm)
   @assert Iqm isa Vector{<: Integer}
   @assert Vqm isa AbstractCalculator
   @assert Vmm isa AbstractCalculator
   @assert sort(unique(Iqm)) == sort(Iqm)
   @assert minimum(Iqm) > 0
   @assert maximum(Iqm) <= length(at)

   # create an atoms object for the QM calculator
   # => TODO: add buffer, carve out region, etc...

   # store the QM region in `at`
   at["qmmm_e_Iqm"] = Iqm
   # store the new calculator
   set_calculator!(at, EnergyMixing(at, Vqm, Vmm))

   return at
end


text:
"A `DomainPoint` is a point which is an element of a domain by construction."
abstract type DomainPoint{T} end

in(p::DomainPoint, d::Domain) = domain(p) == d || in(point(p), d)

## Points on a sphere

"A point on the unit sphere."
abstract type SpherePoint{T} <: DomainPoint{T} end

domain(p::SpherePoint{T}) where {T<:StaticTypes} = UnitSphere{T}()
domain(p::SpherePoint{T}) where {T<:AbstractVector} = UnitSphere{T}(length(point(p)))

"A point on the unit sphere represented by a standard Euclidean vector."
struct EuclideanSpherePoint{T} <: SpherePoint{T}
    x   ::  T
end
point(p::EuclideanSpherePoint) = p.x

"A point on the unit sphere represented in spherical coordinates."
struct SphericalCoordinate{T} <: SpherePoint{SVector{3,T}}
    θ   ::  T   # inclination or polar angle
    ϕ   ::  T   # azimuthal angle
end

point(p::SphericalCoordinate) = SVector(sin(p.θ)*cos(p.ϕ), sin(p.θ)*sin(p.ϕ), cos(p.θ))


text:
using Plots
using Random

include_implicit_dynamics()
include_ddp()
include(joinpath(pwd(), "examples/implicit_dynamics/models/double_pendulum/model.jl"))

# problem setup
T = 101
h = 0.05

# build implicit dynamics

# impact model
data = dynamics_data(model, h, T,
        r_func, rz_func, rθ_func, rz_array, rθ_array;
        idx_ineq = idx_ineq,
        z_subset_init = 0.1 * ones(4),
        # diff_idx = -1, # 630, 67.62
        # diff_idx = 6, # 582, 133.3
        # diff_idx = 5, # 206, 9.6
        diff_idx = 4, # 128.6, 9.2
        # diff_idx = 3, # failure
        # diff_idx = 2, # failure
        opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 10.0,
						diff_sol = false))

model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model.dim.q, model.dim.u, 0, data)

# no impact model
data = dynamics_data(model_no_impact, h, T,
        r_no_impact_func, rz_no_impact_func, rθ_no_impact_func,
		rz_no_impact_array, rθ_no_impact_array;
        idx_ineq = collect(1:0),
		z_subset_init = zeros(0),
        diff_idx = -1,
        opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 0.1,
						κ_init = 0.1,
						diff_sol = false))

model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model_no_impact.dim.q, model_no_impact.dim.u, 0, data)

q0 = [0.0; 0.0]
q1 = [0.0; 0.0]
qT = [π; 0.0]
q_ref = [π; 0.0]

x1 = [q0; q1]
xT = [qT; qT]

n = model_implicit.n
nq = model_implicit.dynamics.m.dim.q
m = model_implicit.m

Random.seed!(1) # reset random
ū = [1.0e-3 * randn(model_implicit.m) for t = 1:T-1]
w = [zeros(model_implicit.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model_implicit, x1, ū, w, h, T)

# Objective
V = 0.1 * Diagonal(ones(nq))
_Q = [V -V; -V V] ./ h^2.0
Q = [t < T ? _Q : _Q for t = 1:T]
q = [-2.0 * Q[t] * zeros(n) for t = 1:T]
R = [Diagonal(1.0 * ones(m)) for t = 1:T-1]
r = [zeros(m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = [-1.0]
uu = [1.0]
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		# c .= [ul - u; u - uu]
	else
		c .= x - cons.con[T].info[:xT]
	end
end

prob = problem_data(model_implicit, obj, con_set, copy(x̄), copy(ū), w, h, T,
    analytical_dynamics_derivatives = true)

# Solve
con_tol = 0.005
Random.seed!(1) # reset random
iters = []
con_check = []
# for i = 1:5
ū = [1.0e-3 * randn(model_implicit.m) for t = 1:T-1]

prob = problem_data(model_implicit, obj, con_set, copy(x̄), copy(ū), w, h, T,
	analytical_dynamics_derivatives = true)

@time stats = constrained_ddp_solve!(prob,
    linesearch = :armijo,
    α_min = 1.0e-5,
    obj_tol = 1.0e-3,
    grad_tol = 1.0e-3,
	max_iter = 1000,
    max_al_iter = 15,
	ρ_init = 1.0,
    ρ_scale = 10.0,
	con_tol = con_tol)

# push!(iters, ilqr_iterations(stats))
# push!(con_check, prob.s_data.c_max <= con_tol)
# end

# @show all(con_check)
# @show mean(iters)
# @show std(iters)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)
λ = [data.z_cache[t][model.dim.q .+ (1:nc)] for t = 1:T-1]

# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/acrobot.jld2") x u λ
# @load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/acrobot.jld2") x u λ

@save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/acrobot_no_limits.jld2") x u
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/acrobot_no_limits.jld2") x u

q̄ = state_to_configuration(x)
q2l = -0.5 * π * ones(length(q̄))
q2u = 0.5 * π * ones(length(q̄))
t = range(0, stop = h * (length(q̄) - 2), length = length(q̄)-1)
plt = plot();
plt = plot!(t, q2l[2:end], color = :black, width = 2.0, label = "q2 limit lower")
plt = plot!(t, q2u[2:end], color = :black, width = 2.0, label = "q2 limit upper")
plt = plot!(t, hcat(q̄[2:end]...)', width = 2.0,
	color = [:magenta :orange],
	labels = ["q1" "q2"],
	legend = :topleft,
	xlabel = "time (s)",
	ylabel = "configuration",
	title = "acrobot (w/o joint limits)")
	# title = "acrobot (w/ joint limits)")

λ = [data.z_cache[t][model.dim.q .+ (1:nc)] for t = 1:T-1]
λ_max = maximum(hcat(λ...))
plot!(t, ((1.0 / λ_max) * hcat(λ..., λ[end]) .- 0.0 * π)',
    width = 2.0,
    color = :cyan, linetype = :steppost)

λ_normalized = ((1.0 / λ_max) * hcat(λ..., λ[end]) .- 0.0 * π)

# using PGFPlots
# const PGF = PGFPlots
#
# plt_q1 = PGF.Plots.Linear(t, hcat(q̄[2:end]...)[1, :],
# 	mark="none",style="color=magenta, line width = 2pt",legendentry="q1")
#
# plt_q2 = PGF.Plots.Linear(t, hcat(q̄[2:end]...)[2, :],
# 	mark="none",style="color=orange, line width = 2pt",legendentry="q2")
#
# plt_ql = PGF.Plots.Linear(t, q2l[2:end],
# 	mark="none",style="color=black, line width = 2pt, dashed")#,legendentry="qL")
#
# plt_qu = PGF.Plots.Linear(t, q2u[2:end],
# 	mark="none",style="color=black, line width = 2pt, dashed")#,legendentry="qU")
#
# plt_l1 = PGF.Plots.Linear(t, λ_normalized[1, :],
# 	mark="none",style="const plot, color=green, line width = 2pt",legendentry="l1")
#
# plt_l2 = PGF.Plots.Linear(t, λ_normalized[2, :],
# 	mark="none",style="const plot, color=cyan, line width = 2pt",legendentry="l2")
#
# a = Axis([plt_q1, plt_q2, plt_l1, plt_l2, plt_ql, plt_qu],
#     axisEqualImage=false,
#     hideAxis=false,
# 	ylabel="configuration",
# 	xlabel="time (s)",
# 	# xlim=(0.0, 5.0),
# 	legendStyle="{at={(0.01,0.99)},anchor=north west}")
#
# PGF.save("/home/taylor/Research/implicit_dynamics_manuscript/figures/acrobot_configuration.tikz", a, include_preamble=false)

# show(plt)
# savefig(plt, "/home/taylor/Research/implicit_dynamics_manuscript/figures/acrobot_joint_limits.png")
# savefig(plt, "/home/taylor/Research/implicit_dynamics_manuscript/figures/acrobot_no_joint_limits.png")

# plot(hcat(ū..., ū[end])', linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/models/double_pendulum/visuals.jl"))
vis = Visualizer()
render(vis)
open(vis)
default_background!(vis)
settransform!(vis["/Cameras/default"],
        compose(Translation(0.0, -95.0, -1.0), LinearMap(RotY(0.0 * π) * RotZ(-π / 2.0))))
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 30)
setvisible!(vis["/Grid"], false)

visualize_elbow!(vis, model,
    # q̄,
    [[q̄[1] for i = 1:50]..., q̄..., [q̄[end] for i = 1:50]...],
    # limit_color = RGBA(0.0, 1.0, 0.0, 1.0),
    limit_color = RGBA(1.0, 0.0, 0.0, 1.0),
    Δt = h)
#
# # ghost
# limit_color = [1.0, 0.0, 0.0]
# # limit_color = [0.0, 1.0, 0.0]
#
# t = 1
# id = t
# tl = 0.05
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 10
# id = t
# tl = 0.15
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 20
# tl = 0.25
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 30
# id = t
# tl = 0.35
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 40
# id = t
# tl = 0.45
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 50
# id = t
# tl = 0.55
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 60
# id = t
# tl = 0.65
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 70
# id = t
# tl = 0.75
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 80
# id = t
# tl = 0.85
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 90
# id = t
# tl = 0.95
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# t = 101
# id = t
# tl = 1.0
# _create_acrobot!(vis, model;
#         tl = tl,
#         limit_color = RGBA(limit_color..., tl),
#         i = id)
# _set_acrobot!(vis, model, x[t], i = id)
#
# # line_mat = LineBasicMaterial(color=color=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0), linewidth=10.0)
# line_mat = LineBasicMaterial(color=color=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0), linewidth=10.0)
#
# points = Vector{Point{3,Float64}}()
# for (i, xt) in enumerate(x̄)
#     k = kinematics(model, xt)
# 	push!(points, Point(k[1], 0.0, k[2]))
#
#     setobject!(vis["ee_vertex_$i"], Sphere(Point3f0(0),
#         convert(Float32, 0.001)),
#         MeshPhongMaterial(color = RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0)))
#         settransform!(vis["ee_vertex_$i"], Translation(points[i]))
# end
# setobject!(vis[:ee_traj], MeshCat.Line(points, line_mat))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using PyPlot
using LinearAlgebra
using SparseArrays
using Arpack
#using BenchmarkTools
#using Formatting
using FileIO
using DelimitedFiles
using FEMSparse
using Flux:batch
using IterativeSolvers

t0 = time()
#readfiles
println("opening file... \n")
a = open("/home/billykon/programing/2d-elasticity-quadrilateral/plaka15x1.unv", "r")
b = open("/home/billykon/programing/2d-elasticity-quadrilateral/plaka15x1.dat", "r")
include("readfile.jl") #returns array of: nodes, edge, elem, bc
#include("mesh.jl")
include("elastic_functions.jl") #defines basic functions

#start building global array
n = length(elem[:, 1])
dofs = 2 * length(nodes[:, 1])
v = .3
Y = 200e9
global E = Ematrix(Y, v, "stress")

println("calculating stiffness/mass matrix...\n")

K = SparseMatrixCOO()
for i = 1:n
    i1 = elem[i, :]'
    xs = nodes[i1, 1]
    ys = nodes[i1, 2]
    i1 = [i1 (i1 .+ Int64(0.5 * dofs))][:]
    kre = ke(E, xs, ys)

    add!(K,i1,i1,kre)

    if i % 1000 ==0
        println(i)
    end
end
K = SparseMatrixCSC(K)
uv = [0.5 0.5]

#start boundary conditions
println("calculating right-hand-side...\n")
F = zeros(dofs, 1)

#---------------------------------------------------
P=1e6
global f(x) = [[0 P*(x[2]-x[2]^2)*6]', [0 0]']
global f(x) = [[0 P]', [0 0]']
#---------------------------------------------------
println("finding IDs...")
(nbc, ni, ix) = IDs(bc) #nbc has number of conditions
#ix has edges nodes ids for each condition
global bc, edge, elem, nodes, i, F
i = 1

print("main loop to assign rhs...")
while true
    global i, ix, elem
    if i > nbc
        break
    end
    rows = [j for k=1:ni[i] for j=1:size(elem)[1] if issubset(ix[i][k, 1], elem[j, :]) &&
               issubset(ix[i][k, 2], elem[j, :]) &&
               issubset(ix[i][k, 3], elem[j, :])] #contains superset id element
    xedge = nodes[elem[rows, :], :][:, :, 1]
    yedge = nodes[elem[rows, :], :][:, :, 2] #contains x,y of edges

    superelem = elem[rows, :]
    mapx = [sum(dims = 2, superelem[k, :] .== ix[i][k, :]')  for k = 1:length(superelem[:, 1])]#zero the non edge indexes of element
    for k = 1:length(mapx)
        g(x) = f(x)[i]
        a = BC(g, xedge[k, :]', yedge[k, :]', mapx[k])
        id = [ix[i][k, :]' (ix[i][k, :]' .+ Int64(0.5 * dofs))]
        id = sort(id, dims = 2)

        idlocal = [mapx[k]' mapx[k]'] #has the global degrees of freedoms
        idlocal = [idlocal[i] .== 1 for i = 1:18]
        F[id] .= F[id] .+ a[idlocal]'
    end
    i = i + 1
end
#gib traction nodes where (u,v) are known
#reminder, ix contains boundary condition nodes INDEXES
println("solving for U....\n")
freedofs = zeros(dofs, 1)
fixedofs = ix[2][:]'#take first's b.c. nodes

fixedofs = [fixedofs Int.(fixedofs .+0.5*dofs)]
freedofs[fixedofs] .= 1
freedofs = [i for i = 1:dofs if freedofs[i] == 0]

GC.gc() #removes garbage RAM
println("initializing U with 0....")
U = zeros(dofs, 1)
a =1
b = 1

println("solving the Kfree for UFree system...")
varinfo()
Ufree = K[freedofs, freedofs] \ F[freedofs]
#Ufree = minres(K[freedofs,freedofs],F[freedofs],verbose=true)
println("assigning Ufree to U")
U[freedofs] .= Ufree
Ufree = 1
freedofs = 1
fixedofs = 1
K = 1
Kfree = 1

include("stress.jl")

println("saving nodes & stresses")
writedlm("/home/billykon/programing/2d-elasticity-quadrilateral/data.csv",["x" "y" "z" "sx" "sy" "txy" ; nodes fstress' ./ 1e6],',')
println("saving Us....")
writedlm("/home/billykon/programing/2d-elasticity-quadrilateral/U.csv",U,',')

println("max U ",maximum(U))
println(time()-t0,'s')

#post processing of bending
#below analytical solution is compared to the found one
amplify = 50
u = amplify .*U
nodex = nodes[:,1]
nodey = nodes[:,2]
nodedx = u[1:Int(0.5*dofs)].+nodex
nodedy = u[Int(0.5*dofs)+1:Int(dofs)].+nodey

figure(1)
scatter(nodedx,nodedy,marker="x")
axis("equal")

#analytical bending, I = 1/12, D=1, L=15
uexact = -(2*P/Y).*(nodey.-0.5).*(3. .*nodex.*(30. .-nodex).+2*v.*nodey.*(nodey.-1.))
vexact = (2*P/Y).*(nodex.^2 .*(45. .-nodex)+3*v*(15. .-nodex).*(nodey.-0.5).^2 +(4+5*v).*nodex./4.)
nodedxact = nodex .+ uexact.*amplify
nodedyact = nodey .+ vexact.*amplify

scatter(nodedxact,nodedyact,marker="o",edgecolors="red",facecolor="none")
axis("equal")

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using Conda
using PyCall

begin # See learn_PyCall.jl
    path = Conda.bin_dir(dirname(PyCall.pyprogramname))
    if isnothing(findfirst(path, ENV["PATH"])) && Sys.iswindows()
        ENV["PATH"] = path * ";" * ENV["PATH"]
    end
end

# Must load in Python before loading in Julia

py"""
import wx
"""
#pygui_start(:wx)
wx = pyimport("wx")

# https://wiki.wxpython.org/Getting%20Started
function learn_hello()
    app = wx.App(false)
    frame = wx.Frame(nothing, wx.ID_ANY, "Hello World")
    frame.Show(true)
    app.MainLoop()

end

# https://wiki.wxpython.org/Getting%20Started#Building_a_simple_text_editor
function learn_building_a_simple_text_editor()
    # py"""
    # import wx
    # class MyFrame(wx.Frame):
    #     def __init__(self, parent, title):
    #         wx.Frame.__init__(self, parent, title=title, size=(200,100))
    #         self.control = wx.TextCtrl(self, style=wx.TE_MULTILINE)
    #         self.Show(True)
    # """
    # app = wx.App(false)
    # frame = py"MyFrame"(nothing, "Small editor")
    # app.MainLoop()

    # Use my own package
    lwx = pyimport("learnall.wx")
    app = wx.App(false)
    frame = lwx.MyFrame(nothing, "Small editor")
    #app.MainLoop()

end

#learn_hello()
#learn_building_a_simple_text_editor()

nothing


text:
# This code was due to Jeffrey Sarnoff. See
# https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/372

using IntervalArithmetic
import SpecialFunctions
import ArbNumerics
using ArbNumerics: ArbFloat, ArbReal, setinterval, lowerbound, upperbound

# number of bits that suffices to alter the least significant digit
#   use: prevfloat(x, DigitBits), nextfloat(x, DigitBits)
const DigitBits = ceil(Int, log2(10))

# Interval{T} --> ArbReal, thin enclosure
function arbreal(x::Interval{T}) where T
     setprecision(ArbFloat, T)
     lo = prevfloat(ArbFloat(x.lo), DigitBits) # assure least significant
     hi = nextfloat(ArbFloat(x.hi), DigitBits) #   decimal digit changes
     return setinterval(lo, hi)
end

# T, ArbReal --> Interval{T}
#    P is the working precision from ArbNumerics
function tinterval(::Type{T}, x::ArbReal{P}) where {T, P}
     lo, hi = lowerbound(x), upperbound(x)
     return Interval{T}(T(lo), T(hi))
end

function specialfunction(fn, a::Interval{T}) where T
     x = fn(arbreal(a))
     return tinterval(T, x)
end

function specialfunction(fn, a::Interval{T}, b::Interval{T}) where T
     x = fn(arbreal(a), arbreal(b))
     return tinterval(T, x)
end

function specialfunction(fn, a::Interval{T}, b::Interval{T},
                             c::Interval{T}) where T
     x = fn(arbreal(a), arbreal(b), arbreal(c))
     return tinterval(T, x)
end

function specialfunction(fn, a::Interval{T}, b::Interval{T},
                             c::Interval{T}, d::Interval{T}) where T
     x = fn(arbreal(a), arbreal(b), arbreal(c), arbreal(d))
     return tinterval(T, x)
end

Base.setprecision(::Type{ArbFloat}, ::Type{Float64}) =
    setprecision(ArbFloat, 53)
Base.setprecision(::Type{ArbFloat}, ::Type{Float32}) =
    setprecision(ArbFloat, 24)
Base.setprecision(::Type{ArbFloat}, ::Type{BigFloat}) =
    setprecision(ArbFloat, precision(BigFloat))

# Don't import these functions
# functions of 1 arg from Base.Math (for Interval{BigFloat})
#for F in (:acos, :acosh, :acot, :acoth, :acsc, :acsch, :asec,
#          :asech, :asin, :asinh, :atan, :atanh, :cbrt, :cos,
#          :cosc, :cosh, :cospi, :cot, :coth, :csc, :csch,
#          :exp, :expm1, :log, :log10, :log1p, :log2, :sec,
#          :sech, :sin, :sinc, :sinh, :sinpi, :sqrt, :tan, :tanh)
#    @eval Base.$F(x::A) where {T, A<:Interval{BigFloat}} =
#          specialfunction($F,x)
#end

# functions of 2 arg from Base.Math (for Interval{BigFloat})
for F in (:^, :atan, :hypot)
  @eval begin
    Base.$F(x::A, y::A) where {T, A<:Interval{BigFloat}} =
        specialfunction($F, x, y)
    Base.$F(x::A, y::Real) where {T, A<:Interval{BigFloat}} =
        specialfunction($F, x, A(y))
    Base.$F(x::Real, y::A) where {T, A<:Interval{BigFloat}} =
        specialfunction($F, A(x), y)
  end
end

# functions of 1 arg exported by SpecialFunctions and ArbNumerics
for F in (:airyai, :airyaiprime, :airybi, :airybiprime,
          :besselj0, :besselj1, :bessely0, :bessely1,
          #:erf, :erfc, :erfi,
          :gamma, :lgamma, :digamma, :trigamma,
          :eta, :zeta)
    @eval SpecialFunctions.$F(x::A) where {T, A<:Interval{T}} =
          specialfunction($F,x)
end

# functions of 1 arg exported by both with distinct names
for (SF, AF) in ((:cosint, :ci), (:sinint, :si),
                 (:ellipk, :elliptic_k), (:ellipe, :elliptic_e))
    @eval SpecialFunctions.$SF(x::A) where {T, A<:Interval{T}} =
          specialfunction($AF,x)
end

# functions of 1 arg exported by ArbNumerics only
for F in (:chi, :shi, :ei, :xi,  :lambertw)
    @eval ArbNumerics.$F(x::A) where {T, A<:Interval{T}} =
          specialfunction($F,x)
end

# functions of a 2 args exported by ArbNumerics only
for F in (:elliptic_pi, :elliptic_pi2, # the 2 arg versions
          :polylog,
          :hypergeometric_0F1, :regular_hypergeometric_0F1)
    @eval ArbNumerics.$F(x::A, y::A) where {T, A<:Interval{T}} =
          specialfunction($F, x, y)
end

# functions of a 3 args exported by ArbNumerics only
for F in (:elliptic_pi, :elliptic_pi2, # the 3 arg versions
          :hypergeometric_1F1, :regular_hypergeometric_1F1)
    @eval ArbNumerics.$F(a::A, b::A, z::A) where {T, A<:Interval{T}} =
          specialfunction($F, a, b, z)
end

# functions of a 4 args exported by ArbNumerics only
for F in (:hypergeometric_2F1, :regular_hypergeometric_2F1)
    @eval ArbNumerics.$F(a::A, b::A, c::A, z::A) where
                                   {T, A<:Interval{T}} =
          specialfunction($F, a, b, c, z)
end

function __init__()
    # provide padding for prevfloat, nextfloat in arbreal()
    extrabits = ArbNumerics.ExtraBits[] + DigitBits
    ArbNumerics.setextrabits(extrabits)
    return nothing
end


text:
using DSGE
using Test

path = dirname(@__FILE__)

### Model
model = Model1002("ss2")

### Parameters

@testset "Test parameter bounds checking" begin
    for θ in model.parameters
        @test isa(θ, AbstractParameter)
        if !θ.fixed
            (left, right) = θ.valuebounds
            @test left < θ.value < right
        end
    end
end

### Model indices

@testset "Checking model indices" begin
    # Endogenous states
    endo = model.endogenous_states
    @test length(endo) == 68
    @test endo[:Ez_t] == 60

    # Exogenous shocks
    exo = model.exogenous_shocks
    @test length(exo) == 24
    @test exo[:corepce_sh] == 16

    # Expectation shocks
    ex = model.expected_shocks
    @test length(ex) == 13
    @test ex[:EL_f_sh] == 12

    # Equations
    eq = model.equilibrium_conditions
    @test length(eq) == 68
    @test eq[:eq_Ez] == 59

    # Additional states
    endo_new = model.endogenous_states_augmented
    @test length(endo_new) == 16
    @test endo_new[:y_t1] == 69

    # Observables
    obs = model.observables
    @test length(obs) == 19
    @test obs[:obs_tfp] == 12
end

### Equilibrium conditions
Γ0, Γ1, C, Ψ, Π = eqcond(model)

# Transition and measurement equations
TTT, RRR, CCC = solve(model)
meas = measurement(model, TTT, RRR, CCC)

### Pseudo-measurement equation
pseudo_meas = pseudo_measurement(model, TTT, RRR, CCC)

@testset "Check eqcond and state-space dimensions" begin
    # Matrices are of expected dimensions
    @test size(Γ0) == (68, 68)
    @test size(Γ1) == (68, 68)
    @test size(C) == (68,)
    @test size(Ψ) == (68, 24)
    @test size(Π) == (68, 13)

    @test size(meas[:ZZ]) == (19,84)
    @test size(meas[:DD]) == (19,)
    @test size(meas[:QQ]) == (24,24)
    @test size(meas[:EE]) == (19,19)

    @test size(TTT) == (84,84)
    @test size(RRR) == (84,24)
    @test size(CCC) == (84,)
end


text:
module Perf
# compare performance with plain assertion
using BenchmarkTools
using ArgCheck

truthy(x) = true

function fallback_argcheck(x)
    @argcheck x
end
function comparison_argcheck(x)
    @argcheck x == x
end
function call_argcheck(x)
    @argcheck truthy(x)
end
function fallback_assert(x)
    @assert x
end
function comparison_assert(x)
    @assert x == x
end
function call_assert(x)
    @assert truthy(x)
end

benchmarks =[
    (fallback_assert, fallback_argcheck, true),
    (call_assert, call_argcheck, 42),
    (comparison_assert, comparison_argcheck, 42),
    ]

for (f_argcheck, f_assert, arg) in benchmarks
    println(f_argcheck)
    @btime ($f_argcheck)($arg)
    println(f_assert)
    @btime ($f_assert)($arg)
end

end#module


text:
#! /usr/bin/env julia

io = open("output.txt","w")

a = "Hello, World!"

serialize(io,a)

close(io)

f = open("output.txt","r")

println(deserialize(f))

close(f)

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
using Documenter, NeuralQuantum

makedocs(
    modules   = [NeuralQuantum],
    format    = Documenter.HTML(),
    sitename  = "NeuralQuantum.jl",
    authors   = "Filippo Vicentini",
    pages     = [
            "Home"          => "index.md",
            "Manual"        => Any[
                "Basics"        => "basics.md",
                "Problems"      => "problems.md",
                "Algorithms"    => "algorithms.md",
                "Networks"      => "networks.md",
                "Optimizers"    => "optimizers.md"
            ],
            "Internals"     => Any[
                "States"        => "states.md",
            ]
    ]
)

deploydocs(
    repo   = "github.com/PhilipVinc/NeuralQuantum.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
function sub_project_by_sample(project, sample_names)
	columns = append!([:PFAM_Model], reverse(map(symbol, sample_names)))
	new_read_counts_table = project["read_counts_table"][columns]
	new_rpoB_equiv_table = project["rpoB_equiv_table"][columns]
	new_project_table = DataFrame()
	for sample in sample_names
		new_project_table = vcat(new_project_table, project["project_table"][project["project_table"][:SampleID] .== sample,:])
	end
	new_seq_info_table = DataFrame()
	if size(project["seq_info_table"]) != (0,0)
		for sample in sample_names
			new_seq_info_table = vcat(new_seq_info_table, project["seq_info_table"][project["seq_info_table"][:Sample] .== sample,:])
		end
	else
		new_seq_info_table = project["seq_info_table"]
	end

	return Dict("project_table"=>new_project_table,"read_counts_table"=>new_read_counts_table,"hmm_database_info"=>project["hmm_database_info"],"rpoB_equiv_table"=>new_rpoB_equiv_table, "seq_info_table"=>new_seq_info_table)

end

function sub_project_by_model(project, model_names)
	new_read_counts_table = DataFrame()
	new_rpoB_equiv_table = DataFrame()
	new_seq_info_table = DataFrame()
	for model in model_names
		new_read_counts_table = vcat(new_read_counts_table, project["read_counts_table"][project["read_counts_table"][:PFAM_Model] .== model,:])
		new_rpoB_equiv_table = vcat(new_rpoB_equiv_table, project["rpoB_equiv_table"][project["rpoB_equiv_table"][:PFAM_Model] .== model,:])
		if size(project["seq_info_table"]) != (0,0)
			new_seq_info_table = vcat(new_seq_info_table, project["seq_info_table"][project["seq_info_table"][:model] .== model,:])
		else
			new_seq_info_table = DataFrame()
		end
	end

	return Dict("project_table"=>project["project_table"],"read_counts_table"=>new_read_counts_table,"hmm_database_info"=>project["hmm_database_info"],"rpoB_equiv_table"=>new_rpoB_equiv_table, "seq_info_table"=>new_seq_info_table)

end

function sub_project_by_metadata(project, parameter, operator, value)
	expression = Expr(:call,operator,project["project_table"][parameter],value)
	new_sample_IDs = project["project_table"][eval(expression),:][:SampleID]
	return sub_project_by_sample(project, new_sample_IDs)
end

function merge_projects(project1, project2)
	sample_intersect = intersect(project1["project_table"][:SampleID],project2["project_table"][:SampleID])
	if length(sample_intersect) > 0
		error("Both projects share the following sample names and therefore cannot be merged:\n"*sample_intersect)
	end

	if project1["hmm_database_info"] != project2["hmm_database_info"]
		error("It looks like these projects were generated using different HMM databases. Merging them probably isn't the best idea.")
	end

	new_read_counts_table = join(project1["read_counts_table"], project2["read_counts_table"], on = :PFAM_Model, kind = :outer)
	new_rpoB_equiv_table = join(project1["rpoB_equiv_table"], project2["rpoB_equiv_table"], on = :PFAM_Model, kind = :outer)
	new_seq_info_table = vcat(project1["seq_info_table"],project2["seq_info_table"])
	new_project_table = vcat(project1["project_table"],project2["project_table"])
	new_hmm_database_info = project1["hmm_database_info"]

	return Dict("project_table"=>new_project_table,"read_counts_table"=>new_read_counts_table,"hmm_database_info"=>new_hmm_database_info,"rpoB_equiv_table"=>new_rpoB_equiv_table, "seq_info_table"=>new_seq_info_table)

end

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
# TODO: define global verbosity throughout the package

# function for normal generation

"""
    sample_normal(;prompt="I hate the customs", tokenizer=nothing, model=nothing, max_length=100, method="top_k", k=50, t=1.2, p=0.5, add_eos_start=true)

Function to generate normal Sentences with `model` and `tokenizer` provided. In case not provided, function itself create instance of GPT2-small tokenizer and LM Head Model. The sentences are started with the provided `prompt`, and generated till token length reaches `max_length`.

Two sampling methods of generation are provided with this function:
1) method='top_k'
2) method='nucleus'

Any of these methods can be used provided with either k or p.

"""
function sample_normal(;prompt="I hate the customs", tokenizer=nothing, model=nothing, max_length=100, method="top_k", k=50, t=1.2, p=0.5, add_eos_start=true)
    global device
    if tokenizer==nothing
        tokenizer = load_pretrained_tokenizer(GPT2)
        model = hgf"gpt2:lmheadmodel"
    end
    if add_eos_start==true
        input_ = [tokenizer.eos_token_id; tokenizer(prompt)]
    else
        input_ = tokenizer(prompt)
    end
    model = model |> device
    @showprogress "Generating..." for i in 1:max_length
        input_ids = reshape(input_[:], :, 1) |> device
        outputs = model(input_ids; output_attentions=false,
                        output_hidden_states=false,
                        use_cache=false)
        logits = outputs.logits[:, end, 1]
        probs = temp_softmax(logits; t=1.2)  |> cpu
        if method=="top_k"
            new_token = top_k_sample(probs; k=k)[1]
        elseif method == "nucleus"
            new_token = nucleus_sample(probs; p=p)[1]
        end
        push!(input_, new_token)
    end
    text = detokenize(tokenizer, input_)
    if add_eos_start==true
        text = text[length(tokenizer.eos_token)+1:end]
    end
    print("Unpurturbed generated text: ", text, "\n")
    return text
end

# function for perturbed generation
"""
    function sample_pplm(pplm_args;tokenizer=nothing, model=nothing, prompt="I hate the customs", sample_method="top_k", add_eos_start=true)

Function for PPLM model based generation. Generate perturbed sentence using `pplm_args`, tokenizer and model (GPT2, in case not provided), starting with `prompt`. In this function the generation is based on the arguments/parameters provided in `pplm_args`, which is an instance of `pplm` struct.

"""
function sample_pplm(pplm_args;tokenizer=nothing, model=nothing, prompt="I hate the customs", sample_method="top_k", add_eos_start=true)
    if tokenizer == nothing
        if pplm_args.embed_size==1024 || pplm_args.discrim in ["sentiment", "clickbait"]
            tokenizer, model = get_gpt2_medium()
        else
            tokenizer, model = get_gpt2()
        end
    end

    if pplm_args.cuda[1]==true
        global device
        device = gpu
        set_device(pplm_args.cuda[2])
    end

    if add_eos_start==true
        input_ = [tokenizer.eos_token_id; tokenizer(prompt)]
    else
        input_ = tokenizer(prompt)
    end
    args=pplm_args
    model = model |> device
    @showprogress "Generating..." for i in 1:args.max_length
        input_ids = reshape(input_[:], :, 1) |> device
        outputs = model(input_ids; output_attentions=false,
                            output_hidden_states=true,
                            use_cache=false);
        original_logits = outputs.logits[:, end, 1]
        original_probs = PPLM.temp_softmax(original_logits; t=args.temperature)
        if pplm_args.method=="BoW"
            if pplm_args.perturb=="probs"
                pert_probs = perturb_probs(original_probs, tokenizer, pplm_args)
            elseif pplm_args.perturb == "hidden"
                hidden = outputs.hidden_states[end]
                modified_hidden = perturb_hidden_bow(hidden, model, tokenizer, pplm_args)
                pert_logits = model.lm_head(modified_hidden)[:, end, 1]
                pert_probs = temp_softmax(pert_logits; t=pplm_args.temperature)
            else
                inp = input_ids[1:end-1, :]
                prev = input_ids[end:end, :]
                outputs = model(inp; output_attentions=false,
                                    output_hidden_states=true,
                                    use_cache=true);
                past = outputs.past_key_values;
                original_logits_ = outputs.logits[:, end, 1]
                original_probs_ = temp_softmax(original_logits_; t=pplm_args.temperature)
                past = outputs.past_key_values;
                new_past = perturb_past_bow(model, prev, past, original_probs_, tokenizer, pplm_args)
                output_new = model(prev; past_key_values=new_past,
                                                    output_attentions=false,
                                                    output_hidden_states=true,
                                                    use_cache=true);
                pert_logits = output_new.logits[:, end, 1]
                pert_probs = temp_softmax(pert_logits; t=args.temperature)
            end
        else
            if pplm_args.perturb == "hidden"
                hidden = outputs.hidden_states[end]
                modified_hidden = perturb_hidden_discrim(hidden, model, tokenizer, pplm_args)
                pert_logits = model.lm_head(modified_hidden)[:, end, 1]
                pert_probs = temp_softmax(pert_logits; t=pplm_args.temperature)
            else
                inp = input_ids[1:end-1, :]
                prev = input_ids[end:end, :]
                outputs = model(inp; output_attentions=false,
                                    output_hidden_states=true,
                                    use_cache=true);
                past = outputs.past_key_values;
                original_logits_ = outputs.logits[:, end, 1]
                original_probs_ = temp_softmax(original_logits_; t=pplm_args.temperature)
                past = outputs.past_key_values;
                new_past = perturb_past_discrim(model, prev, past, original_probs_, pplm_args)
                output_new = model(prev; past_key_values=new_past,
                                                    output_attentions=false,
                                                    output_hidden_states=true,
                                                    use_cache=true);
                pert_logits = output_new.logits[:, end, 1]
                pert_probs = temp_softmax(pert_logits; t=args.temperature)
            end
        end
        gm_scale=pplm_args.fusion_gm_scale
        pert_probs = Float32.((original_probs.^(1-gm_scale)).*(pert_probs.^gm_scale)) |> cpu
        if sample_method=="top_k"
            new_token = top_k_sample(pert_probs; k=args.top_k)[1]
        elseif sample_method == "nucleus"
            new_token = nucleus_sample(pert_probs; p=args.top_p)[1]
        end

        push!(input_, new_token)
    end

    text = detokenize(tokenizer, input_)
    if add_eos_start==true
        text = text[length(tokenizer.eos_token)+1:end]
    end
    print("PPLM (",pplm_args.method,"-",pplm_args.perturb,") generated text: ", text, "\n")
    return text
end


text:
## %PREDICTMD_GENERATED_BY%

import PredictMDExtra
PredictMDExtra.import_all()

import PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "breast_cancer_biopsy", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin model comparison code

Kernel = LIBSVM.Kernel

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

smoted_training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_features_df.csv",
    )
smoted_training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_labels_df.csv",
    )
smoted_training_features_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_features_df_filename;
        type_detect_rows = 100,
        )
    )
smoted_training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_labels_df_filename;
        type_detect_rows = 100,
        )
    )

logistic_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "logistic_classifier.jld2",
    )
random_forest_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "random_forest_classifier.jld2",
    )
c_svc_svm_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "c_svc_svm_classifier.jld2",
    )
nu_svc_svm_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "nu_svc_svm_classifier.jld2",
    )
knet_mlp_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "knet_mlp_classifier.jld2",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
logistic_classifier = nothing
Test.@test isnothing(logistic_classifier)

random_forest_classifier = nothing
Test.@test isnothing(random_forest_classifier)

c_svc_svm_classifier = nothing
Test.@test isnothing(c_svc_svm_classifier)

nu_svc_svm_classifier = nothing
Test.@test isnothing(nu_svc_svm_classifier)

knet_mlp_classifier = nothing
Test.@test isnothing(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

logistic_classifier =
    PredictMD.load_model(logistic_classifier_filename)
random_forest_classifier =
    PredictMD.load_model(random_forest_classifier_filename)
c_svc_svm_classifier =
    PredictMD.load_model(c_svc_svm_classifier_filename)
nu_svc_svm_classifier =
    PredictMD.load_model(nu_svc_svm_classifier_filename)
knet_mlp_classifier =
    PredictMD.load_model(knet_mlp_classifier_filename)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

all_models = PredictMD.AbstractFittable[
    logistic_classifier,
    random_forest_classifier,
    c_svc_svm_classifier,
    nu_svc_svm_classifier,
    knet_mlp_classifier,
    ]

single_label_name = :Class
negative_class = "benign"
positive_class = "malignant"

single_label_levels = [negative_class, positive_class]

categorical_label_names = Symbol[single_label_name]
continuous_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

println(
    logger_stream, string(
        "Single label binary classification metrics, training set, ",
        "fix sensitivity",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        training_features_df,
        training_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, training set, ",
        "fix specificity",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        training_features_df,
        training_labels_df,
        single_label_name,
        positive_class;
        specificity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, training set, ",
        "maximize F1 score",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        training_features_df,
        training_labels_df,
        single_label_name,
        positive_class;
        maximize = :f1score,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, training set, ",
        "maximize Cohen's kappa",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        training_features_df,
        training_labels_df,
        single_label_name,
        positive_class;
        maximize = :cohen_kappa,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, testing set, ",
        "fix sensitivity",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
metrics = PredictMD.singlelabelbinaryclassificationmetrics(all_models,
                                                 testing_features_df,
                                                 testing_labels_df,
                                                 single_label_name,
                                                 positive_class;
                                                 sensitivity = 0.95)
auprc_row = first(
    findall(
        strip.(metrics[:metric]) .== "AUPRC"
        )
    )
Test.@test(
    strip(metrics[auprc_row, :metric]) == "AUPRC"
    )
Test.@test(
    metrics[auprc_row, Symbol("Logistic regression")] > 0.910
    )
Test.@test(
    metrics[auprc_row, Symbol("Random forest")] > 0.910
    )
Test.@test(
    metrics[auprc_row, Symbol("SVM (C-SVC)")] > 0.910
    )
Test.@test(
    metrics[auprc_row, Symbol("SVM (nu-SVC)")] > 0.910
    )
Test.@test(
    metrics[auprc_row, Symbol("Knet MLP")] > 0.910
    )
aurocc_row = first(
    findall(
        strip.(metrics[:metric]) .== "AUROCC"
        )
    )
Test.@test(
    metrics[aurocc_row, :metric] == "AUROCC"
    )
Test.@test(
    metrics[aurocc_row, Symbol("Logistic regression")] > 0.910
    )
Test.@test(
    metrics[aurocc_row, Symbol("Random forest")] > 0.910
    )
Test.@test(
    metrics[aurocc_row, Symbol("SVM (C-SVC)")] > 0.910
    )
Test.@test(
    metrics[aurocc_row, Symbol("SVM (nu-SVC)")] > 0.910
    )
Test.@test(
    metrics[aurocc_row, Symbol("Knet MLP")] > 0.910
    )
avg_precision_row = first(
    findall(
        strip.(metrics[:metric]) .== "Average precision"
        )
    )
Test.@test(
    strip(metrics[avg_precision_row, :metric]) == "Average precision"
    )
Test.@test(
    metrics[avg_precision_row, Symbol("Logistic regression")] > 0.910
    )
Test.@test(
    metrics[avg_precision_row, Symbol("Random forest")] > 0.910
    )
Test.@test(
    metrics[avg_precision_row, Symbol("SVM (C-SVC)")] > 0.910
    )
Test.@test(
    metrics[avg_precision_row, Symbol("SVM (nu-SVC)")] > 0.910
    )
Test.@test(
    metrics[avg_precision_row, Symbol("Knet MLP")] > 0.910
    )
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

println(
    logger_stream, string(
        "Single label binary classification metrics, testing set, ",
        "fix specificity",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        specificity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, testing set, ",
        "maximize F1 score",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        maximize = :f1score,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(
    logger_stream, string(
        "Single label binary classification metrics, testing set, ",
        "maximize Cohen's kappa",
        )
    )
show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        all_models,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        maximize = :cohen_kappa,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

rocplottesting = PredictMD.plotroccurve(
    all_models,
    testing_features_df,
    testing_labels_df,
    single_label_name,
    positive_class,
    );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "rocplottesting",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, rocplottesting)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(rocplottesting)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "rocplottesting.pdf",
        ),
    rocplottesting,
    )

prplottesting = PredictMD.plotprcurve(
    all_models,
    testing_features_df,
    testing_labels_df,
    single_label_name,
    positive_class,
    );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "prplottesting",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, prplottesting)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(prplottesting)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "prplottesting.pdf",
        ),
    prplottesting,
    )

### End model comparison code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "breast_cancer_biopsy", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
module  TestIO

using PopGenCore
using Test

#cats_gen = normpath(joinpath(@__DIR__,"data/", "nancycats.gen"))
#sharks_csv = normpath(joinpath(@__DIR__,"data/", "gulfsharks.csv"))
example_vcf = normpath(joinpath(@__DIR__,"../src/precompile", "precompile.vcf.gz"))
#example_plink = normpath(joinpath(@__DIR__,"../src/precompile", "precompile.vcf.gz"))

@testset "File IO" begin
    @testset "Genepop io" begin
        @test typeof(@nancycats) == PopData
    end

    @testset "delimited io" begin
        @test typeof(@gulfsharks) == PopData
    end

    @testset "VCF io" begin
        @test typeof(vcf(example_vcf, silent = true)) == PopData
    end
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ cfea9340-3fd6-11eb-26cc-57e76c7353cd
using DSJulia

# ╔═╡ d00b3712-3fd6-11eb-354e-3182c3cb8eb1
md"""## Some examples on how to use DSJulia"""

# ╔═╡ 4e4c000e-3fda-11eb-3b32-8fad55386401
md"At the start of the notebook we need the name and email of the students"

# ╔═╡ eb3259ee-3fdb-11eb-39d1-6d8b2f5a29ed
student = (name="Bram De Jaegher", email="Bram.DeJaegher@gmail.com")

# ╔═╡ 28e74fe2-3fdc-11eb-28b5-1fc8d35e65cf
md"Next step is to initialise the tracker"

# ╔═╡ 1a6a2844-3fdd-11eb-2931-dd7fca88e484
tracker = ProgressTracker(student.name, student.email);

# ╔═╡ 4707f4ae-3fdc-11eb-3ebd-2b483041d459
md"-----"

# ╔═╡ e27e6aa0-2dab-11eb-3ccc-43c68f37114b
myclamp(x) = max.(0.0, min.(x, 1.0))

# ╔═╡ 9bbee332-4170-11eb-05a6-4998f14d307e
begin
	q₁ = Question(;
		description=md"""
		Complete the function `myclamp(x)` that clamps a number `x` between 0 and 1.

		Open assignments always return `missing`.
		""",
		validators= @safe[myclamp(-1)==0, myclamp(0.3)==0.3, myclamp(1.1)==1.0]
	)

	q₂ = QuestionOptional{Easy}(;
		description=md"""
		Try to make the clamping also work for arrays.
		""",
		validators= @safe[myclamp([2.0, 0.3, -10])==[1.0, 0.3, 0.0]]
	)

	q₃ = QuestionOptional{Intermediate}(;
		description=md"""
		This is an intermediate question. Surely you can complete this
		""",
		validators= @safe[true]
	)

	q₄ = QuestionOptional{Hard}(;
		description=md"""
		I admit, this one is definitely harder
		""",
		validators= @safe[false]
	)

   qb = QuestionBlock(;
	title=md"### Question 1.0: What a crazy exercise",
	description=md"""
		Some additional general kind off description and all.
		You can even add other admonitions.
		$(fyi(md"I'm here!"))

		Anything markdowny. Just make sure to use the triple accolades `\"\"\"`.

		""",
	questions = [q₁, q₂, q₃, q₄],
	hints=[	hint(md"Have you tried this?"),
			hint(md"Have you tried switching it on and off again?")]
	);

	validate(qb)
end

# ╔═╡ 87e6c2a8-2dac-11eb-33d3-77a35fc13d71
myclamp(1.1)

# ╔═╡ e05ee714-3fde-11eb-0209-e30f44655997
md"The automatic grading happens in the second argument of `validate()` and by passing the tracker-instance in the first argument the state of the question (correct/incorrect) is tracked"

# ╔═╡ 221b27f8-3fdf-11eb-06c1-f97060a68bca
md"**Hints:**"

# ╔═╡ 542d5fa6-2da9-11eb-1037-3b35a5b22bd5
hint(md"Did you think of this?")

# ╔═╡ 2e7109a0-3fdf-11eb-28af-0d99a01c5066
md"**FYI's:**"

# ╔═╡ 3cd3e814-3fdf-11eb-1918-75993ba2e121
fyi(md"""Did you know that scientists recently discovered pink UV-fluorescent squirrels?""")

# ╔═╡ c64b613a-49bc-11eb-3689-93a41d001eb1
md"**Recursive FYI's**"

# ╔═╡ d25f6fac-49bc-11eb-2379-5fff9c50e09b
fyi(md"""Did you know that scientists recently discovered pink UV-fluorescent squirrels? $(fyi(md"OK, this is a pretty funny accidental feature."))""")

# ╔═╡ Cell order:
# ╟─d00b3712-3fd6-11eb-354e-3182c3cb8eb1
# ╠═cfea9340-3fd6-11eb-26cc-57e76c7353cd
# ╟─4e4c000e-3fda-11eb-3b32-8fad55386401
# ╠═eb3259ee-3fdb-11eb-39d1-6d8b2f5a29ed
# ╠═28e74fe2-3fdc-11eb-28b5-1fc8d35e65cf
# ╠═1a6a2844-3fdd-11eb-2931-dd7fca88e484
# ╠═4707f4ae-3fdc-11eb-3ebd-2b483041d459
# ╠═9bbee332-4170-11eb-05a6-4998f14d307e
# ╠═e27e6aa0-2dab-11eb-3ccc-43c68f37114b
# ╠═87e6c2a8-2dac-11eb-33d3-77a35fc13d71
# ╠═e05ee714-3fde-11eb-0209-e30f44655997
# ╟─221b27f8-3fdf-11eb-06c1-f97060a68bca
# ╠═542d5fa6-2da9-11eb-1037-3b35a5b22bd5
# ╟─2e7109a0-3fdf-11eb-28af-0d99a01c5066
# ╠═3cd3e814-3fdf-11eb-1918-75993ba2e121
# ╟─c64b613a-49bc-11eb-3689-93a41d001eb1
# ╠═d25f6fac-49bc-11eb-2379-5fff9c50e09b


text:
module JLLEnvs

using Pkg
using Pkg.Artifacts
using Downloads

include("utils.jl")

const JLL_ENV_SHARDS_URL = "https://raw.githubusercontent.com/JuliaPackaging/BinaryBuilderBase.jl/master/Artifacts.toml"
const JLL_ENV_SHARDS = Dict{String,Any}()

function __init__()
    if haskey(ENV, "JULIA_CLANG_SHARDS_URL") && !isempty(get(ENV, "JULIA_CLANG_SHARDS_URL", ""))
        merge!(JLL_ENV_SHARDS, Artifacts.load_artifacts_toml(ENV["JULIA_CLANG_SHARDS_URL"]))
    else
        merge!(JLL_ENV_SHARDS, Artifacts.load_artifacts_toml(Downloads.download(JLL_ENV_SHARDS_URL)))
    end
end

const JLL_ENV_HOST_TRIPLE = "x86_64-linux-musl"
const JLL_ENV_GCC_SHARD_NAME = "GCCBootstrap"
const JLL_ENV_SYSTEM_SHARD_NAME = "PlatformSupport"

const JLL_ENV_GCC_VERSIONS = VersionNumber[
    v"4.8.5",
    v"5.2.0",
    v"6.1.0",
    v"7.1.0",
    v"8.1.0",
    v"9.1.0",
    v"10.2.0",
    v"11.0.0-iains",
]

# External API.
# The triples used for cross-platform configuration
const JLL_ENV_TRIPLES = String[
    "aarch64-apple-darwin20",
    "aarch64-linux-gnu",           # Tier 1
    "aarch64-linux-musl",
    "armv7l-linux-gnueabihf",
    "armv7l-linux-musleabihf",
    "i686-linux-gnu",              # Tier 1
    "i686-linux-musl",
    "i686-w64-mingw32",            # Tier 2
    "powerpc64le-linux-gnu",
    "x86_64-apple-darwin14",       # Tier 1
    "x86_64-linux-gnu",            # Tier 1
    "x86_64-linux-musl",
    "x86_64-unknown-freebsd",      # Tier 1
    "x86_64-w64-mingw32",          # Tier 1
]

const JLL_ENV_CLANG_TARGETS_MAPPING = Dict(
    "aarch64-apple-darwin20"=>"aarch64-apple-darwin20",
    "aarch64-linux-gnu"=>"aarch64-unknown-linux-gnu",
    "aarch64-linux-musl"=>"aarch64-unknown-linux-musl",
    "armv7l-linux-gnueabihf"=>"armv7l-unknown-linux-gnueabihf",
    "armv7l-linux-musleabihf"=>"armv7l-unknown-linux-musleabihf",
    "i686-linux-gnu"=>"i686-unknown-linux-gnu",
    "i686-linux-musl"=>"i686-unknown-linux-musl",
    "i686-w64-mingw32"=>"i686-w64-windows-gnu",
    "powerpc64le-linux-gnu"=>"powerpc64le-unknown-linux-gnu",
    "x86_64-apple-darwin14"=>"x86_64-apple-darwin14",
    "x86_64-linux-gnu"=>"x86_64-unknown-linux-gnu",
    "x86_64-linux-musl"=>"x86_64-unknown-linux-musl",
    "x86_64-unknown-freebsd"=>"x86_64-unknown-freebsd12.2",
    "x86_64-w64-mingw32"=>"x86_64-w64-windows-gnu",
)

triple2target(triple::String) = get(JLL_ENV_CLANG_TARGETS_MAPPING, triple, "unknown")

const JLL_ENV_CLANG_TRIPLE_ABBREVATION_MAPPING = Dict(
    "x86_64-unknown-freebsd"=>"x86_64-unknown-freebsd12.2"
)

# expand freebsd to freebsd12.2
expand_triple(triple::AbstractString) = get(JLL_ENV_CLANG_TRIPLE_ABBREVATION_MAPPING, triple, triple)

# Internal use, triples used in artifacts
const JLL_ENV_ARTIFACT_TRIPLES = expand_triple.(JLL_ENV_TRIPLES)

function get_gcc_shard_key(triple::String, version::VersionNumber=v"4.8.5")
    @assert version ∈ JLL_ENV_GCC_VERSIONS "Wrong JLL gcc version: $version. Please choose a version listed in `JLL_ENV_GCC_VERSIONS`."
    @assert triple ∈ JLL_ENV_ARTIFACT_TRIPLES "Wrong JLL target triple: $triple. Please choose a triple listed in `JLL_ENV_ARTIFACT_TRIPLES`."
    # ignore gcc version for Apple Silicon
    if triple == "aarch64-apple-darwin20"
        return "$JLL_ENV_GCC_SHARD_NAME-aarch64-apple-darwin20.v11.0.0-iains.$JLL_ENV_HOST_TRIPLE.unpacked"
    else
        return "$JLL_ENV_GCC_SHARD_NAME-$triple.v$version.$JLL_ENV_HOST_TRIPLE.unpacked"
    end
end

function get_system_shard_key(triple::String)
    @assert triple ∈ JLL_ENV_ARTIFACT_TRIPLES "Wrong JLL target triple: $triple. Please choose a triple listed in `JLL_ENV_ARTIFACT_TRIPLES`."
    platform_keys = filter(collect(keys(JLL_ENV_SHARDS))) do key
        startswith(key, "$JLL_ENV_SYSTEM_SHARD_NAME-$triple.") &&
        endswith(key, "$JLL_ENV_HOST_TRIPLE.unpacked")
    end
    return platform_keys[]
end

function get_environment_info(triple::String, version::VersionNumber=v"4.8.5")
    @assert triple ∈ JLL_ENV_ARTIFACT_TRIPLES "Wrong JLL target triple: $triple. Please choose a triple listed in `JLL_ENV_ARTIFACT_TRIPLES`."
    gcc = JLL_ENV_SHARDS[get_gcc_shard_key(triple, version)][]
    sys = JLL_ENV_SHARDS[get_system_shard_key(triple)][]
    gcc_download = gcc["download"][]
    sys_download = sys["download"][]
    info = [
        (id=gcc["git-tree-sha1"], url=gcc_download["url"], chk=gcc_download["sha256"]),
        (id=sys["git-tree-sha1"], url=sys_download["url"], chk=sys_download["sha256"]),
    ]
    return info
end

function get_system_dirs(triple::String, version::VersionNumber=v"4.8.5")
    triple = expand_triple(triple)
    info = get_environment_info(triple, version)
    gcc_info = info[1]
    sys_info = info[2]

    # download shards
    if haskey(ENV, "JULIA_CLANG_SHARDS_URL") && !isempty(get(ENV, "JULIA_CLANG_SHARDS_URL", ""))
        @info "Downloading artifact($(gcc_info.id)) from $(gcc_info.url) ..."
    end
    Artifacts.download_artifact(Base.SHA1(gcc_info.id), gcc_info.url, gcc_info.chk)
    # Artifacts.download_artifact(Base.SHA1(sys_info.id), sys_info.url, sys_info.chk)

    # -isystem paths
    gcc_triple_path = Artifacts.artifact_path(Base.SHA1(gcc_info.id))
    # sys_triple_path = Artifacts.artifact_path(Base.SHA1(sys_info.id))
    isys = String[]
    if triple == "x86_64-apple-darwin14"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, triple, "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "usr", "include"))
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "System", "Library", "Frameworks"))
    elseif triple == "aarch64-apple-darwin20"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, "11.0.0", "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, "11.0.0", "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, triple, "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "usr", "include"))
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "System", "Library", "Frameworks"))
    elseif triple == "x86_64-w64-mingw32" || triple == "i686-w64-mingw32"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, triple, "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "include"))
    elseif triple == "i686-linux-gnu" || triple == "x86_64-linux-gnu" ||
            triple == "aarch64-linux-gnu" || triple == "powerpc64le-linux-gnu" ||
            triple == "x86_64-unknown-freebsd12.2"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, triple, "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "usr", "include"))
    elseif triple == "i686-linux-musl" || triple == "x86_64-linux-musl" ||
            triple == "aarch64-linux-musl"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, triple, "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, triple, "sys-root", "usr", "include"))
    elseif triple == "armv7l-linux-gnueabihf"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", "arm-linux-gnueabihf", string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", "arm-linux-gnueabihf", string(version), "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, "arm-linux-gnueabihf", "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, "arm-linux-gnueabihf", "sys-root", "usr", "include"))
    elseif triple == "armv7l-linux-musleabihf"
        # compiler
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", "arm-linux-musleabihf", string(version), "include"))
        push!(isys, joinpath(gcc_triple_path, "lib", "gcc", "arm-linux-musleabihf", string(version), "include-fixed"))
        push!(isys, joinpath(gcc_triple_path, "arm-linux-musleabihf", "include"))
        # sys-root
        push!(isys, joinpath(gcc_triple_path, "arm-linux-musleabihf", "sys-root", "usr", "include"))
    else
        error("Platform $triple is not supported.")
    end

    @assert all(isdir, isys) "failed to setup environment due to missing dirs, please file an issue at https://github.com/JuliaInterop/Clang.jl/issues."

    return normpath.(isys)
end

function get_pkg_artifact_dir(pkg::Module, target::String)
    afts = first(values(Artifacts.load_artifacts_toml(Artifacts.find_artifacts_toml(Pkg.pathof(pkg)))))
    target_arch, target_os, target_libc = get_arch_os_libc(target)
    candidates = Dict[]
    for info in afts
        if info isa Dict
            arch = get(info, "arch", "")
            os = get(info, "os", "")
            libc = get(info, "libc", "")
            if arch == target_arch && os == target_os && libc == target_libc
                push!(candidates, info)
            end
        else
            # this could be an "Any"-platform JLL package
            push!(candidates, afts)
            break
        end
    end
    isempty(candidates) && return ""
    length(candidates) > 1 && @warn "found more than one candidate artifacts, only use the first one: $(first(candidates))"
    info = first(candidates)
    download_info = info["download"][]
    id, url, chk = info["git-tree-sha1"], download_info["url"], download_info["sha256"]
    Artifacts.download_artifact(Base.SHA1(id), url, chk)
    return normpath(Artifacts.artifact_path(Base.SHA1(id)))
end

function get_pkg_include_dir(pkg::Module, target::String)
    artifact_dir = get_pkg_artifact_dir(pkg, target)
    if isempty(artifact_dir)
        return ""
    else
        joinpath(artifact_dir, "include")
    end
end

end # module


text:
include("blocksys.jl")
include("matrixgen.jl")
using .matrixgen
using LinearAlgebra

function benchmark(n, l, repeats, func)
    totalTime = 0
    totalMemory = 0
    for i in 1:repeats
        matrixgen.blockmat(n, l, 10.0, "test.txt")
        A = Blocksys.readMatrix("test.txt")
        b = Blocksys.computeRightSideVector(A[1], n, l)
        (_, time, memory) = @timed func(A[1], n, l, b)
        totalTime += time
        totalMemory += memory
    end
    println(n, ";",totalTime/repeats, ";", totalMemory/repeats)
end

function lu(A, n, l, b)
    L, U = Blocksys.lu(A, n, l)
    return Blocksys.solveLu(U, L, n, l, b)
end

function luChoice(A, n, l, b)
    U, L, perm = Blocksys.luWithChoice(A, n, l)
    return Blocksys.solveLuWithChoice(U, L, n, l, b, perm)
end

function classic(A, n, l, b)
    return \(Array(A), b)
end

functions = [classic, Blocksys.solveGauss, Blocksys.solveGaussWithChoice, lu, luChoice]

for j in 1:length(functions)
    println("Algorytm $j")
    # for k in 100:100:1000
    #     benchmark(k, 4, 1, functions[j])
    # end

    benchmark(1000, 4, 3, functions[j])

    for k in 5000:5000:50000
        benchmark(k, 4, 3, functions[j])
    end
end


text:
using Random

function RLCore.Experiment(::Val{:JuliaRL}, ::Val{:Minimax}, ::Val{:OpenSpiel}, game;)
    env = OpenSpielEnv(string(game))
    agents = (
        Agent(policy = MinimaxPolicy(), role = 0),
        Agent(policy = MinimaxPolicy(), role = 1),
    )
    hooks = (TotalRewardPerEpisode(), TotalRewardPerEpisode())
    description = """
      # Play `$game` in OpenSpiel with Minimax
      """
    Experiment(agents, env, StopAfterEpisode(1), hooks, description)
end

function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:TabularCFR},
    ::Val{:OpenSpiel},
    game;
    n_iter = 300,
    seed = 123,
)
    env = OpenSpielEnv(
        game;
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    rng = MersenneTwister(seed)
    π = TabularCFRPolicy(; n_iter = n_iter, env = env, rng = rng)

    agents = map(get_players(env)) do p
        if p == get_chance_player(env)
            Agent(; policy = RandomPolicy(), role = p)
        else
            Agent(; policy = π, role = p)
        end
    end

    hooks = [
        p == get_chance_player(env) ? EmptyHook() : TotalRewardPerEpisode()
        for p in get_players(env)
    ]
    description = """
      # Play `$game` in OpenSpiel with TabularCFRPolicy
      """
    Experiment(agents, env, StopAfterEpisode(100_000), hooks, description)
end


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct MultiStartCost <: OperationalCost
        variable::VariableCost
        no_load::Float64
        fixed::Float64
        start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
        shut_down::Float64
    end

Data Structure Operational Cost Data which includes fixed, variable cost, multiple start up cost and stop costs.

# Arguments
- `variable::VariableCost`: variable cost
- `no_load::Float64`: no load cost
- `fixed::Float64`: fixed cost
- `start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}`: start-up cost
- `shut_down::Float64`: shut-down cost, validation range: `(0, nothing)`, action if invalid: `warn`
"""
mutable struct MultiStartCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "no load cost"
    no_load::Float64
    "fixed cost"
    fixed::Float64
    "start-up cost"
    start_up::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
    "shut-down cost"
    shut_down::Float64
end

function MultiStartCost(; variable, no_load, fixed, start_up, shut_down, )
    MultiStartCost(variable, no_load, fixed, start_up, shut_down, )
end

# Constructor for demo purposes; non-functional.
function MultiStartCost(::Nothing)
    MultiStartCost(;
        variable=VariableCost((0.0, 0.0)),
        no_load=0.0,
        fixed=0.0,
        start_up=(hot = START_COST, warm = START_COST, cold = START_COST),
        shut_down=0.0,
    )
end

"""Get [`MultiStartCost`](@ref) `variable`."""
get_variable(value::MultiStartCost) = value.variable
"""Get [`MultiStartCost`](@ref) `no_load`."""
get_no_load(value::MultiStartCost) = value.no_load
"""Get [`MultiStartCost`](@ref) `fixed`."""
get_fixed(value::MultiStartCost) = value.fixed
"""Get [`MultiStartCost`](@ref) `start_up`."""
get_start_up(value::MultiStartCost) = value.start_up
"""Get [`MultiStartCost`](@ref) `shut_down`."""
get_shut_down(value::MultiStartCost) = value.shut_down

"""Set [`MultiStartCost`](@ref) `variable`."""
set_variable!(value::MultiStartCost, val) = value.variable = val
"""Set [`MultiStartCost`](@ref) `no_load`."""
set_no_load!(value::MultiStartCost, val) = value.no_load = val
"""Set [`MultiStartCost`](@ref) `fixed`."""
set_fixed!(value::MultiStartCost, val) = value.fixed = val
"""Set [`MultiStartCost`](@ref) `start_up`."""
set_start_up!(value::MultiStartCost, val) = value.start_up = val
"""Set [`MultiStartCost`](@ref) `shut_down`."""
set_shut_down!(value::MultiStartCost, val) = value.shut_down = val