text:
export BrickREDCreateSession
struct BrickREDCreateSession
    error_code::Integer
    session_id::Integer
end

export BrickREDAllocateString
struct BrickREDAllocateString
    error_code::Integer
    string_id::Integer
end

export BrickREDStringLength
struct BrickREDStringLength
    error_code::Integer
    length::Integer
end

export BrickREDStringChunk
struct BrickREDStringChunk
    error_code::Integer
    buffer::String
end

export BrickREDAllocateList
struct BrickREDAllocateList
    error_code::Integer
    list_id::Integer
end

export BrickREDListLength
struct BrickREDListLength
    error_code::Integer
    length::Integer
end

export BrickREDListItem
struct BrickREDListItem
    error_code::Integer
    item_object_id::Integer
    type::Integer
end

export BrickREDOpenFile
struct BrickREDOpenFile
    error_code::Integer
    file_id::Integer
end

export BrickREDCreatePipe
struct BrickREDCreatePipe
    error_code::Integer
    file_id::Integer
end

export BrickREDFileInfo
struct BrickREDFileInfo
    error_code::Integer
    type::Integer
    name_string_id::Integer
    flags::Integer
    permissions::Integer
    uid::Integer
    gid::Integer
    length::Integer
    access_timestamp::Integer
    modification_timestamp::Integer
    status_change_timestamp::Integer
end

export BrickREDReadFile
struct BrickREDReadFile
    error_code::Integer
    buffer::Vector{Integer}
    length_read::Integer
end

export BrickREDWriteFile
struct BrickREDWriteFile
    error_code::Integer
    length_written::Integer
end

export BrickREDSetFilePosition
struct BrickREDSetFilePosition
    error_code::Integer
    position::Integer
end

export BrickREDFilePosition
struct BrickREDFilePosition
    error_code::Integer
    position::Integer
end

export BrickREDFileEvents
struct BrickREDFileEvents
    error_code::Integer
    events::Integer
end

export BrickREDOpenDirectory
struct BrickREDOpenDirectory
    error_code::Integer
    directory_id::Integer
end

export BrickREDDirectoryName
struct BrickREDDirectoryName
    error_code::Integer
    name_string_id::Integer
end

export BrickREDNextDirectoryEntry
struct BrickREDNextDirectoryEntry
    error_code::Integer
    name_string_id::Integer
    type::Integer
end

export BrickREDProcesses
struct BrickREDProcesses
    error_code::Integer
    processes_list_id::Integer
end

export BrickREDSpawnProcess
struct BrickREDSpawnProcess
    error_code::Integer
    process_id::Integer
end

export BrickREDProcessCommand
struct BrickREDProcessCommand
    error_code::Integer
    executable_string_id::Integer
    arguments_list_id::Integer
    environment_list_id::Integer
    working_directory_string_id::Integer
end

export BrickREDProcessIdentity
struct BrickREDProcessIdentity
    error_code::Integer
    pid::Integer
    uid::Integer
    gid::Integer
end

export BrickREDProcessStdio
struct BrickREDProcessStdio
    error_code::Integer
    stdin_file_id::Integer
    stdout_file_id::Integer
    stderr_file_id::Integer
end

export BrickREDProcessState
struct BrickREDProcessState
    error_code::Integer
    state::Integer
    timestamp::Integer
    exit_code::Integer
end

export BrickREDPrograms
struct BrickREDPrograms
    error_code::Integer
    programs_list_id::Integer
end

export BrickREDDefineProgram
struct BrickREDDefineProgram
    error_code::Integer
    program_id::Integer
end

export BrickREDProgramIdentifier
struct BrickREDProgramIdentifier
    error_code::Integer
    identifier_string_id::Integer
end

export BrickREDProgramRootDirectory
struct BrickREDProgramRootDirectory
    error_code::Integer
    root_directory_string_id::Integer
end

export BrickREDProgramCommand
struct BrickREDProgramCommand
    error_code::Integer
    executable_string_id::Integer
    arguments_list_id::Integer
    environment_list_id::Integer
    working_directory_string_id::Integer
end

export BrickREDProgramStdioRedirection
struct BrickREDProgramStdioRedirection
    error_code::Integer
    stdin_redirection::Integer
    stdin_file_name_string_id::Integer
    stdout_redirection::Integer
    stdout_file_name_string_id::Integer
    stderr_redirection::Integer
    stderr_file_name_string_id::Integer
end

export BrickREDProgramSchedule
struct BrickREDProgramSchedule
    error_code::Integer
    start_mode::Integer
    continue_after_error::Bool
    start_interval::Integer
    start_fields_string_id::Integer
end

export BrickREDProgramSchedulerState
struct BrickREDProgramSchedulerState
    error_code::Integer
    state::Integer
    timestamp::Integer
    message_string_id::Integer
end

export BrickREDLastSpawnedProgramProcess
struct BrickREDLastSpawnedProgramProcess
    error_code::Integer
    process_id::Integer
    timestamp::Integer
end

export BrickREDCustomProgramOptionNames
struct BrickREDCustomProgramOptionNames
    error_code::Integer
    names_list_id::Integer
end

export BrickREDCustomProgramOptionValue
struct BrickREDCustomProgramOptionValue
    error_code::Integer
    value_string_id::Integer
end

export BrickREDIdentity
struct BrickREDIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickRED
"""
Executes user programs and controls other Bricks/Bricklets standalone
"""
mutable struct BrickRED <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickRED(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.brick_red")
        deviceInternal = package.BrickRED(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export create_session
"""
    $(SIGNATURES)

"""
function create_session(device::BrickRED, lifetime)
    return device.deviceInternal.create_session(lifetime)
end

export expire_session
"""
    $(SIGNATURES)

"""
function expire_session(device::BrickRED, session_id)
    return device.deviceInternal.expire_session(session_id)
end

export expire_session_unchecked
"""
    $(SIGNATURES)

"""
function expire_session_unchecked(device::BrickRED, session_id)
    device.deviceInternal.expire_session_unchecked(session_id)
end

export keep_session_alive
"""
    $(SIGNATURES)

"""
function keep_session_alive(device::BrickRED, session_id, lifetime)
    return device.deviceInternal.keep_session_alive(session_id, lifetime)
end

export release_object
"""
    $(SIGNATURES)

Decreases the reference count of an object by one and returns the resulting
error code. If the reference count reaches zero the object gets destroyed.
"""
function release_object(device::BrickRED, object_id, session_id)
    return device.deviceInternal.release_object(object_id, session_id)
end

export release_object_unchecked
"""
    $(SIGNATURES)

"""
function release_object_unchecked(device::BrickRED, object_id, session_id)
    device.deviceInternal.release_object_unchecked(object_id, session_id)
end

export allocate_string
"""
    $(SIGNATURES)

Allocates a new string object, reserves ``length_to_reserve`` bytes memory
for it and sets up to the first 60 bytes. Set ``length_to_reserve`` to the
length of the string that should be stored in the string object.

Returns the object ID of the new string object and the resulting error code.
"""
function allocate_string(device::BrickRED, length_to_reserve, buffer, session_id)
    return device.deviceInternal.allocate_string(length_to_reserve, buffer, session_id)
end

export truncate_string
"""
    $(SIGNATURES)

Truncates a string object to ``length`` bytes and returns the resulting
error code.
"""
function truncate_string(device::BrickRED, string_id, length)
    return device.deviceInternal.truncate_string(string_id, length)
end

export get_string_length
"""
    $(SIGNATURES)

Returns the length of a string object and the resulting error code.
"""
function get_string_length(device::BrickRED, string_id)
    return device.deviceInternal.get_string_length(string_id)
end

export set_string_chunk
"""
    $(SIGNATURES)

Sets a chunk of up to 58 bytes in a string object beginning at ``offset``.

Returns the resulting error code.
"""
function set_string_chunk(device::BrickRED, string_id, offset, buffer)
    return device.deviceInternal.set_string_chunk(string_id, offset, buffer)
end

export get_string_chunk
"""
    $(SIGNATURES)

Returns a chunk up to 63 bytes from a string object beginning at ``offset`` and
returns the resulting error code.
"""
function get_string_chunk(device::BrickRED, string_id, offset)
    return device.deviceInternal.get_string_chunk(string_id, offset)
end

export allocate_list
"""
    $(SIGNATURES)

Allocates a new list object and reserves memory for ``length_to_reserve``
items. Set ``length_to_reserve`` to the number of items that should be stored
in the list object.

Returns the object ID of the new list object and the resulting error code.

When a list object gets destroyed then the reference count of each object in
the list object is decreased by one.
"""
function allocate_list(device::BrickRED, length_to_reserve, session_id)
    return device.deviceInternal.allocate_list(length_to_reserve, session_id)
end

export get_list_length
"""
    $(SIGNATURES)

Returns the length of a list object in items and the resulting error code.
"""
function get_list_length(device::BrickRED, list_id)
    return device.deviceInternal.get_list_length(list_id)
end

export get_list_item
"""
    $(SIGNATURES)

Returns the object ID and type of the object stored at ``index`` in a list
object and returns the resulting error code.

Possible object types are:

* String = 0
* List = 1
* File = 2
* Directory = 3
* Process = 4
* Program = 5
"""
function get_list_item(device::BrickRED, list_id, index, session_id)
    return device.deviceInternal.get_list_item(list_id, index, session_id)
end

export append_to_list
"""
    $(SIGNATURES)

Appends an object to a list object and increases the reference count of the
appended object by one.

Returns the resulting error code.
"""
function append_to_list(device::BrickRED, list_id, item_object_id)
    return device.deviceInternal.append_to_list(list_id, item_object_id)
end

export remove_from_list
"""
    $(SIGNATURES)

Removes the object stored at ``index`` from a list object and decreases the
reference count of the removed object by one.

Returns the resulting error code.
"""
function remove_from_list(device::BrickRED, list_id, index)
    return device.deviceInternal.remove_from_list(list_id, index)
end

export open_file
"""
    $(SIGNATURES)

Opens an existing file or creates a new file and allocates a new file object
for it.

FIXME: name has to be absolute

The reference count of the name string object is increased by one. When the
file object gets destroyed then the reference count of the name string object is
decreased by one. Also the name string object is locked and cannot be modified
while the file object holds a reference to it.

The ``flags`` parameter takes a ORed combination of the following possible file
flags (in hexadecimal notation):

* ReadOnly = 0x0001 (O_RDONLY)
* WriteOnly = 0x0002 (O_WRONLY)
* ReadWrite = 0x0004 (O_RDWR)
* Append = 0x0008 (O_APPEND)
* Create = 0x0010 (O_CREAT)
* Exclusive = 0x0020 (O_EXCL)
* NonBlocking = 0x0040 (O_NONBLOCK)
* Truncate = 0x0080 (O_TRUNC)
* Temporary = 0x0100
* Replace = 0x0200

FIXME: explain *Temporary* and *Replace* flag

The ``permissions`` parameter takes a ORed combination of the following
possible file permissions (in octal notation) that match the common UNIX
permission bits:

* UserRead = 00400
* UserWrite = 00200
* UserExecute = 00100
* GroupRead = 00040
* GroupWrite = 00020
* GroupExecute = 00010
* OthersRead = 00004
* OthersWrite = 00002
* OthersExecute = 00001

Returns the object ID of the new file object and the resulting error code.
"""
function open_file(device::BrickRED, name_string_id, flags, permissions, uid, gid, session_id)
    return device.deviceInternal.open_file(name_string_id, flags, permissions, uid, gid, session_id)
end

export create_pipe
"""
    $(SIGNATURES)

Creates a new pipe and allocates a new file object for it.

The ``flags`` parameter takes a ORed combination of the following possible
pipe flags (in hexadecimal notation):

* NonBlockingRead = 0x0001
* NonBlockingWrite = 0x0002

The length of the pipe buffer can be specified with the ``length`` parameter
in bytes. If length is set to zero, then the default pipe buffer length is used.

Returns the object ID of the new file object and the resulting error code.
"""
function create_pipe(device::BrickRED, flags, length, session_id)
    return device.deviceInternal.create_pipe(flags, length, session_id)
end

export get_file_info
"""
    $(SIGNATURES)

Returns various information about a file and the resulting error code.

Possible file types are:

* Unknown = 0
* Regular = 1
* Directory = 2
* Character = 3
* Block = 4
* FIFO = 5
* Symlink = 6
* Socket = 7
* Pipe = 8

If the file type is *Pipe* then the returned name string object is invalid,
because a pipe has no name. Otherwise the returned name string object was used
to open or create the file object, as passed to :func:`Open File`.

The returned flags were used to open or create the file object, as passed to
:func:`Open File` or :func:`Create Pipe`. See the respective function for a list
of possible file and pipe flags.

FIXME: everything except flags and length is invalid if file type is *Pipe*
"""
function get_file_info(device::BrickRED, file_id, session_id)
    return device.deviceInternal.get_file_info(file_id, session_id)
end

export read_file
"""
    $(SIGNATURES)

Reads up to 62 bytes from a file object.

Returns the bytes read, the actual number of bytes read and the resulting
error code.

If there is not data to be read, either because the file position reached
end-of-file or because there is not data in the pipe, then zero bytes are
returned.

If the file object was created by :func:`Open File` without the *NonBlocking*
flag or by :func:`Create Pipe` without the *NonBlockingRead* flag then the
error code *NotSupported* is returned.
"""
function read_file(device::BrickRED, file_id, length_to_read)
    return device.deviceInternal.read_file(file_id, length_to_read)
end

export read_file_async
"""
    $(SIGNATURES)

Reads up to 2\\ :sup:`63`\\  - 1 bytes from a file object asynchronously.

Reports the bytes read (in 60 byte chunks), the actual number of bytes read and
the resulting error code via the :cb:`Async File Read` callback.

If there is not data to be read, either because the file position reached
end-of-file or because there is not data in the pipe, then zero bytes are
reported.

If the file object was created by :func:`Open File` without the *NonBlocking*
flag or by :func:`Create Pipe` without the *NonBlockingRead* flag then the error
code *NotSupported* is reported via the :cb:`Async File Read` callback.
"""
function read_file_async(device::BrickRED, file_id, length_to_read)
    device.deviceInternal.read_file_async(file_id, length_to_read)
end

export abort_async_file_read
"""
    $(SIGNATURES)

Aborts a :func:`Read File Async` operation in progress.

Returns the resulting error code.

On success the :cb:`Async File Read` callback will report *OperationAborted*.
"""
function abort_async_file_read(device::BrickRED, file_id)
    return device.deviceInternal.abort_async_file_read(file_id)
end

export write_file
"""
    $(SIGNATURES)

Writes up to 61 bytes to a file object.

Returns the actual number of bytes written and the resulting error code.

If the file object was created by :func:`Open File` without the *NonBlocking*
flag or by :func:`Create Pipe` without the *NonBlockingWrite* flag then the
error code *NotSupported* is returned.
"""
function write_file(device::BrickRED, file_id, buffer, length_to_write)
    return device.deviceInternal.write_file(file_id, buffer, length_to_write)
end

export write_file_unchecked
"""
    $(SIGNATURES)

Writes up to 61 bytes to a file object.

Does neither report the actual number of bytes written nor the resulting error
code.

If the file object was created by :func:`Open File` without the *NonBlocking*
flag or by :func:`Create Pipe` without the *NonBlockingWrite* flag then the
write operation will fail silently.
"""
function write_file_unchecked(device::BrickRED, file_id, buffer, length_to_write)
    device.deviceInternal.write_file_unchecked(file_id, buffer, length_to_write)
end

export write_file_async
"""
    $(SIGNATURES)

Writes up to 61 bytes to a file object.

Reports the actual number of bytes written and the resulting error code via the
:cb:`Async File Write` callback.

If the file object was created by :func:`Open File` without the *NonBlocking*
flag or by :func:`Create Pipe` without the *NonBlockingWrite* flag then the
error code *NotSupported* is reported via the :cb:`Async File Write` callback.
"""
function write_file_async(device::BrickRED, file_id, buffer, length_to_write)
    device.deviceInternal.write_file_async(file_id, buffer, length_to_write)
end

export set_file_position
"""
    $(SIGNATURES)

Set the current seek position of a file object relative to ``origin``.

Possible file origins are:

* Beginning = 0
* Current = 1
* End = 2

Returns the resulting absolute seek position and error code.

If the file object was created by :func:`Create Pipe` then it has no seek
position and the error code *InvalidSeek* is returned.
"""
function set_file_position(device::BrickRED, file_id, offset, origin)
    return device.deviceInternal.set_file_position(file_id, offset, origin)
end

export get_file_position
"""
    $(SIGNATURES)

Returns the current seek position of a file object and returns the
resulting error code.

If the file object was created by :func:`Create Pipe` then it has no seek
position and the error code *InvalidSeek* is returned.
"""
function get_file_position(device::BrickRED, file_id)
    return device.deviceInternal.get_file_position(file_id)
end

export set_file_events
"""
    $(SIGNATURES)

"""
function set_file_events(device::BrickRED, file_id, events)
    return device.deviceInternal.set_file_events(file_id, events)
end

export get_file_events
"""
    $(SIGNATURES)

"""
function get_file_events(device::BrickRED, file_id)
    return device.deviceInternal.get_file_events(file_id)
end

export open_directory
"""
    $(SIGNATURES)

Opens an existing directory and allocates a new directory object for it.

FIXME: name has to be absolute

The reference count of the name string object is increased by one. When the
directory object is destroyed then the reference count of the name string
object is decreased by one. Also the name string object is locked and cannot be
modified while the directory object holds a reference to it.

Returns the object ID of the new directory object and the resulting error code.
"""
function open_directory(device::BrickRED, name_string_id, session_id)
    return device.deviceInternal.open_directory(name_string_id, session_id)
end

export get_directory_name
"""
    $(SIGNATURES)

Returns the name of a directory object, as passed to :func:`Open Directory`, and
the resulting error code.
"""
function get_directory_name(device::BrickRED, directory_id, session_id)
    return device.deviceInternal.get_directory_name(directory_id, session_id)
end

export get_next_directory_entry
"""
    $(SIGNATURES)

Returns the next entry in a directory object and the resulting error code.

If there is not next entry then error code *NoMoreData* is returned. To rewind
a directory object call :func:`Rewind Directory`.

Possible directory entry types are:

* Unknown = 0
* Regular = 1
* Directory = 2
* Character = 3
* Block = 4
* FIFO = 5
* Symlink = 6
* Socket = 7
"""
function get_next_directory_entry(device::BrickRED, directory_id, session_id)
    return device.deviceInternal.get_next_directory_entry(directory_id, session_id)
end

export rewind_directory
"""
    $(SIGNATURES)

Rewinds a directory object and returns the resulting error code.
"""
function rewind_directory(device::BrickRED, directory_id)
    return device.deviceInternal.rewind_directory(directory_id)
end

export create_directory
"""
    $(SIGNATURES)

FIXME: name has to be absolute
"""
function create_directory(device::BrickRED, name_string_id, flags, permissions, uid, gid)
    return device.deviceInternal.create_directory(name_string_id, flags, permissions, uid, gid)
end

export get_processes
"""
    $(SIGNATURES)

"""
function get_processes(device::BrickRED, session_id)
    return device.deviceInternal.get_processes(session_id)
end

export spawn_process
"""
    $(SIGNATURES)

"""
function spawn_process(device::BrickRED, executable_string_id, arguments_list_id, environment_list_id, working_directory_string_id, uid, gid, stdin_file_id, stdout_file_id, stderr_file_id, session_id)
    return device.deviceInternal.spawn_process(executable_string_id, arguments_list_id, environment_list_id, working_directory_string_id, uid, gid, stdin_file_id, stdout_file_id, stderr_file_id, session_id)
end

export kill_process
"""
    $(SIGNATURES)

Sends a UNIX signal to a process object and returns the resulting error code.

Possible UNIX signals are:

* Interrupt = 2
* Quit = 3
* Abort = 6
* Kill = 9
* User1 = 10
* User2 = 12
* Terminate = 15
* Continue =  18
* Stop = 19
"""
function kill_process(device::BrickRED, process_id, signal)
    return device.deviceInternal.kill_process(process_id, signal)
end

export get_process_command
"""
    $(SIGNATURES)

Returns the executable, arguments, environment and working directory used to
spawn a process object, as passed to :func:`Spawn Process`, and the resulting
error code.
"""
function get_process_command(device::BrickRED, process_id, session_id)
    return device.deviceInternal.get_process_command(process_id, session_id)
end

export get_process_identity
"""
    $(SIGNATURES)

Returns the process ID and the user and group ID used to spawn a process object,
as passed to :func:`Spawn Process`, and the resulting error code.

The process ID is only valid if the state is *Running* or *Stopped*, see
:func:`Get Process State`.
"""
function get_process_identity(device::BrickRED, process_id)
    return device.deviceInternal.get_process_identity(process_id)
end

export get_process_stdio
"""
    $(SIGNATURES)

Returns the stdin, stdout and stderr files used to spawn a process object, as
passed to :func:`Spawn Process`, and the resulting error code.
"""
function get_process_stdio(device::BrickRED, process_id, session_id)
    return device.deviceInternal.get_process_stdio(process_id, session_id)
end

export get_process_state
"""
    $(SIGNATURES)

Returns the current state, timestamp and exit code of a process object, and
the resulting error code.

Possible process states are:

* Unknown = 0
* Running = 1
* Error = 2
* Exited = 3
* Killed = 4
* Stopped = 5

The timestamp represents the UNIX time since when the process is in its current
state.

The exit code is only valid if the state is *Error*, *Exited*, *Killed* or
*Stopped* and has different meanings depending on the state:

* Error: error code for error occurred while spawning the process (see below)
* Exited: exit status of the process
* Killed: UNIX signal number used to kill the process
* Stopped: UNIX signal number used to stop the process

Possible exit/error codes in *Error* state are:

* InternalError = 125
* CannotExecute = 126
* DoesNotExist = 127

The *CannotExecute* error can be caused by the executable being opened for
writing.
"""
function get_process_state(device::BrickRED, process_id)
    return device.deviceInternal.get_process_state(process_id)
end

export get_programs
"""
    $(SIGNATURES)

"""
function get_programs(device::BrickRED, session_id)
    return device.deviceInternal.get_programs(session_id)
end

export define_program
"""
    $(SIGNATURES)

"""
function define_program(device::BrickRED, identifier_string_id, session_id)
    return device.deviceInternal.define_program(identifier_string_id, session_id)
end

export purge_program
"""
    $(SIGNATURES)

"""
function purge_program(device::BrickRED, program_id, cookie)
    return device.deviceInternal.purge_program(program_id, cookie)
end

export get_program_identifier
"""
    $(SIGNATURES)

"""
function get_program_identifier(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_identifier(program_id, session_id)
end

export get_program_root_directory
"""
    $(SIGNATURES)

FIXME: root directory is absolute: <home>/programs/<identifier>
"""
function get_program_root_directory(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_root_directory(program_id, session_id)
end

export set_program_command
"""
    $(SIGNATURES)

FIXME: working directory is relative to <home>/programs/<identifier>/bin
"""
function set_program_command(device::BrickRED, program_id, executable_string_id, arguments_list_id, environment_list_id, working_directory_string_id)
    return device.deviceInternal.set_program_command(program_id, executable_string_id, arguments_list_id, environment_list_id, working_directory_string_id)
end

export get_program_command
"""
    $(SIGNATURES)

FIXME: working directory is relative to <home>/programs/<identifier>/bin
"""
function get_program_command(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_command(program_id, session_id)
end

export set_program_stdio_redirection
"""
    $(SIGNATURES)

FIXME: stdio file names are relative to <home>/programs/<identifier>/bin
"""
function set_program_stdio_redirection(device::BrickRED, program_id, stdin_redirection, stdin_file_name_string_id, stdout_redirection, stdout_file_name_string_id, stderr_redirection, stderr_file_name_string_id)
    return device.deviceInternal.set_program_stdio_redirection(program_id, stdin_redirection, stdin_file_name_string_id, stdout_redirection, stdout_file_name_string_id, stderr_redirection, stderr_file_name_string_id)
end

export get_program_stdio_redirection
"""
    $(SIGNATURES)

FIXME: stdio file names are relative to <home>/programs/<identifier>/bin
"""
function get_program_stdio_redirection(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_stdio_redirection(program_id, session_id)
end

export set_program_schedule
"""
    $(SIGNATURES)

"""
function set_program_schedule(device::BrickRED, program_id, start_mode, continue_after_error, start_interval, start_fields_string_id)
    return device.deviceInternal.set_program_schedule(program_id, start_mode, continue_after_error, start_interval, start_fields_string_id)
end

export get_program_schedule
"""
    $(SIGNATURES)

"""
function get_program_schedule(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_schedule(program_id, session_id)
end

export get_program_scheduler_state
"""
    $(SIGNATURES)

FIXME: message is currently valid in error-occurred state only
"""
function get_program_scheduler_state(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_program_scheduler_state(program_id, session_id)
end

export continue_program_schedule
"""
    $(SIGNATURES)

"""
function continue_program_schedule(device::BrickRED, program_id)
    return device.deviceInternal.continue_program_schedule(program_id)
end

export start_program
"""
    $(SIGNATURES)

"""
function start_program(device::BrickRED, program_id)
    return device.deviceInternal.start_program(program_id)
end

export get_last_spawned_program_process
"""
    $(SIGNATURES)

"""
function get_last_spawned_program_process(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_last_spawned_program_process(program_id, session_id)
end

export get_custom_program_option_names
"""
    $(SIGNATURES)

"""
function get_custom_program_option_names(device::BrickRED, program_id, session_id)
    return device.deviceInternal.get_custom_program_option_names(program_id, session_id)
end

export set_custom_program_option_value
"""
    $(SIGNATURES)

"""
function set_custom_program_option_value(device::BrickRED, program_id, name_string_id, value_string_id)
    return device.deviceInternal.set_custom_program_option_value(program_id, name_string_id, value_string_id)
end

export get_custom_program_option_value
"""
    $(SIGNATURES)

"""
function get_custom_program_option_value(device::BrickRED, program_id, name_string_id, session_id)
    return device.deviceInternal.get_custom_program_option_value(program_id, name_string_id, session_id)
end

export remove_custom_program_option
"""
    $(SIGNATURES)

"""
function remove_custom_program_option(device::BrickRED, program_id, name_string_id)
    return device.deviceInternal.remove_custom_program_option(program_id, name_string_id)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Brick is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position is the position in the stack from '0' (bottom) to '8' (top).

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickRED)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickRED, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
println(DEPOT_PATH)
# Local Singularity Image
# ["/home/koki/.julia", "/usr/local/julia/local/share/julia", "/usr/local/julia/share/julia"]

# Job Schedular (Slurm)
# ["/home/koki/Dev/mesh-workflow/.julia", "/usr/local/julia/local/share/julia", "/usr/local/julia/share/julia"]

println(LOAD_PATH)
# Local Singularity Image
# ["@", "@v#.#", "@stdlib"]

# Job Schedular (Slurm)
# ["@", "@v#.#", "@stdlib"]

println(Base.load_path())
# Local Singularity Image
# ["/home/koki/.julia/environments/v1.6/Project.toml", "/usr/local/julia/share/julia/stdlib/v1.6"]

# Job Schedular (Slurm)
# ["/home/koki/Dev/mesh-workflow/.julia/environments/v1.6/Project.toml", "/usr/local/julia/share/julia/stdlib/v1.6"]

using GZip
println(pathof(GZip))
# Local Singularity Image
# /home/koki/.julia/packages/GZip/JNmGn/src/GZip.jl
# Job Schedular (Slurm)
# ERROR: LoadError: ArgumentError: Package GZip not found in current

taxid = ARGS[1]
outfile = ARGS[2]

regex = "^" * "$taxid" * "\t" * "(.*?)" * "\t" * "(.*?)" * "\t" * "(.*?)" * "\t" * "(.*?)" * "\t" * "(.*?)" * "\t" * "(.*?)" * "\t.*"

open(outfile, "w") do fp
	GZip.open("data/gene/test.gz", "r") do fp2
		for line in eachline(fp2)
		  	m = match(Regex(regex), line)
			if !isnothing(m)
				if m[1] != "-" && m[6] != "-"
					write(fp, "$(m[1])\t$(m[6])\n")
				end
			end
		end
	end
end

ext: jl

max_stars_repo_name: manabuishii/mesh-workflow

max_issues_repo_name: manabuishii/mesh-workflow

max_forks_repo_name: manabuishii/mesh-workflow


text:
# script that wants to numerically test chatteris diaconis (misspelled with 99% prob)
# for the estimates of beta, fitness,ergm (many names..) in the DirBin1 case
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase#,DynNets
using PyCall; pygui(:qt); using PyPlot

Nsample = 100

# load my data
halfPeriod = false
## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

matY_T = YeMidWeekly_T
strIO_T = [sumSq(matY_T,2);sumSq(matY_T,1)]
N,~,T = size(matY_T)

t=1
(sum(strIO_T[1:N,t]) - sum(strIO_T[N+1:end,t]) )< 10^Base.eps()? S =sum(strIO_T[N+1:end,t]) :error()
SperLink  = S/(N*N-N)
tmp = -0.5log(SperLink)

1/exp(2tmp)
1/exp(tmp)
plot(strIO_T')
lstrI = log.(strIO_T[1:N,t])
lstrO = log.(strIO_T[1+N:end,t])
1./exp(-lstrI)
expMat = expMatrix(fooErgmDirW1Afixed,-[lstrI;lstrO],trues(N,N))

prevDegEqFlag = (degsIO_T[:,1:end-1] .== degsIO_T[:,2:end])
prevDegEqFlagI,prevDegEqFlagO = prevDegEqFlag[1:N,:], prevDegEqFlag[1+1N:end,:]
tmp = (sum(prevDegEqFlagI,1) - N).*(sum(prevDegEqFlagO,1) - N)/(N*(N-1))
plot(tmp')

prevDegDiff = abs.(degsIO_T[:,1:end-1] - degsIO_T[:,2:end])
sum(prevDegDiff.<=1)
degDiffMax = 20
plt[:hist](prevDegDiff[prevDegDiff.<=degDiffMax],degDiffMax,density = true);grid()
 title("Absolute Values of One Step Degree Variation ")
 ylabel("Fraction of Occurences")
# plot(1:degDiffMax,cumsum(counts(prevDegDiff[:]))[1:degDiffMax]/sum(counts(prevDegDiff[:])))

prevDegEqFlagI,prevDegEqFlagO = prevDegEqFlag[1:N,:], prevDegEqFlag[1+1N:end,:]

tmp = (sum(prevDegEqFlagI,1) - N).*(sum(prevDegEqFlagO,1) - N)/(N*(N-1))
plot(tmp')
allFitSS =  StaticNetsW.estimate( StaticNetsW.SnapSeqNetDirW1(strIO_T./10000000); identPost = false,identIter= true )
dgpPar = meanSq(allFitSS,2 )
expDgpMat = StaticNets.expMatrix(fooErgmDirBin1,dgpPar)
dgpDegs = [sumSq(expDgpMat,2); sumSq(expDgpMat,1)]
N = length(matA_T[1,:,1])
Nsample = 1000
L = maximum(abs.(dgpPar))
C = L
errBnd = C*sqrt(log(N)/N )
prob = 1 - C/(N^2)
dgpMod = ErgmDirBin1(zeros(dgpDegs),dgpPar,Int.(1:N))
findmax(dgpDegs)
#sample the dgp model
sam = sampl(dgpMod,Nsample )

est = zeros(2N,Nsample)
degsIO_Sample = [sumSq(sam,2);sumSq(sam,1)]
estPar = zeros(2N,Nsample)
estDegs = zeros(2N,Nsample)
#estimate
for i=1:Nsample
    println(i)
    estPar[:,i]= estimate(fooErgmDirBin1;degIO =degsIO_Sample[:,i])[1]
    expEstMat = StaticNets.expMatrix(fooErgmDirBin1,estPar[:,i])
    estDegs[:,i] = [sumSq(expEstMat,2); sumSq(expEstMat,1)]

end

estPar[.!isfinite(estPar)] = -5

diffdat = estDegs .- dgpDegs # abs.(estPar .- dgpPar)
plt[:hist](diffdat[:]  ,100,density=true )
 title("Sample  Static Model, Run MLE and Conpare Degrees")
 grid(which = "both" , linewidth=1,alpha = 0.5)
 xlim(-3, 3)
 xlabel("deg_dgp - deg_MLE")

sampMeanDegs = meanSq(degsIO_Sample,2)
ind1= findmin(sampMeanDegs)[2]
ind2 = findmax(sampMeanDegs)[2]
diffdat = abs.(estPar .- dgpPar)./abs.(dgpPar)
plt[:hist]([diffdat[ind1,:] diffdat[ind2,:]]  ,50,density=true )
 title("DGP fit = mean sing snap estimates  ")
 legend(["Sample Mean Deg = $(round(sampMeanDegs[ind1],2)) " ; "Sample Mean Deg = $(round(sampMeanDegs[ind2],2)) "])
  xlabel("θ_est - θ_dgp  ")
 grid(which = "both" , linewidth=1,alpha = 0.5)

degsIO_Sample[1,:]
degsIO_Sample[ind2,:]

funProb(x::Array{T,1} where T<: Real ) = 1./(1+exp.(-x))
derfunProb(x::Array{T,1} where T<: Real )  =  exp.(-x)./((1+exp.(-x)).^2)
x = Vector(-5:0.001:5)
othPar = -3:1.5:3
for i = 1:length(othPar) plot(x, funProb(x.+othPar[i] )); end
 legTex = ["θ_2 =  $(othPar[i])" for i = 1:length(othPar)]
 title( "1/(1+ exp( - θ - θ_2 ) )")
 legend(legTex)
 xlabel("θ")
 grid()

#


text:
@inline permute(i, N) = isodd(i) ? floor(Int, i/2) + 1 : N - floor(Int, (i-1)/2)

@inline unpermute(i, N) = i <= ceil(N/2) ? 2i-1 : 2(N-i+1)

function dct_makhoul_1d(A::CuArray)
    B = similar(A)
    N = length(A)

    for k in 1:N
        B[permute(k, N)] = A[k]
    end

    B = CUDA.CUFFT.fft(B)

    for k in 1:N
        B[k] = 2 * ω(4N, k-1) * B[k]
    end

    return real(B)
end

function idct_makhoul_1d(A::CuArray)
    B = similar(A, complex(eltype(A)))
    N = length(A)

    B[1] = 1/2 * ω(4N, 0) * A[1]
    for k in 2:N
        B[k] = ω(4N, 1-k) * A[k]
    end

    B = CUDA.CUFFT.ifft(B)

    C = similar(A)
    for k in 1:N
        C[unpermute(k, N)] = real(B[k])
    end

    return C
end

function dct_makhoul_2d(A::CuArray)
    Nx, Ny = size(A)

    # DCT along dimension 1

    B = similar(A)

    for j in 1:Ny, i in 1:Nx
        B[permute(i, Nx), j] = A[i, j]
    end

    B = CUDA.CUFFT.fft(B, 1)

    for j in 1:Ny, i in 1:Nx
        B[i, j] = 2 * ω(4Nx, i-1) * B[i, j]
    end

    B = real(B)

    # DCT along dimension 2

    C = similar(A)

    for j in 1:Ny, i in 1:Nx
        C[i, permute(j, Ny)] = B[i, j]
    end

    C = CUDA.CUFFT.fft(C, 2)

    for j in 1:Ny, i in 1:Nx
        C[i, j] = 2 * ω(4Ny, j-1) * C[i, j]
    end

    return real(C)
end

function idct_makhoul_2d(A::CuArray)
    Nx, Ny = size(A)

    # IDCT along dimension 1

    B = similar(A, complex(eltype(A)))

    for j in 1:Ny
        B[1, j] = 1/2 * ω(4Nx, 0) * A[1, j]
    end

    for j in 1:Ny, i in 2:Nx
        B[i, j] = ω(4Nx, 1-i) * A[i, j]
    end

    B = CUDA.CUFFT.ifft(B, 1)

    C = similar(A)
    for j in 1:Ny, i in 1:Nx
        C[unpermute(i, Nx), j] = real(B[i, j])
    end

    # IDCT along dimension 2

    D = similar(A, complex(eltype(A)))

    for i in 1:Nx
        D[i, 1] = 1/2 * ω(4Ny, 0) * C[i, 1]
    end

    for j in 2:Ny, i in 1:Nx
        D[i, j] = ω(4Ny, 1-j) * C[i, j]
    end

    D = CUDA.CUFFT.ifft(D, 2)

    E = similar(A)
    for j in 1:Ny, i in 1:Nx
        E[i, unpermute(j, Ny)] = real(D[i, j])
    end

    return E
end


text:
export Sin

"""
`Sin([domainType=Float64::Type,] dim_in::Tuple)`

Creates a sinusoid non-linear operator with input dimensions `dim_in`:
```math
\\sin( \\mathbf{x} ).
```

"""
struct Sin{T,N} <: NonLinearOperator
	dim::NTuple{N,Int}
end

function Sin(DomainType::Type, DomainDim::NTuple{N,Int}) where {N}
	Sin{DomainType,N}(DomainDim)
end

Sin(DomainDim::NTuple{N,Int}) where {N} = Sin{Float64,N}(DomainDim)
Sin(DomainDim::Vararg{Int}) = Sin{Float64,length(DomainDim)}(DomainDim)

function mul!(y::AbstractArray{T,N}, L::Sin{T,N}, x::AbstractArray{T,N}) where {T,N}
	y .= sin.(x)
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{A,TT}},
              b::AbstractArray) where {T,N, A<: Sin{T,N}, TT <: AbstractArray{T,N}}
    L = J.A
    y .= conj.(cos.(L.x)).*b
end

fun_name(L::Sin) = "sin"

size(L::Sin) = (L.dim, L.dim)

domainType(L::Sin{T,N}) where {T,N} = T
codomainType(L::Sin{T,N}) where {T,N} = T


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

typealias ContactElements3D Union{Tri3, Tri6, Quad4, Quad8, Quad9}

function create_orthogonal_basis(n)
    I = eye(3)
    k = indmax([norm(cross(n,I[:,k])) for k in 1:3])
    t1 = cross(n, I[:,k])/norm(cross(n, I[:,k]))
    t2 = cross(n, t1)
    return t1, t2
end

"""
Frictionless 2d small sliding contact.

problem
time
dimension
finite_sliding
friction
use_forwarddiff
"""
function assemble!(problem::Problem{Contact}, time::Float64,
                   ::Type{Val{2}}, ::Type{Val{false}},
                   ::Type{Val{false}}, ::Type{Val{false}}; debug=true)

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals = calculate_normals(slave_elements, time, Val{2};
                                rotate_normals=props.rotate_normals)
    update!(slave_elements, "normal", normals)

    # 2. loop all slave elements
    for (slave_num, slave_element) in enumerate(slave_elements)

        nsl = length(slave_element)
        X1 = slave_element("geometry", time)
        u1 = slave_element("displacement", time)
        la = slave_element("reaction force", time)
        n1 = slave_element("normal", time)
        t11, t21 = create_orthogonal_basis(n1[1])
        t12, t22 = create_orthogonal_basis(n1[2])
        t13, t23 = create_orthogonal_basis(n1[3])
        Q1_ = [n1[1] t11 t21]
        Q2_ = [n1[2] t12 t22]
        Q3_ = [n1[3] t13 t23]
        Z = zeros(3, 3)
        Q3 = [Q1_ Z Z; Z Q2_ Z; Z Z Q3_]
        contact_area = 0.0
        contact_error = 0.0

        element_area = 0.0
        for ip in get_integration_points(slave_element)
            detJ = slave_element(ip, time, Val{:detJ})
            w = ip.weight*detJ
            element_area += w
        end
        if "element area" in props.store_fields
            update!(slave_element, "element area", time => element_area)
        end

        if slave_num == 1
            info("First slave element area = $element_area")
            info("NT basis of first slave element")
            dump(Q3)
        end

        # project slave nodes to auxiliary plane (x0, Q)
        #xi = get_reference_element_midpoint(slave_element)
        xi = [1/3, 1/3]
        N = vec(get_basis(slave_element, xi, time))
        x0 = N*X1
        n0 = N*n1
        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]

        # 3. loop all master elements
        for master_element in slave_element("master elements", time)

            nm = length(master_element)
            X2 = master_element("geometry", time)
            u2 = master_element("displacement", time)
            x2 = X2 + u2

            #=
            norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
            norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
            norm(mean(X1) - X2[3]) / norm(X1[2] - X1[1]) < props.distval || continue
            =#

            # 3.1 project master nodes to auxiliary plane and create polygon clipping
            M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
            P = get_polygon_clip(S, M, n0)
            length(P) < 3 && continue # no clipping or shared edge (no volume)
            check_orientation!(P, n0)
            C0 = calculate_centroid(P)

            De = zeros(nsl, nsl)
            Me = zeros(nsl, nm)
            ge = zeros(field_dim*nsl)

            # 4. loop integration cells
            for cell in get_cells(P, C0)
                virtual_element = Element(Tri3)
                update!(virtual_element, "geometry", cell)

                # 5. loop integration point of integration cell
                for ip in get_integration_points(virtual_element, 3)

                    # project gauss point from auxiliary plane to master and slave element
                    x_gauss = virtual_element("geometry", ip, time)
                    if isnan(x_gauss[1])
                        info("is nan")
                        info("x_gauss = $x_gauss")
                        info("cell = $cell")
                        info("C0 = $C0")
                        info("P = $P")
                        info("S = $S")
                        info("M = $M")
                        info("n0 = $n0")
                        error("nan, unable to continue")
                    end
                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
                    xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)

                    detJ = virtual_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ
                    # add contributions
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    N2 = vec(get_basis(master_element, xi_m, time))
                    De += w*N1*N1'
                    Me += w*N1*N2'

                    x_s = N1*(X1+u1)
                    x_m = N2*(X2+u2)
                    ge += w*vec((x_m-x_s)*N1')
                    contact_area += w
                    n_s = N1*n1
                    contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
                end # integration points done

            end # integration cells done

            # 6. add contribution to contact virtual work
            sdofs = get_gdofs(problem, slave_element)
            mdofs = get_gdofs(problem, master_element)
            nsldofs = length(sdofs)
            nmdofs = length(mdofs)
            D3 = zeros(nsldofs, nsldofs)
            M3 = zeros(nmdofs, nmdofs)
            for i=1:field_dim
                D3[i:field_dim:end, i:field_dim:end] += De
                M3[i:field_dim:end, i:field_dim:end] += Me
            end

            add!(problem.assembly.C1, sdofs, sdofs, D3)
            add!(problem.assembly.C1, sdofs, mdofs, -M3)
            add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
            add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
            add!(problem.assembly.g, sdofs, Q3'*ge)

        end # master elements done

        if "contact area" in props.store_fields
            update!(slave_element, "contact area", time => contact_area)
        end

    end # slave elements done, contact virtual work ready

    S = sort(collect(keys(normals))) # slave element nodes
    weighted_gap = Dict{Int64, Vector{Float64}}()
    contact_pressure = Dict{Int64, Vector{Float64}}()
    complementarity_condition = Dict{Int64, Vector{Float64}}()
    is_active = Dict{Int64, Int}()
    is_inactive = Dict{Int64, Int}()
    is_slip = Dict{Int64, Int}()
    is_stick = Dict{Int64, Int}()

    g = full(problem.assembly.g)
    la = problem.assembly.la

    # active / inactive node detection
    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        weighted_gap[j] = g[dofs]
        if length(la) != 0
            normal = normals[j]
            tangent1, tangent2 = create_orthogonal_basis(normal)
            p = dot(normal, la[dofs])
            t1 = dot(tangent1, la[dofs])
            t2 = dot(tangent2, la[dofs])
            contact_pressure[j] = [p, t1, t2]
        else
            contact_pressure[j] = [0.0, 0.0, 0.0]
        end
        complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
        if complementarity_condition[j][1] < 0
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        else
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
        end
    end

    if "weighted gap" in props.store_fields
        update!(slave_elements, "weighted gap", time => weighted_gap)
    end
    if "contact pressure" in props.store_fields
        update!(slave_elements, "contact pressure", time => contact_pressure)
    end
    if "complementarity condition" in props.store_fields
        update!(slave_elements, "complementarity condition", time => complementarity_condition)
    end
    if "active nodes" in props.store_fields
        update!(slave_elements, "active nodes", time => is_active)
    end
    if "inactive nodes" in props.store_fields
        update!(slave_elements, "inactive nodes", time => is_inactive)
    end
    if "stick nodes" in props.store_fields
        update!(slave_elements, "stick nodes", time => is_stick)
    end
    if "slip nodes" in props.store_fields
        update!(slave_elements, "slip nodes", time => is_slip)
    end

    info("# | active | inactive | stick | slip | gap | pres | comp")
    for j in S
        str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
        str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
        info(str1 * str2)
    end

    # solve variational inequality

    C1 = sparse(problem.assembly.C1)
    ndofs = size(C1, 1)
    C2 = sparse(problem.assembly.C2)
    D = spzeros(ndofs, ndofs)

    # constitutive modelling in tangent direction, frictionless contact
    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        tdofs = dofs[[2,3]]
        if (is_active[j] == 1) && (is_slip[j] == 1)
            info("$j is in active/slip, removing tangential constraints $tdofs")
            C2[tdofs,:] = 0.0
            g[tdofs] = 0.0
            normal = normals[j]
            tangent1, tangent2 = create_orthogonal_basis(normal)
            D[tdofs[1], dofs] = tangent1
            D[tdofs[2], dofs] = tangent2
        end
    end

    # remove inactive nodes from assembly
    for j in S
        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
        if is_inactive[j] == 1
            info("$j is inactive, removing dofs $dofs")
            C1[dofs,:] = 0.0
            C2[dofs,:] = 0.0
            D[dofs,:] = 0.0
            g[dofs,:] = 0.0
        end
    end

    problem.assembly.C1 = C1
    problem.assembly.C2 = C2
    problem.assembly.D = D
    problem.assembly.g = g

end


text:
##
# LowRankFun represents f(x,y) by r column and row slices stored as Funs.
##

export LowRankFun

"""
`LowRankFun` gives an approximation to a bivariate function in low rank form.
"""
mutable struct LowRankFun{S<:Space,M<:Space,SS<:AbstractProductSpace,T<:Number} <: BivariateFun{T}
    A::Vector{VFun{S,T}}
    B::Vector{VFun{M,T}}
    space::SS

    function LowRankFun{S,M,SS,T}(A::Vector{VFun{S,T}},
                                  B::Vector{VFun{M,T}},
                                  space::SS) where {S,M,SS,T}
        @assert length(A) == length(B)
        @assert length(A) > 0
        new{S,M,SS,T}(A,B,space)
    end
end

LowRankFun(A::Vector{VFun{S,T}},B::Vector{VFun{M,T}},space::SS) where {S,M,SS,T} =
    LowRankFun{S,M,SS,T}(A,B,space)
LowRankFun(A::Vector{VFun{S,T}},B::Vector{VFun{M,T}}) where {S,M,T} =
    LowRankFun(A,B,space(first(A))⊗space(first(B)))
LowRankFun(A::Vector{VFun{S,T}},B::Vector{VFun{M,V}}) where {S,M,T,V} =
    LowRankFun(convert(Vector{VFun{S,promote_type(T,V)}},A),
               convert(Vector{VFun{M,promote_type(T,V)}},B),
               space(first(A))⊗space(first(B)))
rank(f::LowRankFun) = length(f.A)
size(f::LowRankFun,k::Integer) = k==1 ? mapreduce(length,max,f.A) : mapreduce(length,max,f.B)
size(f::LowRankFun) = size(f,1),size(f,2)

## Construction via a Matrix of coefficients

function LowRankFun(X::Array{T},dx::S,dy::M) where {S<:Space,M<:Space,T<:Number}
    U,Σ,V=svd(X)
    m=max(1,count(s->s>10eps(T),Σ))

    A=VFun{S,T}[Fun(dx,U[:,k].*sqrt(Σ[k])) for k=1:m]
    B=VFun{M,T}[Fun(dy,conj(V[:,k]).*sqrt(Σ[k])) for k=1:m]

    LowRankFun(A,B)
end

## Construction in a TensorSpace via a Vector of Funs

function LowRankFun(X::Vector{VFun{S,T}},d::TensorSpace{SV,DD}) where {S,T,DD<:EuclideanDomain{2},SV}
    @assert d[1] == space(X[1])
    LowRankFun(X,d[2])
end

function LowRankFun(X::Vector{VFun{S,T}},dy::Space) where {S,T}
    m=mapreduce(ncoefficients,max,X)
    M=zeros(T,m,length(X))
    for k=1:length(X)
        M[1:ncoefficients(X[k]),k]=X[k].coefficients
    end

    LowRankFun(M,space(X[1]),dy)
end

## Adaptive constructor selector

function LowRankFun(fin::Function,dx::Space,dy::Space;
                    method::Symbol=:standard,tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,
                    retmax::Bool=false,gridx::Integer=64,gridy::Integer=64,maxrank::Integer=100)
    f = dynamic(fin)
    if method == :standard
        F,maxabsf=standardLowRankFun(f,dx,dy;tolerance=tolerance,gridx=gridx,gridy=gridy,maxrank=maxrank)
    elseif method == :Cholesky
        @assert domain(dx) == domain(dy)
        if dx == dy
            F,maxabsf=CholeskyLowRankFun(f,dx;tolerance=tolerance,grid=max(gridx,gridy),maxrank=maxrank)
        else
            G,maxabsf=CholeskyLowRankFun(f,dx;tolerance=tolerance,grid=max(gridx,gridy),maxrank=maxrank)
            F=LowRankFun(copy(G.A),map(b->Fun(b,dy),G.B))
        end
    end
    retmax ? (F,maxabsf) : F
end

## Standard adaptive construction

function standardLowRankFun(f::Function,dx::Space,dy::Space;tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,gridx::Integer=64,gridy::Integer=64,maxrank::Integer=100)
    xy = checkpoints(dx⊗dy)
    T = promote_type(eltype(f(first(xy)...)),prectype(dx),prectype(dy))

    # We start by sampling on the given grid, find the approximate maximum and create the first rank-one approximation.
    ptsx,ptsy=points(dx,gridx),points(dy,gridy)
    X = zeros(T,gridx,gridy)
    maxabsf,r=findapproxmax!(f,X,ptsx,ptsy,gridx,gridy)
    if maxabsf < eps(zero(T))/eps(T) return LowRankFun([Fun(dx,[zero(T)])],[Fun(dy,[zero(T)])]),maxabsf end
    a,b=Fun(x->f(x,r[2]),dx),Fun(y->f(r[1],y),dy)

    # If necessary, we resize the grid to be at least as large as the
    # lengths of the first row and column Funs and we recompute the values of X.
    if gridx < ncoefficients(a) || gridy < ncoefficients(b)
        gridx,gridy = max(gridx,ncoefficients(a)),max(gridy,ncoefficients(b))
        ptsx,ptsy=points(dx,gridx),points(dy,gridy)
        X = zeros(T,gridx,gridy)
        maxabsf,r=findapproxmax!(f,X,ptsx,ptsy,gridx,gridy)
        a,b=Fun(x->f(x,r[2]),dx),Fun(y->f(r[1],y),dy)
    end

    A,B=typeof(a)[],typeof(b)[]
    if tolerance == :relative
        tol = 100maxabsf*eps(T)
    elseif tolerance[1] == :absolute
        tol = 100*tolerance[2]*eps(T)
    end
    tol10 = tol/10
    Avals,Bvals = zeros(T,gridx),zeros(T,gridy)
    p₁,p₂ = plan_transform(dx,Avals),plan_transform(dy,Bvals)

    # Eat, drink, subtract rank-one, repeat.
    for k=1:maxrank
        if (norm(a.coefficients,Inf) < tol || norm(b.coefficients,Inf) < tol)
            return LowRankFun(A,B),maxabsf
        end
        A,B =[A;a/sqrt(abs(a(r[1])))],[B;b/(sqrt(abs(b(r[2])))*sign(b(r[2])))]
        r=findapproxmax!(A[k],B[k],X,ptsx,ptsy,gridx,gridy)
        Ar,Br=evaluate(A,r[1]),evaluate(B,r[2])
        for i=1:gridx
            @inbounds Avals[i] = f(ptsx[i],r[2])
        end
        for j=1:gridy
            @inbounds Bvals[j] = f(r[1],ptsy[j])
        end
        a,b = Fun(dx,p₁*Avals) - dotu(Br,A),Fun(dy,p₂*Bvals) - dotu(Ar,B)
        chop!(a,tol10),chop!(b,tol10)
    end
    @warn "Maximum rank of " * string(maxrank) * " reached"
    return LowRankFun(A,B),maxabsf
end

## Adaptive Cholesky decomposition, when f is Hermitian positive (negative) definite

function CholeskyLowRankFun(f::Function,dx::Space;tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,grid::Integer=64,maxrank::Integer=100)
    xy = checkpoints(dx⊗dx)
    T = promote_type(eltype(f(first(xy)...)),prectype(dx))

    # We start by sampling on the given grid, find the approximate maximum and create the first rank-one approximation.
    pts=points(dx,grid)
    X = zeros(T,grid)
    maxabsf,r=findcholeskyapproxmax!(f,X,pts,grid)
    if maxabsf < eps(zero(T))/eps(T) return LowRankFun([Fun(dx,[zero(T)])],[Fun(dx,[zero(T)])]),maxabsf end
    a=Fun(x->f(x,r),dx)

    # If necessary, we resize the grid to be at least as large as the
    # ncoefficients of the first row/column Fun and we recompute the values of X.
    if grid < ncoefficients(a)
        grid = max(grid,ncoefficients(a))
        pts=points(dx,grid)
        X = zeros(T,grid)
        maxabsf,r=findcholeskyapproxmax!(f,X,pts,grid)
        a=Fun(x->f(x,r),dx)
    end

    A,B=typeof(a)[],typeof(a)[]
    if tolerance == :relative
        tol = 100maxabsf*eps(T)
    elseif tolerance[1] == :absolute
        tol = 100*tolerance[2]*eps(T)
    end
    tol10 = tol/10
    Avals = zeros(T,grid)
    p₁ = plan_transform(dx,Avals)

    # Eat, drink, subtract rank-one, repeat.
    for k=1:maxrank
        if norm(a.coefficients,Inf) < tol return LowRankFun(A,B),maxabsf end
        A,B = [A;a/sqrt(abs(a(r)))],[B;a/(sqrt(abs(a(r)))*sign(a(r)))]
        r=findcholeskyapproxmax!(A[k],B[k],X,pts,grid)
        Br=evaluate(B,r)
        for i=1:grid
            @inbounds Avals[i] = f(pts[i],r)
        end
        a = Fun(dx,p₁*Avals) - dotu(Br,A)
        chop!(a,tol10)
    end
    @warn "Maximum rank of " * string(maxrank) * " reached"
    return LowRankFun(A,B),maxabsf
end

## Construction via TensorSpaces and ProductDomains

LowRankFun(f::Function,S::TensorSpace{SV,DD,RR};kwds...) where {SV,DD<:EuclideanDomain{2},RR} =
    LowRankFun(dynamic(f),factor(S,1),factor(S,2);kwds...)
LowRankFun(f::Function,dx::Domain,dy::Domain;kwds...) =
    LowRankFun(dynamic(f),Space(dx),Space(dy);kwds...)
LowRankFun(f::Function,d::ProductDomain;kwds...) =
    LowRankFun(dynamic(f),factors(d)...;kwds...)

LowRankFun(f::Function;kwds...) = LowRankFun(dynamic(f),ChebyshevInterval(),ChebyshevInterval();kwds...)

## Construction from values

LowRankFun(A::Array{T}) where {T<:Number} = LowRankFun(A,Interval{T}(),Interval{T}())
LowRankFun(c::Number,etc...) = LowRankFun((x,y)->c,etc...)

## Construction from other LowRankFuns

LowRankFun(f::LowRankFun,d1::IntervalOrSegment,d2::IntervalOrSegment) =
    LowRankFun(map(g->Fun(d1,g.coefficients),f.A),
               map(g->Fun(d2,g.coefficients),f.B))
LowRankFun(f::LowRankFun) = LowRankFun(f,ChebyshevInterval(),ChebyshevInterval())

## Utilities

function findapproxmax!(f::Function,X::AbstractMatrix,ptsx::AbstractVector,ptsy::Vector,gridx,gridy)
    for j=1:gridy
        ptsyj = ptsy[j]
        @simd for k=1:gridx
            @inbounds X[k,j]+=f(ptsx[k],ptsyj)
        end
    end
    maxabsf,impt = findmaxabs(X)
    imptple = CartesianIndices((gridx,gridy))[impt]
    maxabsf,[ptsx[imptple[1]],ptsy[imptple[2]]]
end

function findapproxmax!(A::Fun,B::Fun,X::AbstractMatrix,ptsx::Vector,ptsy::Vector,gridx,gridy)
    Ax,By = A.(ptsx),B.(ptsy)
    subtractrankone!(Ax,By,X,gridx,gridy)
    maxabsf,impt = findmaxabs(X)
    imptple = CartesianIndices((gridx,gridy))[impt]
    [ptsx[imptple[1]],ptsy[imptple[2]]]
end

function findcholeskyapproxmax!(f::Function,X::AbstractVector,pts::Vector,grid)
    @simd for k=1:grid
        @inbounds X[k]+=f(pts[k],pts[k])
    end
    maxabsf,impt = findmaxabs(X)
    maxabsf,pts[impt]
end

function findcholeskyapproxmax!(A::Fun,B::Fun,X::Vector,pts::Vector,grid)
    Ax,By = A.(pts),B.(pts)
    subtractrankone!(Ax,By,X,grid)
    maxabsf,impt = findmaxabs(X)
    pts[impt]
end

function subtractrankone!(A::AbstractVector,B::AbstractVector,X::AbstractMatrix,gridx::Int,gridy::Int)
    for j=1:gridy
        @inbounds Bj = B[j]
        @simd for k=1:gridx
            @inbounds X[k,j] -= A[k]*Bj
        end
    end
end

function subtractrankone!(A::AbstractVector,B::AbstractVector,X::AbstractVector,grid::Int)
    @simd for k=1:grid
        @inbounds X[k] -= A[k]*B[k]
    end
end

## TODO: in Julia base?
function findmaxabs(a)
    if isempty(a)
        throw(ArgumentError("collection must be non-empty"))
    end
    m = abs(a[1])
    mi = 1
    for i in eachindex(a)
        ai = abs(a[i])
        if ai > m || m!=m
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

(f::LowRankFun)(x,y)=evaluate(f,x,y)

domain(f::LowRankFun,k::Integer) = k==1 ? domain(first(f.A)) : domain(first(f.B))
space(f::LowRankFun,k::Integer) = k==1 ? space(first(f.A)) : space(first(f.B))
space(f::LowRankFun)=f.space

transpose(f::LowRankFun{S,M,SS,T}) where {S,M,SS,T}=LowRankFun(f.B,f.A,transpose(space(f)))

function values(f::LowRankFun)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=values(pad(f.A[k],xm))*transpose(values(pad(f.B[k],ym)))
    end
    ret
end

#TODO: this is inconsistent with 1D where it does canonical
function coefficients(f::LowRankFun)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=pad(f.A[k].coefficients,xm)*transpose(pad(f.B[k].coefficients,ym))
    end
    ret
end

function coefficients(f::LowRankFun,n::Space,m::Space)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=pad(coefficients(f.A[k],n),xm)*transpose(pad(coefficients(f.B[k],m),ym))
    end
    ret
end

function vecpoints(f::LowRankFun,k::Integer)
    if k==1
        xm=mapreduce(ncoefficients,max,f.A)
        points(space(first(f.A)),xm)
    else
        ym=mapreduce(ncoefficients,max,f.B)
        points(space(first(f.B)),ym)
    end
end

evaluate(A::Vector{T},B::Vector{M},x,y) where {T<:Fun,M<:Fun}=dotu(evaluate(A,x),evaluate(B,y))
evaluate(A::Vector{T},B::Vector{M},x::AbstractVector,y::AbstractVector) where {T<:Fun,M<:Fun}=evaluate.(transpose(A),x)*evaluate.(B,transpose(y))

evaluate(f::LowRankFun,x,y)=evaluate(f.A,f.B,x,y)
evaluate(f::LowRankFun,::Colon,::Colon)=f
evaluate(f::LowRankFun,x::Number,::Colon)=dotu(f.B,evaluate(f.A,x))

function evaluate(f::LowRankFun,::Colon,y::Number)
    m = maximum(map(ncoefficients,f.A))
    r=rank(f)
    ret = zeros(m)

    for k=1:r
        for j=1:ncoefficients(f.A[k])
            @inbounds ret[j] += f.A[k].coefficients[j]*f.B[k](y)
        end
    end

    Fun(first(f.A).space,ret)
end

"""
    f.(x::AbstractVector, y':::AbstractMatrix)
Fast evaluation of a LowRankFun on a cartesian grid x ⨂ y.
"""
function broadcast(f::LowRankFun, x::AbstractVector, y::AbstractMatrix)
    ret = zeros(length(y), length(x))
    #println(InRed * "FROM BROADCAST" * InDefault)
    for k = 1:rank(f)
        ret +=  (f.A[k].(x) * f.B[k].(y))'
    end
    ret
end

## Truncate
#TODO: should reduce rank if needed
chop(f::LowRankFun,tol)=LowRankFun(map(g->chop(g,tol),f.A),map(g->chop(g,tol),f.B),f.space)
function pad(f::LowRankFun,m::Integer,n::Integer)
    A,B = deepcopy(f.A),deepcopy(f.B)

    for k=1:rank(f)
        pad!(A[k],m);pad!(B[k],n)
    end

    LowRankFun(A,B,f.space)
end

## Algebra

for op = (:*,:/)
    @eval ($op)(A::Array{T,1},c::Number) where {T<:Fun}=map(f->($op)(f,c),A)
    @eval ($op)(f::LowRankFun,c::Number) = LowRankFun(($op)(f.A,c),f.B)
    @eval ($op)(c::Number,f::LowRankFun) = LowRankFun(($op)(c,f.A),f.B)
end

# Let K be a LowRankFun and f be a Fun.
# op(f,K) acts as operating in the x variable, and
# op(K,f) acts as operating in the y variable.

for op = (:*,:/)
    @eval ($op)(f::Fun,K::LowRankFun) = LowRankFun(($op).(f,K.A),K.B)
    @eval ($op)(K::LowRankFun,f::Fun) = LowRankFun(K.A,($op).(K.B,f))
end

+(f::LowRankFun,g::LowRankFun) = LowRankFun([f.A;g.A],[f.B;g.B])
-(f::LowRankFun) = LowRankFun(-f.A,f.B)
-(f::LowRankFun,g::LowRankFun) = f+(-g)

## QR factorization of a LowRankFun

function qr(f::LowRankFun)
    sp,r = space(f),rank(f)
    Q,R = qr(coefficients(f.A))
    BR = coefficients(f.B)*transpose(R)
    LowRankFun(map(i->Fun(sp[1],Q[:,i]),1:r),map(i->Fun(sp[2],BR[:,i]),1:r),sp)
end

## Special functions

real(u::LowRankFun)=LowRankFun([map(real,u.A),map(imag,u.A)],[map(real,u.B),-map(imag,u.B)])
imag(u::LowRankFun)=LowRankFun([map(real,u.A),map(imag,u.A)],[map(imag,u.B),map(real,u.B)])

## Calculus

sum(g::LowRankFun)=dotu(map(sum,g.A),map(sum,g.B))
sum(g::LowRankFun,n::Integer)=(n==1) ? dotu(g.B,map(sum,g.A)) : dotu(g.A,map(sum,g.B))
cumsum(g::LowRankFun,n::Integer)=(n==1) ? LowRankFun(map(cumsum,g.A),copy(g.B)) : LowRankFun(copy(g.A),map(cumsum,g.B))
differentiate(g::LowRankFun,n::Integer)=(n==1) ? LowRankFun(map(differentiate,g.A),copy(g.B)) : LowRankFun(copy(g.A),map(differentiate,g.B))
integrate(g::LowRankFun,n::Integer)=(n==1) ? LowRankFun(map(integrate,g.A),copy(g.B)) : LowRankFun(copy(g.A),map(integrate,g.B))


text:
mutable struct KrigingTuner
    counter::Int
    g_inter::Int
    l_inter::Int
    g_opt::KrigingModel.GlobalOptimizer
    l_opt::KrigingModel.LocalOptimizer
    kernbounds
    meanbounds
end

function KrigingTuner(kernbounds,
    meanbounds;
    global_every::Integer = 50,
    global_opt::KrigingModel.GlobalOptimizer = ISRESOptimizer(maxeval = 10000),
    local_every::Integer = 5,
    local_opt::KrigingModel.LocalOptimizer = LBFGSOptimizer(maxeval = 1000))
    global_every < 0 && throw(ArgumentError("`global_every` must be non-negative"))
    local_every < 0 && throw(ArgumentError("`local_every` must be non-negative"))
    if local_every==global_every==0
        @warn "Both `global_every` and `local_every` are zeros, never tune hyperparameters!"
    end
    if (local_every != 0 && global_every % local_every != 0)
        throw(ArgumentError("`global_every` must be a multiple of `local_every`"))
    end
    KrigingTuner(0, global_every, local_every, global_opt, local_opt,kernbounds, meanbounds)
end

function tunekrg!(krg::Kriging, tuner::KrigingTuner)
    nobs = length(krg.gps[1].y)
    if tuner.g_inter > 0 && ceil(nobs/tuner.g_inter) > ceil(tuner.counter/tuner.g_inter)
        tune!(krg, tuner.kernbounds, tuner.meanbounds, tuner.g_opt)
    end
    if tuner.l_inter > 0 && ceil(nobs/tuner.l_inter) > ceil(tuner.counter/tuner.l_inter)
        tune!(krg, tuner.kernbounds, tuner.meanbounds, tuner.l_opt)
    end
    tuner.counter = nobs
end

function forcetune!(krg::Kriging, tuner::KrigingTuner)
    tune!(krg, tuner.kernbounds, tuner.meanbounds, tuner.g_opt)
    tune!(krg, tuner.kernbounds, tuner.meanbounds, tuner.l_opt)
    tuner.counter = length(krg.gps[1].y)
end


text:
module KLDivergence
	export kl_divergence, kl_sample_divergence, skl_divergence, skl_sample_divergence, js_divergence, js_sample_divergence, mc_divergence_estimate

	include("kl.jl")
end


text:
"Holds the tableau of an explicit Runge-Kutta method."
struct TableauERK{T} <: AbstractTableauERK{T}
    @HeaderTableau

    q::CoefficientsRK{T}

    function TableauERK{T}(q) where {T}
        @assert q.c[1] == 0
        @assert istrilstrict(q.a)
        @assert !(q.s==1 && q.a[1,1] ≠ 0)
        new(q.name, q.o, q.s, q)
    end
end

function TableauERK(q::CoefficientsRK{T}) where {T}
    TableauERK{T}(q)
end

function TableauERK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    TableauERK{T}(CoefficientsRK(name, order, a, b, c))
end

"Read explicit Runge-Kutta tableau from file."
function readTableauERKFromFile(dir::AbstractString, name::AbstractString)
    file = string(dir, "/", name, ".tsv")

    o, s, T = readTableauRKHeaderFromFile(file)

    # TODO Read data in original format (e.g., Rational).
    #      For this we need to save tableaus as jld or hdf5.
#    tab_array = readdlm(file, T)
    tab_array = readdlm(file, comments=true)

    if s == 0
        s = size(tab_array, 1)
    end

    @assert s == size(tab_array, 1) == size(tab_array, 2)-2

    c = tab_array[1:s, 1]
    b = tab_array[1:s, 2]
    a = tab_array[1:s, 3:s+2]

    get_config(:verbosity) > 1 ? @info("Reading explicit Runge-Kutta tableau $(name) with $(s) stages and order $(o) from file\n$(file)") : nothing

    TableauERK(Symbol(name), o, a, b, c)
end

"Explicit Runge-Kutta integrator."
struct IntegratorERK{DT,TT,FT} <: DeterministicIntegrator{DT,TT}
    equation::ODE{DT,TT,FT}
    tableau::TableauERK{TT}
    Δt::TT

    function IntegratorERK{DT,TT,FT}(equation, tableau, Δt) where {DT,TT,FT}
        new(equation, tableau, Δt)
    end
end

function IntegratorERK(equation::ODE{DT,TT,FT}, tableau::TableauERK{TT}, Δt::TT) where {DT,TT,FT}
    IntegratorERK{DT,TT,FT}(equation, tableau, Δt)
end

equation(int::IntegratorERK) = int.equation
timestep(int::IntegratorERK) = int.Δt

"Explicit Runge-Kutta integrator cache."
mutable struct IntegratorCacheERK{DT,TT,D,S} <: ODEIntegratorCache{DT,D}
    n::Int
    t::TT
    t̅::TT
    q::Vector{TwicePrecision{DT}}
    q̅::Vector{TwicePrecision{DT}}
    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}

    function IntegratorCacheERK{DT,TT,D,S}() where {DT,TT,D,S}
        q = zeros(TwicePrecision{DT}, D)
        q̅ = zeros(TwicePrecision{DT}, D)
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        new(0, zero(TT), zero(TT), q, q̅, Q, V)
    end
end

function create_integrator_cache(int::IntegratorERK{DT,TT}) where {DT,TT}
    IntegratorCacheERK{DT, TT, ndims(int), int.tableau.s}()
end

function CommonFunctions.reset!(cache::IntegratorCacheERK{DT,TT}, Δt::TT) where {DT,TT}
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.set_solution!(cache::IntegratorCacheERK, sol, n=0)
    t, q = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
end

"Integrate ODE with explicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorERK{DT,TT}, cache::IntegratorCacheERK{DT,TT}) where {DT,TT}
    local tᵢ::TT
    local yᵢ::DT

    # reset cache
    reset!(cache, int.Δt)

    # compute internal stages
    for i in 1:int.tableau.q.s
        @inbounds for k in eachindex(cache.Q[i], cache.V[i])
            yᵢ = 0
            for j in 1:i-1
                yᵢ += int.tableau.q.a[i,j] * cache.V[j][k]
            end
            cache.Q[i][k] = cache.q̅[k] + int.Δt * yᵢ
        end
        tᵢ = cache.t̅ + int.Δt * int.tableau.q.c[i]
        int.equation.v(tᵢ, cache.Q[i], cache.V[i])
    end

    # compute final update
    update_solution!(cache.q, cache.V, int.tableau.q.b, int.Δt)

    # take care of periodic solutions
    cut_periodic_solution!(cache.q, int.equation.periodicity)
end


text:
# Pre-process the dataset before use in discriminator training
using CSV
using ProgressMeter

"""
    pad_seq(batch::AbstractVector{T}, pad_token::Integer=0)

Function to add pad tokens in shorter sequence, to make the length of each sequence equal to the `max_length` ( calculated as `max(map(length, batch))`) in the batch. Pad token defaults to `0`.
"""
function pad_seq(batch::AbstractVector{T}, pad_token::Integer=0) where T	# send a list of sequence
    length_vec = map(length, batch)
    max_len = maximum(length_vec)
    padded_seq = zeros(max_len, length(batch))
    for (index,tokens) in enumerate(batch)
        len = length_vec[index]
        padded_seq[1:len, index] = tokens
        padded_seq[len+1:end, index] .= pad_token
    end
    return padded_seq
end

"""
    get_mask(seq::AbstractMatrix{T}, pad_token::Integer=0, embed_size::Integer=768)

Function to create mask for sequences against padding, so as to inform the model, that some part of sequenece is padded and hence to be ignored.
"""
function get_mask(seq::AbstractMatrix{T}, pad_token::Integer=0, embed_size::Integer=768) where T
    seq_len, batch_size = size(seq)
    mask = reshape(repeat(seq.!=pad_token, embed_size),  seq_len, embed_size, batch_size)
    permutedims(mask, (2,1,3))	# shape (embed_size, seq_length, batch_size)
end

# expect a tuple of (X, Y) where each is a list
"""
    data_preprocess(data_x, data_y, classification_type::String="Binary", num_classes::Integer=2; args=nothing)

Function to preprocess `data_x` and `data_y` along with creating mask for the data_x.

Preprocessing for `data_x` consist of padding with pad token (expected to be provided as `args.pad_token`).

Preprocessing for `data_y` consist of creating `onehotbach` for `data_y` (if `classification_type` is not "Binary"), for `1:num_classes` else reshape the data as `(1, length(data_y))`

Returns `data_x`, `data_y`, `mask` after pre-processing.

"""
function data_preprocess(data_x, data_y, classification_type::String="Binary", num_classes::Integer=2; args=nothing)
    if classification_type!="Binary"
        data_y = Float32.(onehotbatch(data_y, 1:num_classes))
    else
        data_y = Float32.(reshape(data_y, 1, length(data_y)))
    end
    if args==nothing
        pad_token = 0
        embed_size = 768
    else
        pad_token = args.pad_token
        embed_size = args.embed_size
    end
    data_x = pad_seq(data_x, pad_token)
    mask_src = get_mask(data_x, pad_token, embed_size)
    return Int.(data_x), data_y, Float32.(mask_src)
end

"""
    truncate_(x, max_length::Integer)

Truncate the data to minimum of `max_length` and length of `x`.

"""
function truncate_(x, max_length::Integer)
    x[1:min(max_length, length(x))]
end

"""
    load_data(data_x, data_y, tokenizer::PretrainedTokenizer;  batchsize::Integer=8, truncate::Bool=false, max_length::Integer=256, shuffle::Bool=false, drop_last::Bool=false, add_eos_start::Bool=true)

Returns DataLoader for the `data_x` and `data_y` after processing the data_x, with batchsize=`batchsize`. The processing consist of tokenization of data_x and further truncation to `max_len` if `truncate` is set to be true.

If `add_eos_start` is set to true, add EOS token of tokenizer to the start.

"""
function load_data(data_x, data_y, tokenizer::PretrainedTokenizer; batchsize::Integer=8, truncate::Bool=false, max_length::Integer=256, shuffle::Bool=false, drop_last::Bool=false, add_eos_start::Bool=true)
    data_x = map(tokenizer, data_x)
    if add_eos_start==true
        data_x = map(pushfirst!, data_x, repeat([tokenizer.eos_token_id], length(data_x)))
    end
    if truncate==true
        data_x = truncate_.(data_x, max_length)
    end
    loader = DataLoader((data_x, data_y), batchsize=batchsize, partial=drop_last, shuffle=shuffle)
    return loader
end

"""
    load_cached_data(discrim::Union{DiscriminatorV1, DiscriminatorV2}, data_x, data_y, tokenizer::PretrainedTokenizer; truncate::Bool=false, max_length::Integer=256, shuffle::Bool=false, batchsize::Int=4, drop_last::Bool=false, classification_type="Binary", num_classes=2, args=nothing)

Returns a DataLoader with (x, y) which can directly be feeded into classifier layer for training.

The function first loads the data using [`load_data`] (@ref PPLM.load_data) function with batchsize=1, then passes each batch to the transformer model of `discrim` after data preprocessing, and then the average representation of the `hidden_states` are stored in a vector, which are then further loaded into a DataLoader, ready to use for classification training.

**Note**: This functions saves time by cacheing the average representation of hidden states beforehand, avoiding passing the data through model in each epoch of training. This can be done as the model itself is `non-trainable` while training discriminator classifier head.

"""
function load_cached_data(discrim::Union{DiscriminatorV1, DiscriminatorV2}, data_x, data_y, tokenizer::PretrainedTokenizer; truncate::Bool=false, max_length::Integer=256, shuffle::Bool=false, batchsize::Int=4, drop_last::Bool=false, classification_type="Binary", num_classes=2, args=nothing)
    global device
    # if label zero indexed, increase by 1
    if 0 in data_y
        data_y.+=1
    end
    loader = load_data(data_x, data_y, tokenizer, batchsize=1,
                truncate=truncate, max_length=max_length, shuffle=false, drop_last=drop_last)

    xs = Array{Float32}(undef, discrim.embed_size, 0)

    if classification_type == "Binary"
        ys = Array{Int}(undef, 1, 0)
    else
        ys = Array{Float32}(undef, num_classes, 0)
    end

    @showprogress "Computing..." for (x, y) in loader
        x_, y_, mask =  (data_preprocess(x, y, classification_type, num_classes; args=args))|> device
        if typeof(discrim)==DiscriminatorV1
            data_x_ = avg_representation(discrim, x_)
        else
            data_x_ = avg_representation(discrim, x_, mask)
        end
        data_x_ = Float32.(data_x_) |> cpu
        data_y_ = Float32.(y_) |> cpu
        xs = cat(xs, data_x_, dims=2)
        ys = cat(ys, data_y_, dims=2)
    end
    if args!=nothing
        batchsize=args.batchsize
    end
    loader = DataLoader((xs, ys), batchsize=batchsize, partial=drop_last, shuffle=shuffle)
    return loader
end

"""
    load_data_from_csv(path_to_csv; text_col="text", label_col="label", delim=',', header=1)

Load the data from a csv file based on the specified `text_col` column for text and `label_col` for target label. Returns vectors for `text` and `label`.
"""
function load_data_from_csv(path_to_csv; text_col="text", label_col="label", delim=',', header=1)
    df = CSV.File(path_to_csv; delim=delim) |> DataFrame
    return df[!, text_col], df[!, label_col]
end


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

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

# Everything in Julia is an expression.

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

# All of the normal infix operators are available.
1 + 1      # => 2
8 - 1      # => 7
10 * 2     # => 20
35 / 5     # => 7.0
5 / 2      # => 2.5  # dividing integers always results in a Float64
div(5, 2)  # => 2    # for a truncated result, use div
5 \ 35     # => 7.0
2^2        # => 4    # power, not bitwise xor
12 % 10    # => 2

# Enforce precedence with parentheses
(1 + 3) * 2  # => 8

# Bitwise Operators
~2         # => -3 # bitwise not
3 & 5      # => 1  # bitwise and
2 | 4      # => 6  # bitwise or
xor(2, 4)  # => 6  # bitwise xor
2 >>> 1    # => 1  # logical shift right
2 >> 1     # => 1  # arithmetic shift right
2 << 1     # => 4  # logical/arithmetic shift left

# Use the bitstring function to see the binary representation of a number.
bitstring(12345)
# => "0000000000000000000000000000000000000000000000000011000000111001"
bitstring(12345.0)
# => "0100000011001000000111001000000000000000000000000000000000000000"

# Boolean values are primitives
true
false

# Boolean operators
!true   # => false
!false  # => true
1 == 1  # => true
2 == 1  # => false
1 != 1  # => false
2 != 1  # => true
1 < 10  # => true
1 > 10  # => false
2 <= 2  # => true
2 >= 2  # => true
# Comparisons can be chained
1 < 2 < 3  # => true
2 < 3 < 2  # => false

# Strings are created with "
"This is a string."

# Character literals are written with '
'a'

# Strings are UTF8 encoded. Only if they contain only ASCII characters can
# they be safely indexed.
ascii("This is a string")[1]  # => 'T' # Julia indexes from 1
# Otherwise, iterating over strings is recommended (map, for loops, etc).

# $ can be used for string interpolation:
"2 + 2 = $(2 + 2)" # => "2 + 2 = 4"
# You can put any Julia expression inside the parentheses.

# Another way to format strings is the printf macro from the stdlib Printf.
using Printf
@printf "%d is less than %f\n" 4.5 5.3  # => 5 is less than 5.300000

# Printing is easy
println("I'm Julia. Nice to meet you!")

# String can be compared lexicographically
"good" > "bye" # => true
"good" == "good" # => true
"1 + 2 = 3" == "1 + 2 = $(1 + 2)" # => true

####################################################
## 2. Variables and Collections
####################################################

# You don't declare variables before assigning to them.
some_var = 5  # => 5
some_var  # => 5

# Accessing a previously unassigned variable is an error
try
    some_other_var  # => ERROR: UndefVarError: some_other_var not defined
catch e
    println(e)
end

# Variable names start with a letter or underscore.
# After that, you can use letters, digits, underscores, and exclamation points.
SomeOtherVar123! = 6  # => 6

# You can also use certain unicode characters
☃ = 8  # => 8
# These are especially handy for mathematical notation
2 * π # => 6.283185307179586

# A note on naming conventions in Julia:
#
# * Word separation can be indicated by underscores ('_'), but use of
#   underscores is discouraged unless the name would be hard to read
#   otherwise.
#
# * Names of Types begin with a capital letter and word separation is shown
#   with CamelCase instead of underscores.
#
# * Names of functions and macros are in lower case, without underscores.
#
# * Functions that modify their inputs have names that end in !. These
#   functions are sometimes called mutating functions or in-place functions.

# Arrays store a sequence of values indexed by integers 1 through n:
a = Int64[]  # => 0-element Int64 Array

# 1-dimensional array literals can be written with comma-separated values.
b = [4, 5, 6]  # => 3-element Int64 Array: [4, 5, 6]
b = [4; 5; 6]  # => 3-element Int64 Array: [4, 5, 6]
b[1]  # => 4
b[end]  # => 6

# 2-dimensional arrays use space-separated values and semicolon-separated rows.
matrix = [1 2; 3 4]  # => 2x2 Int64 Array: [1 2; 3 4]

# Arrays of a particular type
b = Int8[4, 5, 6]  # => 3-element Int8 Array: [4, 5, 6]

# Add stuff to the end of a list with push! and append!
push!(a, 1)    # => [1]
push!(a, 2)    # => [1,2]
push!(a, 4)    # => [1,2,4]
push!(a, 3)    # => [1,2,4,3]
append!(a, b)  # => [1,2,4,3,4,5,6]

# Remove from the end with pop
pop!(b)  # => 6 and b is now [4,5]

# Let's put it back
push!(b, 6)  # b is now [4,5,6] again.

a[1]  # => 1  # remember that Julia indexes from 1, not 0!

# end is a shorthand for the last index. It can be used in any
# indexing expression
a[end]  # => 6

# we also have popfirst! and pushfirst!
popfirst!(a)  # => 1 and a is now [2,4,3,4,5,6]
pushfirst!(a, 7)  # => [7,2,4,3,4,5,6]

# Function names that end in exclamations points indicate that they modify
# their argument.
arr = [5,4,6]  # => 3-element Int64 Array: [5,4,6]
sort(arr)  # => [4,5,6]; arr is still [5,4,6]
sort!(arr)  # => [4,5,6]; arr is now [4,5,6]

# Looking out of bounds is a BoundsError
try
    a[0]
    # => BoundsError: attempt to access 7-element Array{Int64,1} at index [0]
    a[end + 1]
    # => BoundsError: attempt to access 7-element Array{Int64,1} at index [8]
catch e
    println(e)
end

# Errors list the line and file they came from, even if it's in the standard
# library. You can look in the folder share/julia inside the julia folder to
# find these files.

# You can initialize arrays from ranges
a = [1:5;]  # => 5-element Int64 Array: [1,2,3,4,5]

# You can look at ranges with slice syntax.
a[1:3]  # => [1, 2, 3]
a[2:end]  # => [2, 3, 4, 5]

# Remove elements from an array by index with splice!
arr = [3,4,5]
splice!(arr, 2)  # => 4 ; arr is now [3,5]

# Concatenate lists with append!
b = [1,2,3]
append!(a, b)  # Now a is [1, 2, 3, 4, 5, 1, 2, 3]

# Check for existence in a list with in
in(1, a)  # => true

# Examine the length with length
length(a)  # => 8

# Tuples are immutable.
tup = (1, 2, 3)  # => (1,2,3)  # an (Int64,Int64,Int64) tuple.
tup[1]  # => 1
try
    tup[1] = 3  # => ERROR: no method setindex!((Int64,Int64,Int64),Int64,Int64)
catch e
    println(e)
end

# Many array functions also work on tuples
length(tup)  # => 3
tup[1:2]  # => (1,2)
in(2, tup)  # => true

# You can unpack tuples into variables
a, b, c = (1, 2, 3)  # => (1,2,3)  # a is now 1, b is now 2 and c is now 3

# Tuples are created even if you leave out the parentheses
d, e, f = 4, 5, 6  # => (4,5,6)

# A 1-element tuple is distinct from the value it contains
(1,) == 1  # => false
(1) == 1  # => true

# Look how easy it is to swap two values
e, d = d, e  # => (5,4)  # d is now 5 and e is now 4

# Dictionaries store mappings
empty_dict = Dict()  # => Dict{Any,Any}()

# You can create a dictionary using a literal
filled_dict = Dict("one" => 1, "two" => 2, "three" => 3)
# => Dict{String,Int64}

# Look up values with []
filled_dict["one"]  # => 1

# Get all keys
keys(filled_dict)
# => Base.KeySet for a Dict{String,Int64} with 3 entries. Keys:
# "two", "one", "three"
# Note - dictionary keys are not sorted or in the order you inserted them.

# Get all values
values(filled_dict)
# => Base.ValueIterator{Dict{String,Int64}} with 3 entries. Values: 2, 1, 3
# Note - Same as above regarding key ordering.

# Check for existence of keys in a dictionary with in, haskey
in(("one" => 1), filled_dict)  # => true
in(("two" => 3), filled_dict)  # => false
haskey(filled_dict, "one")     # => true
haskey(filled_dict, 1)         # => false

# Trying to look up a non-existent key will raise an error
try
    filled_dict["four"]  # => KeyError: key "four" not found
catch e
    println(e)
end

# Use the get method to avoid that error by providing a default value
# get(dictionary, key, default_value)
get(filled_dict, "one", 4)   # => 1
get(filled_dict, "four", 4)  # => 4

# Use Sets to represent collections of unordered, unique values
empty_set = Set()  # => Set{Any}()
# Initialize a set with values
filled_set = Set([1, 2, 2, 3, 4])  # => Set([4, 2, 3, 1])

# Add more values to a set
push!(filled_set, 5)  # => Set([4, 2, 3, 5, 1])

# Check if the values are in the set
in(2, filled_set)  # => true
in(10, filled_set)  # => false

# There are functions for set intersection, union, and difference.
other_set = Set([3, 4, 5, 6])          # => Set([4, 3, 5, 6])
intersect(filled_set, other_set)       # => Set([4, 3, 5])
union(filled_set, other_set)           # => Set([4, 2, 3, 5, 6, 1])
setdiff(Set([1,2,3,4]), Set([2,3,5]))  # => Set([4, 1])

####################################################
## 3. Control Flow
####################################################

# Let's make a variable
some_var = 5

# Here is an if statement. Indentation is not meaningful in Julia.
if some_var > 10
    println("some_var is totally bigger than 10.")
elseif some_var < 10    # This elseif clause is optional.
    println("some_var is smaller than 10.")
else                    # The else clause is optional too.
    println("some_var is indeed 10.")
end
# => prints "some var is smaller than 10"

# For loops iterate over iterables.
# Iterable types include Range, Array, Set, Dict, and AbstractString.
for animal = ["dog", "cat", "mouse"]
    println("$animal is a mammal")
    # You can use $ to interpolate variables or expression into strings
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

# You can use 'in' instead of '='.
for animal in ["dog", "cat", "mouse"]
    println("$animal is a mammal")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

for pair in Dict("dog" => "mammal", "cat" => "mammal", "mouse" => "mammal")
    from, to = pair
    println("$from is a $to")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

for (k, v) in Dict("dog" => "mammal", "cat" => "mammal", "mouse" => "mammal")
    println("$k is a $v")
end
# prints:
#    dog is a mammal
#    cat is a mammal
#    mouse is a mammal

# While loops loop while a condition is true
let x = 0
    while x < 4
        println(x)
        x += 1  # Shorthand for x = x + 1
    end
end
# prints:
#   0
#   1
#   2
#   3

for i = 1:10
    z = i
end

z

x = 0;

for i = 1:10
    local x
    x = i + 1
end

x

for i = 1:10
    global z
    z = i
end

z

z=0

for i = 1:10
    z = i
end

z

# Handle exceptions with a try/catch block
try
    error("help")
catch e
    println("caught it $e")
end
# => caught it ErrorException("help")

####################################################
## 4. Functions
####################################################

# The keyword 'function' creates new functions
#function name(arglist)
#  body...
#end
function add(x, y)
    println("x is $x and y is $y")

    # Functions return the value of their last statement
    x + y
end

add(5, 6)  # => 11 after printing out "x is 5 and y is 6"

# Compact assignment of functions
f_add(x, y) = x + y  # => "f (generic function with 1 method)"
f_add(3, 4)  # => 7

# Function can also return multiple values as tuple
fn(x, y) = x + y, x - y
fn(3, 4)  # => (7, -1)

# You can define functions that take a variable number of
# positional arguments
function varargs(args...)
    return args
    # use the keyword return to return anywhere in the function
end
# => varargs (generic function with 1 method)

varargs(1, 2, 3)  # => (1,2,3)

# The ... is called a splat.
# We just used it in a function definition.
# It can also be used in a function call,
# where it will splat an Array or Tuple's contents into the argument list.
add([5,6]...)  # this is equivalent to add(5,6)

x = (5, 6)  # => (5,6)
add(x...)  # this is equivalent to add(5,6)

# You can define functions with optional positional arguments
function defaults(a, b, x=5, y=6)
    return "$a $b and $x $y"
end

defaults('h', 'g')  # => "h g and 5 6"
defaults('h', 'g', 'j')  # => "h g and j 6"
defaults('h', 'g', 'j', 'k')  # => "h g and j k"
try
    defaults('h')  # => ERROR: no method defaults(Char,)
    defaults()  # => ERROR: no methods defaults()
catch e
    println(e)
end

# You can define functions that take keyword arguments
function keyword_args(;k1=4, name2="hello")  # note the ;
    return Dict("k1" => k1, "name2" => name2)
end

keyword_args(name2="ness")  # => ["name2"=>"ness","k1"=>4]
keyword_args(k1="mine")  # => ["k1"=>"mine","name2"=>"hello"]
keyword_args()  # => ["name2"=>"hello","k1"=>4]

# You can combine all kinds of arguments in the same function
function all_the_args(normal_arg, optional_positional_arg=2; keyword_arg="foo")
    println("normal arg: $normal_arg")
    println("optional arg: $optional_positional_arg")
    println("keyword arg: $keyword_arg")
end

all_the_args(1, 3, keyword_arg=4)
# prints:
#   normal arg: 1
#   optional arg: 3
#   keyword arg: 4

# Julia has first class functions
function create_adder(x)
    adder = function (y)
        return x + y
    end
    return adder
end

# This is "stabby lambda syntax" for creating anonymous functions
(x -> x > 2)(3)  # => true

# This function is identical to create_adder implementation above.
function create_adder(x)
    y -> x + y
end

# You can also name the internal function, if you want
function create_adder(x)
    function adder(y)
        x + y
    end
    adder
end

add_10 = create_adder(10)
add_10(3)  # => 13

# There are built-in higher order functions
map(add_10, [1,2,3])  # => [11, 12, 13]
filter(x -> x > 5, [3, 4, 5, 6, 7])  # => [6, 7]

# We can use list comprehensions
[add_10(i) for i = [1, 2, 3]] # => [11, 12, 13]
[add_10(i) for i in [1, 2, 3]] # => [11, 12, 13]
[x for x in [3, 4, 5, 6, 7] if x > 5] # => [6, 7]

####################################################
## 5. Types
####################################################

# Julia has a type system.
# Every value has a type; variables do not have types themselves.
# You can use the `typeof` function to get the type of a value.
typeof(5)  # => Int64

# Types are first-class values
typeof(Int64)  # => DataType
typeof(DataType)  # => DataType
# DataType is the type that represents types, including itself.

# Types are used for documentation, optimizations, and dispatch.
# They are not statically checked.

# Users can define types
# They are like records or structs in other languages.
# New types are defined using the `struct` keyword.

# struct Name
#   field::OptionalType
#   ...
# end
struct Tiger
    taillength::Float64
    coatcolor  # not including a type annotation is the same as `::Any`
end

# The default constructor's arguments are the properties
# of the type, in the order they are listed in the definition
tigger = Tiger(3.5, "orange")  # => Tiger(3.5,"orange")

# The type doubles as the constructor function for values of that type
sherekhan = typeof(tigger)(5.6, "fire")  # => Tiger(5.6,"fire")

# These struct-style types are called concrete types
# They can be instantiated, but cannot have subtypes.
# The other kind of types is abstract types.

# abstract Name
abstract type Cat end  # just a name and point in the type hierarchy

# Abstract types cannot be instantiated, but can have subtypes.
using InteractiveUtils  # defines the subtype and supertype function
# For example, Number is an abstract type
subtypes(Number)  # => 2-element Array{Any,1}:
                 #     Complex{T<:Real}
                 #     Real
subtypes(Cat)  # => 0-element Array{Any,1}

# AbstractString, as the name implies, is also an abstract type
subtypes(AbstractString)  # 4-element Array{Any,1}:
                            #  String
                            #  SubString
                            #  SubstitutionString
                            #  Test.GenericString

# Every type has a super type; use the `supertype` function to get it.
typeof(5)  # => Int64
supertype(Int64)  # => Signed
supertype(Signed)  # => Integer
supertype(Integer)  # => Real
supertype(Real)  # => Number
supertype(Number)  # => Any
supertype(supertype(Signed))  # => Real
supertype(Any)  # => Any
# All of these type, except for Int64, are abstract.
typeof("fire")  # => String
supertype(String)  # => AbstractString
# Likewise here with String
supertype(SubString)  # => AbstractString

# <: is the subtyping operator
struct Lion <: Cat  # Lion is a subtype of Cat
    mane_color
    roar::AbstractString
end

# You can define more constructors for your type
# Just define a function of the same name as the type
# and call an existing constructor to get a value of the correct type
Lion(roar::AbstractString) = Lion("green", roar)
# This is an outer constructor because it's outside the type definition

struct Panther <: Cat  # Panther is also a subtype of Cat
    eye_color
    Panther() = new("green")
    # Panthers will only have this constructor, and no default constructor.
end
# Using inner constructors, like Panther does, gives you control
# over how values of the type can be created.
# When possible, you should use outer constructors rather than inner ones.

####################################################
## 6. Multiple-Dispatch
####################################################

# In Julia, all named functions are generic functions
# This means that they are built up from many small methods
# Each constructor for Lion is a method of the generic function Lion.

# For a non-constructor example, let's make a function meow:

# Definitions for Lion, Panther, Tiger
function meow(animal::Lion)
    animal.roar  # access type properties using dot notation
end

function meow(animal::Panther)
    "grrr"
end

function meow(animal::Tiger)
    "rawwwr"
end

# Testing the meow function
meow(tigger)  # => "rawwr"
meow(Lion("brown", "ROAAR"))  # => "ROAAR"
meow(Panther())  # => "grrr"

# Review the local type hierarchy
Tiger <: Cat # => false
Lion <: Cat # => true
Panther <: Cat # => true

# Defining a function that takes Cats
function pet_cat(cat::Cat)
    println("The cat says $(meow(cat))")
end

pet_cat(Lion("42"))  # => prints "The cat says 42"
try
    pet_cat(tigger)  # => ERROR: no method pet_cat(Tiger,)
catch e
    println(e)
end

# In OO languages, single dispatch is common;
# this means that the method is picked based on the type of the first argument.
# In Julia, all of the argument types contribute to selecting the best method.

# Let's define a function with more arguments, so we can see the difference
function fight(t::Tiger, c::Cat)
    println("The $(t.coatcolor) tiger wins!")
end
# => fight (generic function with 1 method)

fight(tigger, Panther())  # => prints The orange tiger wins!
fight(tigger, Lion("ROAR"))  # => prints The orange tiger wins!

# Let's change the behavior when the Cat is specifically a Lion
fight(t::Tiger, l::Lion) = println("The $(l.mane_color)-maned lion wins!")
# => fight (generic function with 2 methods)

fight(tigger, Panther())  # => prints The orange tiger wins!
fight(tigger, Lion("ROAR"))  # => prints The green-maned lion wins!

# We don't need a Tiger in order to fight
fight(l::Lion, c::Cat) = println("The victorious cat says $(meow(c))")
# => fight (generic function with 3 methods)

fight(Lion("balooga!"), Panther())  # => prints The victorious cat says grrr
try
    fight(Panther(), Lion("RAWR"))
catch e
    println(e)
    # => MethodError(fight, (Panther("green"), Lion("green", "RAWR")),
    #                0x000000000000557b)
end

# Also let the cat go first
fight(c::Cat, l::Lion) = println("The cat beats the Lion")

# This warning is because it's unclear which fight will be called in:
try
    fight(Lion("RAR"), Lion("brown", "rarrr"))
    # => prints The victorious cat says rarrr
catch e
    println(e)
  # => MethodError(fight, (Lion("green", "RAR"), Lion("brown", "rarrr")),
  #                0x000000000000557c)
end
# The result may be different in other versions of Julia

fight(l::Lion, l2::Lion) = println("The lions come to a tie")
fight(Lion("RAR"), Lion("brown", "rarrr"))  # => prints The lions come to a tie

# Under the hood
# You can take a look at the llvm  and the assembly code generated.

square_area(l) = l * l  # square_area (generic function with 1 method)

square_area(5)  # => 25

# What happens when we feed square_area an integer?
code_native(square_area, (Int32,))
    #        .section    __TEXT,__text,regular,pure_instructions
    #    Filename: none
    #    Source line: 1              # Prologue
    #        push    RBP
    #        mov RBP, RSP
    #    Source line: 1
    #        movsxd  RAX, EDI        # Fetch l from memory?
    #        imul    RAX, RAX        # Square l and store the result in RAX
    #        pop RBP                 # Restore old base pointer
    #        ret                     # Result will still be in RAX

code_native(square_area, (Float32,))
    #        .section    __TEXT,__text,regular,pure_instructions
    #    Filename: none
    #    Source line: 1
    #        push    RBP
    #        mov RBP, RSP
    #    Source line: 1
    #        vmulss  XMM0, XMM0, XMM0  # Scalar single precision multiply (AVX)
    #        pop RBP
    #        ret

code_native(square_area, (Float64,))
    #        .section    __TEXT,__text,regular,pure_instructions
    #    Filename: none
    #    Source line: 1
    #        push    RBP
    #        mov RBP, RSP
    #    Source line: 1
    #        vmulsd  XMM0, XMM0, XMM0  # Scalar double precision multiply (AVX)
    #        pop RBP
    #        ret
    #
# Note that julia will use floating point instructions if any of the
# arguments are floats.
# Let's calculate the area of a circle
circle_area(r) = pi * r * r     # circle_area (generic function with 1 method)
circle_area(5)  # 78.53981633974483

code_native(circle_area, (Int32,))
    #        .section    __TEXT,__text,regular,pure_instructions
    #    Filename: none
    #    Source line: 1
    #        push    RBP
    #        mov RBP, RSP
    #    Source line: 1
    #        vcvtsi2sd   XMM0, XMM0, EDI          # Load integer (r) from memory
    #        movabs  RAX, 4593140240              # Load pi
    #        vmulsd  XMM1, XMM0, QWORD PTR [RAX]  # pi * r
    #        vmulsd  XMM0, XMM0, XMM1             # (pi * r) * r
    #        pop RBP
    #        ret
    #

code_native(circle_area, (Float64,))
    #        .section    __TEXT,__text,regular,pure_instructions
    #    Filename: none
    #    Source line: 1
    #        push    RBP
    #        mov RBP, RSP
    #        movabs  RAX, 4593140496
    #    Source line: 1
    #        vmulsd  XMM1, XMM0, QWORD PTR [RAX]
    #        vmulsd  XMM0, XMM1, XMM0
    #        pop RBP
    #        ret
    #