text:
using Test,SphericalHarmonicModes
import SphericalHarmonicModes: NonContiguousError, ModeMissingError, InvalidModeError,
NegativeDegreeError, tRangeError, SHModeRange

@testset "constructors" begin

	@testset "LM" begin
		@test LM(0:3,-1:1) == LM(0,3,-1,1)
		@test LM(0:3,1) == LM(0,3,1,1)
		@test LM(2,1:2) == LM(2,2,1,2)
		@test LM(0:3) == LM(0,3,-3,3)
		@test LM(2) == LM(2,2,-2,2)
	end

	@testset "ML" begin
		@test ML(0:3,-1:1) == ML(0,3,-1,1)
		@test ML(0:3,1) == ML(0,3,1,1)
		@test ML(2,1:2) == ML(2,2,1,2)
		@test ML(0:3) == ML(0,3,-3,3)
		@test ML(2) == ML(2,2,-2,2)
	end

	@testset "L₂L₁ " begin
		Δs_max = rand(1:3)
		l_range = rand(1:3):rand(4:10)
		@test L₂L₁(l_range,Δs_max) == L₂L₁(l_range,LM(0:Δs_max,0))
		@test L₂L₁(l_range,Δs_max) == L₂L₁(l_range,ML(0:Δs_max,0))
	end
end

@testset "length" begin
	s_cutoff=5
	@testset "ML" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
			mr = ML(smin,smax,tmin,tmax)
			@test begin
				res = length(mr) == sum(length(m_range(mr,l))
								for l in l_range(mr))
				if !res
					println(mr)
				end
				res
			end
		end
	end

	@testset "LM" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
			mr = LM(smin,smax,tmin,tmax)
			@test begin
				res = length(mr) == sum(length(l_range(mr,m))
								for m in m_range(mr))
				if !res
					println(mr)
				end
				res
			end
		end
	end

	@testset "LM==ML" begin
	    for smin=0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
			m1 = LM(smin,smax,tmin,tmax)
			m2 = ML(smin,smax,tmin,tmax)
			@test length(m1) == length(m2)
		end
	end

	@testset "L₂L₁ default s′minmax" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,Δs_max=0:s_cutoff
			mr = L₂L₁(smin,smax,Δs_max)
			@test length(mr) == sum(length(l₂_range(mr,l₁)) for l₁ in l₁_range(mr))
		end
	end

	@testset "L₂L₁ all" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,Δs_max=0:s_cutoff
			mr = L₂L₁(smin,smax,Δs_max)
			for spmin=0:smax+Δs_max,spmax=spmin:smax+Δs_max
				m = L₂L₁(smin,smax,Δs_max,spmin,spmax)
				@test begin
					res = length(mr) == sum(length(l₂_range(mr,l₁)) for l₁ in l₁_range(mr))
					if !res
						println(mr)
					end
					res
				end
			end
		end
	end
end

@testset "LM ML modes" begin
	s_cutoff = 1
	for smin in 0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
		l_range = smin:smax
		m_range = tmin:tmax
		m1 = LM(l_range,m_range)
		m2 = ML(l_range,m_range)
		@test sort(collect(m1)) == sort(collect(m2))
	end
end

@testset "modeindex" begin

	function modeindex2(m::ML,s::Integer,t::Integer)
		N_skip = 0
		for si in m.l_min:s-1
			N_skip += length(m_range(m,si))
		end

		N_skip + searchsortedfirst(m_range(m,s),t)
	end

	function modeindex2(m::LM,s::Integer,t::Integer)
		N_skip = 0
		for ti in m.m_min:t-1
			N_skip += length(l_range(m,ti))
		end

		N_skip + searchsortedfirst(l_range(m,t),s)
	end

	function modeindex2(m::L₂L₁,s′::Integer,s::Integer)
		N_skip = 0
		for si in m.l₁_min:s-1
			N_skip += length(l₂_range(m,si))
		end

		N_skip + searchsortedfirst(l₂_range(m,s),s′)
	end

	modeindex2(m::SHModeRange,(s,t)::Tuple) = modeindex(m,s,t)

	s_cutoff = 1
	@testset "LM" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
			m1 = LM(smin,smax,tmin,tmax)
			for (s,t) in m1
				@test modeindex(m1,s,t) == modeindex2(m1,s,t)
			end
			m1c = collect(m1)
			for t in m_range(m1), s1 in l_range(m1,t),
					s2 in l_range(m1,t)

				smin,smax=minmax(s1,s2)
				@test modeindex(m1,smin:smax,t) ==
				findfirst(isequal((smin,t)),m1c):findfirst(isequal((smax,t)),m1c)
			end
		end
	end

	@testset "ML" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,tmin=-smax:smax,tmax=tmin:smax
			m2 = ML(smin,smax,tmin,tmax)
			for (s,t) in m2
				@test modeindex(m2,s,t) == modeindex2(m2,s,t)
			end
			m2c = collect(m2)
			for s in l_range(m2), t1 in  m_range(m2,s),
					t2 in m_range(m2,s)

				tmin,tmax = minmax(t1,t2)
				@test modeindex(m2,s,tmin:tmax) ==
				findfirst(isequal((s,tmin)),m2c):findfirst(isequal((s,tmax)),m2c)
			end
		end
	end

	@testset "L₂L₁" begin
		for smin=0:s_cutoff,smax=smin:s_cutoff,
			Δs_max=0:s_cutoff,s′min=0:s_cutoff,s′max=s′min:s_cutoff

			m3 = L₂L₁(smin,smax,Δs_max,s′min,s′max)
			for (s′,s) in m3
				@test begin
					res = modeindex(m3,s′,s) == modeindex2(m3,s′,s)
					if !res
						println(m3)
					end
					res
				end
			end
			m3c = collect(m3)
			for s in l₁_range(m3), s′1 in l₂_range(m3,s),
				s′2 in l₂_range(m3,s)

				s′min,s′max = minmax(s′1,s′2)
				@test begin
					res = modeindex(m3,s′min:s′max,s) ==
					findfirst(isequal((s′min,s)),m3c):findfirst(isequal((s′max,s)),m3c)
					if !res
						println(m3)
					end
					res
				end
			end
		end
	end

	@testset "ranges" begin
	    @testset "LM" begin
	        m=LM(0:3,-2:1)
	        sr = 1:3; tr = 0:1
	        ind1 = modeindex(m,minimum(sr),minimum(tr))
	        ind2 = modeindex(m,maximum(sr),maximum(tr))
	        @test modeindex(m,sr,tr) == ind1:ind2
	        @test_throws ModeMissingError modeindex(m,0:5,0:2)
	        @test_throws NonContiguousError modeindex(m,1:3,-1:0)
	        @test_throws InvalidModeError modeindex(m,l_range(m),m_range(m))
	    end
	    @testset "ML" begin
	        m=ML(0:5,-2:1)
	        sr = 2:5; tr = -2:1
	        ind1 = modeindex(m,minimum(sr),minimum(tr))
	        ind2 = modeindex(m,maximum(sr),maximum(tr))
	        @test modeindex(m,sr,tr) == ind1:ind2
	        @test_throws ModeMissingError modeindex(m,0:5,0:2)
	        @test_throws NonContiguousError modeindex(m,4:5,0:1)
	        @test_throws InvalidModeError modeindex(m,l_range(m),m_range(m))
	    end
	    @testset "L₂L₁" begin
	        m=L₂L₁(0:3,2,0:3);
	        s′r = 0:2; sr = 0:1;
	        ind1 = modeindex(m,minimum(s′r),minimum(sr))
	        ind2 = modeindex(m,maximum(s′r),maximum(sr))
	        @test modeindex(m,s′r,sr) == ind1:ind2
	        @test_throws ModeMissingError modeindex(m,0:4,0:4)
	        @test_throws NonContiguousError modeindex(m,1:3,1:2)
	        @test_throws NonContiguousError modeindex(m,l₂_range(m),l₁_range(m))
	    end
	end

	@testset "ModeRange" begin
	    @testset "LM" begin
	        m=LM(0:2)
	        # Non-Continguous
	        mpart=LM(0:1,0:1)
	        @test modeindex(m,mpart) == 4:7
	        # Contiguous
	        mpart=LM(0:2,0:0)
	        @test modeindex(m,mpart) == 4:6
	        @test collect(m)[modeindex(m,mpart)] == collect(mpart)
	    end
	    @testset "ML" begin
	        m=ML(0:2,-1:1)
	        # Non-Continguous
	        mpart=ML(1:2,0:1)
	        @test modeindex(m,mpart) == 3:7
	        # Contiguous
	        mpart=ML(1:2,-1:1)
	        @test modeindex(m,mpart) == 2:7
	        @test collect(m)[modeindex(m,mpart)] == collect(mpart)
	    end
	    @testset "L₂L₁" begin
	        m=L₂L₁(0:2,2);
	        # Non-Continguous
	        mpart=L₂L₁(1:2,1,1:2)
	        @test modeindex(m,mpart) == 5:10
	        # Contiguous
	        mpart=L₂L₁(1:2,2)
	        @test modeindex(m,mpart) == 4:12
	        @test collect(m)[modeindex(m,mpart)] == collect(mpart)
	    end
	end
end

@testset "last" begin
	@testset "ML" begin
		m1 = ML(rand(1:5),rand(6:10))
		@test last(collect(m1)) == last(m1)
	end

	@testset "LM" begin
		m2 = LM(rand(1:5),rand(6:10))
		@test last(collect(m2)) == last(m2)
	end

	@testset "L₂L₁" begin
		m3 = L₂L₁(rand(1:3):rand(4:10),rand(1:5))
		@test last(collect(m3)) == last(m3)
	end
end


text:
using Plots

rotation(x,y, theta)= [cos(theta) -sin(theta); sin(theta) cos(theta)]*[x;y]
points=[[0.0;0.0],[1.0;0.0]]
for step in 1:10
    i=1
    while points[i]!=points[end]
        insert!(points,i+1,(points[i+1]-points[i])/3 + points[i])
        tempo=((points[i+2]-points[i])*(2/3)+points[i])-points[i+1]
        point=rotation(tempo[1], tempo[2], pi/3)+points[i+1]
        insert!(points,i+2,point)
        insert!(points,i+3,tempo+points[i+1])
        i+=4
    end
end

plot(hcat(points...)[1,:] , hcat(points...)[2,:],xlabel="x", ylabel="y")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter2\\report\\Fig\\Koch.png")

ext: jl


text:
function bracket_sign_change(f′, a, b; k=2)
    if a > b; a,b = b,a; end # ensure a < b
    center, half_width = (b+a)/2, (b-a)/2
    while f′(a)*f′(b) > 0
        half_width *= k
        a = center - half_width
        b = center + half_width
    end
    return (a,b)
end

function bisection(f′, a, b, ϵ)
    if a > b; a,b = b,a; end # ensure a < b

    ya, yb = f′(a), f′(b)
    if ya == 0; b = a; end
    if yb == 0; a = b; end

    while b - a > ϵ
        x = (a+b)/2
        y = f′(x)
        if y == 0
            a, b = x, x
        elseif sign(y) == sign(ya)
            a = x
        else
            b = x
        end
    end
    return (a,b)
end

f = x->sin(10x)+cos(3x)
df(x) = 10*cos(10*x)-3*sin(3*x)

as,bs = bracket_sign_change(df,3.4,3.8)
as,df(as),bs,df(bs)

af,bf = bisection(df,as,bs,1.0e-8)

df(af),f(af)


text:
export BrickletBarometerV2AirPressureCallbackConfiguration
struct BrickletBarometerV2AirPressureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletBarometerV2AltitudeCallbackConfiguration
struct BrickletBarometerV2AltitudeCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletBarometerV2TemperatureCallbackConfiguration
struct BrickletBarometerV2TemperatureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletBarometerV2MovingAverageConfiguration
struct BrickletBarometerV2MovingAverageConfiguration
    moving_average_length_air_pressure::Integer
    moving_average_length_temperature::Integer
end

export BrickletBarometerV2Calibration
struct BrickletBarometerV2Calibration
    measured_air_pressure::Integer
    actual_air_pressure::Integer
end

export BrickletBarometerV2SensorConfiguration
struct BrickletBarometerV2SensorConfiguration
    data_rate::Integer
    air_pressure_low_pass_filter::Integer
end

export BrickletBarometerV2SPITFPErrorCount
struct BrickletBarometerV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

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

export BrickletBarometerV2
"""
Measures air pressure and altitude changes
"""
mutable struct BrickletBarometerV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_air_pressure
"""
    $(SIGNATURES)

Returns the measured air pressure.

If you want to get the value periodically, it is recommended to use the
:cb:`Air Pressure` callback. You can set the callback configuration
with :func:`Set Air Pressure Callback Configuration`.
"""
function get_air_pressure(device::BrickletBarometerV2)
    return device.deviceInternal.get_air_pressure()
end

export set_air_pressure_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Air Pressure` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Air Pressure` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_air_pressure_callback_configuration(device::BrickletBarometerV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_air_pressure_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_air_pressure_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Air Pressure Callback Configuration`.
"""
function get_air_pressure_callback_configuration(device::BrickletBarometerV2)
    return device.deviceInternal.get_air_pressure_callback_configuration()
end

export get_altitude
"""
    $(SIGNATURES)

Returns the relative altitude of the air pressure sensor. The value
is calculated based on the difference between the
current air pressure and the reference air pressure that can be set
with :func:`Set Reference Air Pressure`.

If you want to get the value periodically, it is recommended to use the
:cb:`Altitude` callback. You can set the callback configuration
with :func:`Set Altitude Callback Configuration`.
"""
function get_altitude(device::BrickletBarometerV2)
    return device.deviceInternal.get_altitude()
end

export set_altitude_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Altitude` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Altitude` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_altitude_callback_configuration(device::BrickletBarometerV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_altitude_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_altitude_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Altitude Callback Configuration`.
"""
function get_altitude_callback_configuration(device::BrickletBarometerV2)
    return device.deviceInternal.get_altitude_callback_configuration()
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the air pressure sensor.

This temperature is used internally for temperature compensation
of the air pressure measurement. It is not as accurate as the
temperature measured by the :ref:`temperature_v2_bricklet` or the
:ref:`temperature_ir_v2_bricklet`.

If you want to get the value periodically, it is recommended to use the
:cb:`Temperature` callback. You can set the callback configuration
with :func:`Set Temperature Callback Configuration`.
"""
function get_temperature(device::BrickletBarometerV2)
    return device.deviceInternal.get_temperature()
end

export set_temperature_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Temperature` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_temperature_callback_configuration(device::BrickletBarometerV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_temperature_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_temperature_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Temperature Callback Configuration`.
"""
function get_temperature_callback_configuration(device::BrickletBarometerV2)
    return device.deviceInternal.get_temperature_callback_configuration()
end

export set_moving_average_configuration
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the air pressure and temperature measurements.

Setting the length to 1 will turn the averaging off. With less
averaging, there is more noise on the data.

If you want to do long term measurements the longest moving average will give
the cleanest results.
"""
function set_moving_average_configuration(device::BrickletBarometerV2, moving_average_length_air_pressure, moving_average_length_temperature)
    device.deviceInternal.set_moving_average_configuration(moving_average_length_air_pressure, moving_average_length_temperature)
end

export get_moving_average_configuration
"""
    $(SIGNATURES)

Returns the moving average configuration as set by
:func:`Set Moving Average Configuration`.
"""
function get_moving_average_configuration(device::BrickletBarometerV2)
    return device.deviceInternal.get_moving_average_configuration()
end

export set_reference_air_pressure
"""
    $(SIGNATURES)

Sets the reference air pressure for the altitude calculation.
Setting the reference to the
current air pressure results in a calculated altitude of 0mm. Passing 0 is
a shortcut for passing the current air pressure as reference.

Well known reference values are the Q codes
`QNH <https://en.wikipedia.org/wiki/QNH>`__ and
`QFE <https://en.wikipedia.org/wiki/Mean_sea_level_pressure#Mean_sea_level_pressure>`__
used in aviation.
"""
function set_reference_air_pressure(device::BrickletBarometerV2, air_pressure)
    device.deviceInternal.set_reference_air_pressure(air_pressure)
end

export get_reference_air_pressure
"""
    $(SIGNATURES)

Returns the reference air pressure as set by :func:`Set Reference Air Pressure`.
"""
function get_reference_air_pressure(device::BrickletBarometerV2)
    return device.deviceInternal.get_reference_air_pressure()
end

export set_calibration
"""
    $(SIGNATURES)

Sets the one point calibration (OPC) values for the air pressure measurement.

Before the Bricklet can be calibrated any previous calibration has to be removed
by setting ``measured air pressure`` and ``actual air pressure`` to 0.

Then the current air pressure has to be measured using the Bricklet
(``measured air pressure``) and with an accurate reference barometer
(``actual air pressure``) at the same time and passed to this function.

After proper calibration the air pressure measurement can achieve an accuracy
up to 0.2 hPa.

The calibration is saved in the EEPROM of the Bricklet and only needs to be
configured once.
"""
function set_calibration(device::BrickletBarometerV2, measured_air_pressure, actual_air_pressure)
    device.deviceInternal.set_calibration(measured_air_pressure, actual_air_pressure)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the air pressure one point calibration values as set by
:func:`Set Calibration`.
"""
function get_calibration(device::BrickletBarometerV2)
    return device.deviceInternal.get_calibration()
end

export set_sensor_configuration
"""
    $(SIGNATURES)

Configures the data rate and air pressure low pass filter. The low pass filter
cut-off frequency (if enabled) can be set to 1/9th or 1/20th of the configure
data rate to decrease the noise on the air pressure data.

The low pass filter configuration only applies to the air pressure measurement.
There is no low pass filter for the temperature measurement.

A higher data rate will result in a less precise temperature because of
self-heating of the sensor. If the accuracy of the temperature reading is
important to you, we would recommend the 1Hz data rate.
"""
function set_sensor_configuration(device::BrickletBarometerV2, data_rate, air_pressure_low_pass_filter)
    device.deviceInternal.set_sensor_configuration(data_rate, air_pressure_low_pass_filter)
end

export get_sensor_configuration
"""
    $(SIGNATURES)

Returns the sensor configuration as set by :func:`Set Sensor Configuration`.
"""
function get_sensor_configuration(device::BrickletBarometerV2)
    return device.deviceInternal.get_sensor_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletBarometerV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletBarometerV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletBarometerV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletBarometerV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletBarometerV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletBarometerV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletBarometerV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletBarometerV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletBarometerV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletBarometerV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletBarometerV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

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

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

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

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


text:
using Gabor

using PyCall
using PyPlot

t = linspace(0,3,1000)
f = sin(5*ceil(t)*2π.*t)
freq,G = gabor(t, f, 0.1)

fig, axs = subplots(nrows=2,sharex="col")
subplot(axs[1])
plot(t,f, "k")
margins(0,0.1)
ylabel(L"$f(t)$")
title("Signal")

subplot(axs[2])
pcolormesh(t,freq,abs(G), rasterized="true",
           cmap = plt[:get_cmap]("viridis"))
ylim(-20,20)

xlabel(L"$t/T$")
ylabel(L"$\omega/2\pi$")
title("Gabor transform")

tight_layout()
savefig("simple.png", dpi=100)


text:
abstract type AbstractPenaltyterm end
const L2NormConstraint = CI{QuadraticDecisionFunction{Float64,QuadraticPart{Float64}}, MOI.LessThan{Float64}}
const LinearizationConstraint = CI{QuadraticDecisionFunction{Float64,LinearPart{Float64}}, MOI.LessThan{Float64}}
const InfNormConstraint = CI{VectorAffineDecisionFunction{Float64}, MOI.NormInfinityCone}
const ManhattanNormConstraint = CI{VectorAffineDecisionFunction{Float64}, MOI.NormOneCone}

Base.copy(::PT) where PT <: AbstractPenaltyterm = PT()

"""
    Quadratic

Functor object for using a quadratic 2-norm penalty term. Requires an `AbstractMathProgSolver` capable of solving QP problems. Passed by default through `penalty` where applicable.

"""
mutable struct Quadratic <: AbstractPenaltyterm
    t::MOI.VariableIndex
    constraint::L2NormConstraint

    Quadratic() = new(MOI.VariableIndex(0), L2NormConstraint(0))
end

function initialize_penaltyterm!(penalty::Quadratic,
                                 model::MOI.AbstractOptimizer,
                                 α::AbstractFloat,
                                 x::Vector{MOI.VariableIndex},
                                 ξ::Vector{MOI.VariableIndex})
    T = typeof(α)
    n = length(x) + 1
    # Get current objective
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    f = MOI.get(model, MOI.ObjectiveFunction{F}())
    # Check if quadratic constraint is supported
    F = MOI.ScalarQuadraticFunction{Float64}
    quad_support = MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})
    if !quad_support
        throw(MOI.UnsupportedAttribute(MOI.ConstraintFunction(), "Using a quadratic penalty term requires an optimizer that supports quadratic constraints"))
    end
    # Add ℓ₂-norm auxilliary variable
    penalty.t = MOI.add_variable(model)
    t = MOI.SingleVariable(penalty.t)
    # Prepare variable vectors
    x = VectorOfDecisions(x)
    ξ = VectorOfKnowns(ξ)
    # Set name
    MOI.set(model, MOI.VariableName(), penalty.t, "‖x - ξ‖₂²")
    # Add quadratic ℓ₂-norm constraint
    g = MOIU.operate(-, T, x, ξ)
    g = LinearAlgebra.dot(g, g)
    MOIU.operate!(-, T, g, t)
    penalty.constraint =
        MOI.add_constraint(model, g,
                           MOI.LessThan(0.0))
    # Add sense-corrected aux term to objective
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    return nothing
end

function update_penaltyterm!(penalty::Quadratic,
                             model::MOI.AbstractOptimizer,
                             α::AbstractFloat,
                             x::Vector{MOI.VariableIndex},
                             ξ::Vector{MOI.VariableIndex})
    # Update penalty parameter
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    # Update projection targets
    MOI.modify(model,
               penalty.constraint,
               KnownValuesChange())
    return nothing
end

function remove_penalty!(penalty::Quadratic,
                         model::MOI.AbstractOptimizer)
    # Delete ℓ₂-norm constraint
    if !iszero(penalty.constraint.value)
        MOI.delete(model, penalty.constraint)
        penalty.constraint = L2NormConstraint(0)
    end
    # Delete aux variable
    if !iszero(penalty.t.value)
        MOI.delete(model, penalty.t)
        penalty.t = MOI.VariableIndex(0)
    end
    return nothing
end

function remove_penalty_variables!(penalty::AbstractPenaltyterm,
                                   list::Vector{MOI.VariableIndex})
    i = something(findfirst(isequal(penalty.t), list), 0)
    if !iszero(i)
        deleteat!(list, i)
    end
    return nothing
end

function remove_penalty_constraints!(penalty::AbstractPenaltyterm,
                                     list)
    # Nothing to do if constraints do not match
    return nothing
end

function remove_penalty_constraints!(penalty::Quadratic,
                                     list::Vector{<:L2NormConstraint})
    i = something(findfirst(isequal(penalty.constraint), list), 0)
    if !iszero(i)
        deleteat!(list, i)
    end
    return nothing
end

"""
    Linearized

Functor object for using an approximately quadratic penalty term, through linearization. Pass through `penalty` where applicable.

...
# Parameters
- `nbreakpoints::Int`: Number of cutting planes used to approximate quadratic term
...
"""
mutable struct Linearized <: AbstractPenaltyterm
    num_breakpoints::Int
    spacing::Float64
    auxilliary_variables::Vector{MOI.VariableIndex}
    constraints::Vector{LinearizationConstraint}

    function Linearized(num_breakpoints::Integer, spacing::Float64)
        n = num_breakpoints >= 3 ? num_breakpoints : 3
        spacing = spacing >= eps() ? spacing : 1.0
        return new(num_breakpoints,
                   spacing,
                   MOI.VariableIndex[],
                   LinearizationConstraint[])
    end
end
Linearized(; num_breakpoints = 3, spacing = 1.0) = Linearized(num_breakpoints, spacing)
Base.copy(linearized::Linearized) = Linearized(linearized.num_breakpoints, linearized.spacing)

function initialize_penaltyterm!(penalty::Linearized,
                                 model::MOI.AbstractOptimizer,
                                 α::AbstractFloat,
                                 x::Vector{MOI.VariableIndex},
                                 ξ::Vector{MOI.VariableIndex})
    T = typeof(α)
    n = length(x)
    m = penalty.num_breakpoints
    resize!(penalty.auxilliary_variables, n)
    resize!(penalty.constraints, n * m)
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    # Add auxilliary cost variables
    for i in eachindex(x)
        penalty.auxilliary_variables[i] = MOI.add_variable(model)
        var = penalty.auxilliary_variables[i]
        MOI.add_constraint(model, MOI.SingleVariable(var), MOI.GreaterThan(0.0))
        name = add_subscript("‖x - ξ‖₂²", i)
        MOI.set(model, MOI.VariableName(), var, name)
    end
    # Add `num_breakpoints` linearization constraints
    # as tangent planes placed uniform around `penalty.spacing`
    # multiplied by the current incumbent `ξ`
    breakpoints = map(1:m) do i
        2*penalty.spacing*(i-1)/(m-1)
    end .- (penalty.spacing - 1)
    k = 1
    for i in eachindex(x)
        tᵢ = MOI.SingleVariable(penalty.auxilliary_variables[i])
        xᵢ = SingleDecision(x[i])
        ξᵢ = SingleKnown(ξ[i])
        for (j,r) in enumerate(breakpoints)
            # Add linearization constraint
            g = MOIU.operate(-, T, xᵢ, r * ξᵢ)
            ∇f = MOIU.operate(*, T, 2*(r - 1), ξᵢ)
            g = MOIU.operate(*, T, ∇f, g)
            f = MOIU.operate(*, T, r - 1, ξᵢ)
            f = MOIU.operate(*, T, f, f)
            MOIU.operate!(+, T, g, f)
            MOIU.operate!(-, T, g, tᵢ)
            penalty.constraints[k] =
                MOI.add_constraint(model,
                                   g,
                                   MOI.LessThan(0.0))
            k += 1
        end
    end
    # Add sense-corrected aux term to objective
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    for i in eachindex(x)
        MOI.modify(model,
                   MOI.ObjectiveFunction{F}(),
                   MOI.ScalarCoefficientChange(penalty.auxilliary_variables[i],
                                               correction * α))
    end
    return nothing
end

function update_penaltyterm!(penalty::Linearized,
                             model::MOI.AbstractOptimizer,
                             α::AbstractFloat,
                             x::Vector{MOI.VariableIndex},
                             ξ::Vector{MOI.VariableIndex})
    # Update penalty parameter
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    for i in eachindex(x)
        MOI.modify(model,
                   MOI.ObjectiveFunction{F}(),
                   MOI.ScalarCoefficientChange(penalty.auxilliary_variables[i], correction * α))
    end
    # Update projection target
    for constraint in penalty.constraints
        MOI.modify(model,
                   constraint,
                   KnownValuesChange())
    end
    return nothing
end

function remove_penalty!(penalty::Linearized,
                         model::MOI.AbstractOptimizer)
    # Delete linearization constraints
    for (i,constraint) in enumerate(penalty.constraints)
        if !iszero(constraint.value)
            MOI.delete(model, constraint)
            penalty.constraints[i] = LinearizationConstraint(0)
        end
    end
    # Delete aux variables
    for (i,var) in enumerate(penalty.auxilliary_variables)
        if !iszero(var.value)
            MOI.delete(model, var)
            penalty.auxilliary_variables[i] = MOI.VariableIndex(0)
        end
    end
    return nothing
end

function remove_penalty_variables!(penalty::Linearized,
                                   list::Vector{MOI.VariableIndex})
    filter!(vi -> !(vi in penalty.auxilliary_variables), list)
    return nothing
end

function remove_penalty_constraints!(penalty::Linearized,
                                     list::Vector{<:LinearizationConstraint})
    filter!(ci -> !(ci in penalty.constraints), list)
    return nothing
end

"""
    InfNorm

Functor object for using a linear ∞-norm penalty term. Pass through `penalty` where applicable.

"""
mutable struct InfNorm <: AbstractPenaltyterm
    t::MOI.VariableIndex
    constraint::InfNormConstraint

    InfNorm() = new(MOI.VariableIndex(0), InfNormConstraint(0))
end

function initialize_penaltyterm!(penalty::InfNorm,
                                 model::MOI.AbstractOptimizer,
                                 α::AbstractFloat,
                                 x::Vector{MOI.VariableIndex},
                                 ξ::Vector{MOI.VariableIndex})
    T = typeof(α)
    n = length(x) + 1
    # Add ∞-norm auxilliary variable
    penalty.t = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), penalty.t, "||x - ξ||_∞")
    x = VectorOfDecisions(x)
    ξ = VectorOfKnowns(ξ)
    t = MOI.SingleVariable(penalty.t)
    # Add ∞-norm constraint
    f = MOIU.operate(vcat, T, t, x) -
        MOIU.operate(vcat, T, zero(α), ξ)
    penalty.constraint =
        MOI.add_constraint(model, f,
                           MOI.NormInfinityCone(n))
    # Add sense-corrected aux term to objective
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    return nothing
end

function update_penaltyterm!(penalty::InfNorm,
                             model::MOI.AbstractOptimizer,
                             α::AbstractFloat,
                             x::Vector{MOI.VariableIndex},
                             ξ::Vector{MOI.VariableIndex})
    # Update penalty parameter
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    # Update projection targets
    MOI.modify(model,
               penalty.constraint,
               KnownValuesChange())
    return nothing
end

function remove_penalty!(penalty::InfNorm,
                         model::MOI.AbstractOptimizer)
    # Delete ∞-norm constraint
    if !iszero(penalty.constraint.value)
        MOI.delete(model, penalty.constraint)
        penalty.constraint = InfNormConstraint(0)
    end
    # Delete aux variable
    if !iszero(penalty.t.value)
        MOI.delete(model, penalty.t)
        penalty.t = MOI.VariableIndex(0)
    end
    return nothing
end

function remove_penalty_constraints!(penalty::InfNorm,
                                     list::Vector{<:InfNormConstraint})
    i = something(findfirst(isequal(penalty.constraint), list), 0)
    if !iszero(i)
        deleteat!(list, i)
    end
    return nothing
end

"""
    ManhattanNorm

Functor object for using a linear 1-norm penalty term. Pass through `penalty` where applicable.

"""
mutable struct ManhattanNorm <: AbstractPenaltyterm
    t::MOI.VariableIndex
    constraint::ManhattanNormConstraint

    ManhattanNorm() = new(MOI.VariableIndex(0), ManhattanNormConstraint(0))
end

function initialize_penaltyterm!(penalty::ManhattanNorm,
                                 model::MOI.AbstractOptimizer,
                                 α::AbstractFloat,
                                 x::Vector{MOI.VariableIndex},
                                 ξ::Vector{MOI.VariableIndex})
    T = typeof(α)
    n = length(x) + 1
    # Add ∞-norm auxilliary variable
    penalty.t = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), penalty.t, "‖x - ξ‖₁")
    x = VectorOfDecisions(x)
    ξ = VectorOfKnowns(ξ)
    t = MOI.SingleVariable(penalty.t)
    # Add ∞-norm constraint
    f = MOIU.operate(vcat, T, t, x) -
        MOIU.operate(vcat, T, zero(α), ξ)
    penalty.constraint =
        MOI.add_constraint(model, f,
                           MOI.NormOneCone(n))
    # Add sense-corrected aux term to objective
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    return nothing
end

function update_penaltyterm!(penalty::ManhattanNorm,
                             model::MOI.AbstractOptimizer,
                             α::AbstractFloat,
                             x::Vector{MOI.VariableIndex},
                             ξ::Vector{MOI.VariableIndex})
    # Update penalty parameter
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    sense = MOI.get(model, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    MOI.modify(model,
               MOI.ObjectiveFunction{F}(),
               MOI.ScalarCoefficientChange(penalty.t, correction * α))
    # Update projection targets
    MOI.modify(model,
               penalty.constraint,
               KnownValuesChange())
    return nothing
end

function remove_penalty!(penalty::ManhattanNorm,
                         model::MOI.AbstractOptimizer)
    # Delete ∞-norm constraint
    if !iszero(penalty.constraint.value)
        MOI.delete(model, penalty.constraint)
        penalty.constraint = ManhattanNormConstraint(0)
    end
    # Delete aux variable
    if !iszero(penalty.t.value)
        MOI.delete(model, penalty.t)
        penalty.t = MOI.VariableIndex(0)
    end
    return nothing
end

function remove_penalty_constraints!(penalty::ManhattanNorm,
                                     list::Vector{<:ManhattanNormConstraint})
    i = something(findfirst(isequal(penalty.constraint), list), 0)
    if !iszero(i)
        deleteat!(list, i)
    end
    return nothing
end


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

using FEMBase
using FEMBase: assemble_mass_matrix!
using Test

struct Dummy <: FieldProblem end
struct Dummy2 <: FieldProblem end

@testset "test tet4 mass matrix" begin
    X = Dict(
             1 => [2.0, 3.0, 4.0],
             2 => [6.0, 3.0, 2.0],
             3 => [2.0, 5.0, 1.0],
             4 => [4.0, 3.0, 6.0])

    element = Element(Tet4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "density", 5.0)

    p = Problem(Dummy, "test", 1)
    add_elements!(p, [element])
    assemble_mass_matrix!(p, 0.0)
    M_expected = [
                  2 1 1 1
                  1 2 1 1
                  1 1 2 1
                  1 1 1 2
                 ]
    @test isapprox(p.assembly.M, M_expected)
    assemble_mass_matrix!(p, 0.0)
    @test isapprox(p.assembly.M, M_expected)
end

@testset "test tet10 mass matrix" begin
    X = Dict(
             1 => [2.0, 3.0, 4.0],
             2 => [6.0, 3.0, 2.0],
             3 => [2.0, 5.0, 1.0],
             4 => [4.0, 3.0, 6.0])
    X[5] = 1/2*(X[1] + X[2])
    X[6] = 1/2*(X[2] + X[3])
    X[7] = 1/2*(X[3] + X[1])
    X[8] = 1/2*(X[1] + X[4])
    X[9] = 1/2*(X[2] + X[4])
    X[10] = 1/2*(X[3] + X[4])

    element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
    update!(element, "geometry", X)
    update!(element, "density", 105.0)

    p = Problem(Dummy, "test", 1)
    add_elements!(p, [element])
    assemble_mass_matrix!(p, 0.0)

    M_expected = [
        6   1   1   1  -4  -6  -4  -4  -6  -6
        1   6   1   1  -4  -4  -6  -6  -4  -6
        1   1   6   1  -6  -4  -4  -6  -6  -4
        1   1   1   6  -6  -6  -6  -4  -4  -4
       -4  -4  -6  -6  32  16  16  16  16   8
       -6  -4  -4  -6  16  32  16   8  16  16
       -4  -6  -4  -6  16  16  32  16   8  16
       -4  -6  -6  -4  16   8  16  32  16  16
       -6  -4  -6  -4  16  16   8  16  32  16
       -6  -6  -4  -4   8  16  16  16  16  32]

    @test isapprox(p.assembly.M, M_expected)

    X[5] .+= 0.1
    update!(element, "geometry", X)
    empty!(p.assembly.M)
    assemble_mass_matrix!(p, 0.0)
end

@testset "add, empty and compare assembly" begin
    A = Assembly()
    B = Assembly()
    @test isempty(A)
    add!(A.K, [1,2,3,4], [1,2,3,4], [1.0 2.0; 3.0 4.0])
    add!(B.K, [1,2,3,4], [1,2,3,4], [1.0 2.0; 3.0 4.0])
    @test isapprox(A, B)
    @test !isempty(A)
    empty!(A)
    @test isempty(A)
end

function FEMBase.assemble_prehook!(problem::Problem{Dummy}, time::Float64)
    add!(problem.assembly.K, [1, 2], [1, 2], [1.0 -1.0; -1.0 1.0])
end

function FEMBase.assemble_posthook!(problem::Problem{Dummy}, time::Float64)
    add!(problem.assembly.K, [2, 3], [2, 3], [1.0 -1.0; -1.0 1.0])
end

@testset "pre and posthooks, warnings" begin
    p = Problem(Dummy, "test", 1)
    p2 = Problem(Dummy2, "test2", 1)
    assemble!(p, 0.0)
    assemble!(p2, 0.0)
    k = [1.0 -1.0; -1.0 1.0]
    K_expected = zeros(3,3)
    K_expected[1:2,1:2] += k
    K_expected[2:3,2:3] += k
    @test isapprox(p.assembly.K, K_expected)
    assemble!(p, 0.0)
    @test isapprox(p.assembly.K, 2.0*K_expected)
end


text:
# Types
# =====

const Coff_t = Clong

# Incomplete Structs
# ------------------

# These types must be used via pointer (aka. opaque pointer).
const BGZF = Cvoid
const hFILE = Cvoid
const hts_tpool = Cvoid
const knetFile = Cvoid
const faidx_t = Cvoid
const hts_opt = Cvoid
const htsThreadPool = Cvoid
const hts_idx_t = Cvoid
const hts_itr_t = Cvoid
const hts_reglist_t = Cvoid
const errmod_t = Cvoid
const probaln_par_t = Cvoid
const hts_md5_context = Cvoid
const tbx_conf_t = Cvoid
const tbx_t = Cvoid
const regidx_t = Cvoid
const regitr_t = Cvoid
const sam_hdr_t = Cvoid
const bam_pileup1_t = Cvoid
const bam_plp_t = Ptr{Cvoid}
const bam_mplp_t = Ptr{Cvoid}
# cram.h
const cram_file_def = Cvoid
const cram_fd = Cvoid
const cram_container = Cvoid
const cram_block = Cvoid
const cram_slice = Cvoid
const cram_metrics = Cvoid
const cram_block_slice_hdr = Cvoid
const cram_block_compression_hdr = Cvoid
const refs_t = Cvoid

# Function Types
# --------------
const kgets_func = Cvoid
const kgets_func2 = Cvoid
const hts_name2id_f = Ptr{Cvoid}
const hts_id2name_f = Ptr{Cvoid}
const hts_itr_query_func = Cvoid
const hts_readrec_func = Cvoid
const hts_seek_func = Cvoid
const hts_tell_func = Cvoid
const hts_itr_multi_query_func = Cvoid
const regidx_parse_f = Ptr{Cvoid}
const regidx_free_f = Ptr{Cvoid}
const bam_plp_auto_f = Ptr{Cvoid}

# Explicit Types
# --------------

const hts_pos_t = Int64

struct hts_pair_pos_t
    beg::hts_pos_t
    end0::hts_pos_t
end

# TODO: Define an initializer.
struct kstring_t
    l::Csize_t
    m::Csize_t
    s::Ptr{Cchar}
end

@enum fai_format_options begin
    FAI_NONE
    FAI_FASTA
    FAI_FASTQ
end

struct htsFormat
    category::Cint  # enum htsFormatCategory
    format::Cint  # enum htsExactFormat
    version::Tuple{Cshort,Cshort}   # major, minor
    compression::Cint  # enum htsCompression
    compression_level::Cshort
    specific::Ptr{Cvoid}
end

struct htsFile
    bitflags::UInt32  # is_bin, is_write, is_be, is_cram, is_bgzf, and dummy
    lineno::Int64
    line::kstring_t
    fn::Ptr{Cchar}
    fn_aux::Ptr{Cchar}
    fp::Ptr{Cvoid}  # union of bgzf, cram, and hfile
    state::Ptr{Cvoid}
    format::htsFormat
    idx::Ptr{hts_idx_t}
    fnidx::Ptr{Cchar}
    bam_header::Ptr{sam_hdr_t}
end

const samFile = htsFile  # typedef alias

struct bam1_core_t
    pos::hts_pos_t
    tid::Int32
    bin::UInt16
    qual::UInt8
    l_extranul::UInt8
    flag::UInt16
    l_qname::UInt16
    n_cigar::UInt32
    l_qseq::Int32
    mtid::Int32
    mpos::hts_pos_t
    isize::hts_pos_t
end

struct bam1_t
    core::bam1_core_t
    id::UInt64
    data::Ptr{UInt8}
    l_data::Cint
    m_data::UInt32
    mempolicy::UInt32
end

@enum cram_block_method begin
    BM_ERROR = -1
    RAW      = 0
    GZIP     = 1
    BZIP2    = 2
    LZMA     = 3
    #RANS     = 4
    RANS0    = 4
    RANS1    = 10
    GZIP_RLE = 11
end
const RANS = RANS0

@enum cram_content_type begin
    CT_ERROR           = -1
    FILE_HEADER        = 0
    COMPRESSION_HEADER = 1
    MAPPED_SLICE       = 2
    UNMAPPED_SLICE     = 3
    EXTERNAL           = 4
    CORE               = 5
end

# hfile_internal.h
# ----------------

# These are not in the public header file, but they are needed for extending hFILE to support arbitrary IO in julia.

struct hFILE_backend
    read::Ptr{Cvoid}
    write::Ptr{Cvoid}
    seek::Ptr{Cvoid}
    flush::Ptr{Cvoid}
    close::Ptr{Cvoid}
end

# This hFILE_base struct represents the header part of hFILE type for the purpose of accessing
# backend field.
# Note that: "you should imagine that hFILE is an opaque incomplete type. They may change in future releases."

struct hFILE_base
    buffer::Ptr{Cchar}
    beginz::Ptr{Cchar}
    endz::Ptr{Cchar}
    limit::Ptr{Cchar}
    backend::Ptr{hFILE_backend}
    offset::Coff_t
    ## remaining fields (as of 06/2021):
    ##      unsigned at_eof:1, mobile:1, readonly:1;
    ##      int has_errno;
    ## just reserve enough space
    remaining::UInt128
end

# functions
@defun hfile_init(struct_size::Csize_t, mode::Ptr{Cchar}, capacity::Csize_t)::Ptr{hFILE}
@defun hfile_destroy(fp::Ptr{hFILE})::Cvoid


text:
# Formulas for representing and working with linear-model-type expressions
# Harlan D. Harris

# we can use Julia's parser and just write them as expressions

type Formula
    lhs::Vector
    rhs::Vector
end

function _dashtree_to_list(ex::Expr)
    # if the head of the expression is a --, return a list of recursive calls to the LHS and RHS
    # otherwise, return what we got as a cell array
    if ex.args[1] == :(--)
        append(_dashtree_to_list(ex.args[2]), _dashtree_to_list(ex.args[3]))
    else
        [ ex ]
    end
end
_dashtree_to_list(ss::Symbol) = [ ss ]

function Formula(ex::Expr)
    # confirm we've got a call to ~, then recursively build a list out of -- - separated subtrees
    if ex.args[1] != :(~)
        error("Formula must have a ~!")
    else
        lhs = _dashtree_to_list(ex.args[2])
        rhs = _dashtree_to_list(ex.args[3])
        Formula(lhs, rhs)
    end
end

function show(io::IO, f::Formula)
    print(io, strcat("Formula: ", join(f.lhs, ", "), " ~ ", join(f.rhs, " +")))
end

# a ModelFrame is just a wrapper around a DataFrame
# construct with mf::ModelFrame = model_frame(f::Formula, d::DataFrame)
# this unpacks the formula, extracts the columns from the DataFrame, and
# applies the ^, &, *, / operators, and evaluates any other expressions
# in the context of the df. Note that columns remain DataVectors, and NAs are
# preserved through this step.

# minimal first version: support y ~ x1 + x2 + log(x3)
type ModelFrame
    df::AbstractDataFrame
    y_indexes::Vector{Int}
    formula::Formula
end

# Obtain Array of Symbols used in an Expr
function unique_symbols(ex::Expr)
    [[unique_symbols(a) for a in ex.args[2:end]]...]
end
unique_symbols(ex::Array{Expr,1}) = [unique_symbols(ex[1])]
unique_symbols(ex::Symbol) = [ex]
unique_symbols(ex::Array{Symbol,1}) = [ex[1]]

# Create a DataFrame containing only variables required by the Formula
# Include grouped indexes present in original DataFrame that are in Formula
# TODO: Split GroupedDataFrame away from this function?
function model_frame(f::Formula, d::AbstractDataFrame)
    cols = {}
    col_names = Array(ASCIIString,0)
    gidx = Dict{Union(UTF8String,ASCIIString),Array{Union(UTF8String,ASCIIString),1}}()

    lhs = unique_symbols(f.lhs)
    rhs = unique_symbols(f.rhs)

    # so, foreach element in the lhs, evaluate it in the context of the df
    for ll = 1:length(lhs)
      s = lhs[ll]
      if is_group(d, string(s))
        for a = get_groups(d)[string(s)]
          push!(cols, with(d, symbol(a)))
          push!(col_names, a)
        end
        gidx[string(s)] = get_groups(d)[string(s)]
      else
        push!(cols, with(d, s))
        push!(col_names, string(s))
      end
    end
    y_indexes = [1:length(lhs)]

    # and foreach unique symbol in the rhs, do the same
    for rr = 1:length(rhs)
      s = rhs[rr]
      if is_group(d, string(s))
        for a = get_groups(d)[string(s)]
          push!(cols, with(d, symbol(a)))
          push!(col_names, a)
        end
        gidx[string(s)] = get_groups(d)[string(s)]
      else
        push!(cols, with(d, s))
        push!(col_names, string(s))
      end
    end
    # bind together and return
    df = DataFrame(cols, col_names)
    if length(gidx) > 0
        set_groups(df, gidx)
    end
    ModelFrame(df, y_indexes, f)
end

model_frame(ex::Expr, d::AbstractDataFrame) = model_frame(Formula(ex), d)

# a ModelMatrix is a wrapper around a matrix, with column names.
# construct with mm::ModelMatrix = model_matrix(mf::ModelFrame, ...)
# this converts any non-numeric types to contrasts, deals with NAs, etc.

# minimal first version: allow numbers and strings only, complete cases only

_parallel_or(a,b) = [(x[1] || x[2])::Bool for x in zip(a, b)]
_parallel_and(a,b) = [(x[1] && x[2])::Bool for x in zip(a, b)]
complete_cases(df::AbstractDataFrame) = reduce(_parallel_and, colwise(x->!isna(x), df))

type ModelMatrix{T}
    model::Array{Float64}
    response::Array{Float64}
    model_colnames::Array{T}
    response_colnames::Array{T}
end

# Expand dummy variables and equations
function model_matrix(mf::ModelFrame)
    ex = mf.formula.rhs[1]
    # BUG: complete_cases doesn't preserve grouped columns
    df = mf.df#[complete_cases(mf.df),1:ncol(mf.df)]
    rdf = df[mf.y_indexes]
    mdf = expand(ex, df)

    # TODO: Convert to Array{Float64} in a cleaner way
    rnames = colnames(rdf)
    mnames = colnames(mdf)
    r = Array(Float64,nrow(rdf),ncol(rdf))
    m = Array(Float64,nrow(mdf),ncol(mdf))
    for i = 1:nrow(rdf)
      for j = 1:ncol(rdf)
        r[i,j] = float(rdf[i,j])
      end
      for j = 1:ncol(mdf)
        m[i,j] = float(mdf[i,j])
      end
    end

    include_intercept = true
    if include_intercept
      m = hcat(ones(nrow(mdf)), m)
      unshift!(mnames, "(Intercept)")
    end

    ModelMatrix(m, r, mnames, rnames)
    ## mnames = {}
    ## rnames = {}
    ## for c in 1:ncol(rdf)
    ##   r = hcat(r, float(rdf[c]))
    ##   push!(rnames, colnames(rdf)[c])
    ## end
    ## for c in 1:ncol(mdf)
    ##   m = hcat(m, mdf[c])
    ##   push!(mnames, colnames(mdf)[c])
    ## end
end

model_matrix(f::Formula, d::AbstractDataFrame) = model_matrix(model_frame(f, d))
model_matrix(ex::Expr, d::AbstractDataFrame) = model_matrix(model_frame(Formula(ex), d))

# TODO: Make a more general version of these functions
# TODO: Be able to extract information about each column name
function interaction_design_matrix(a::AbstractDataFrame, b::AbstractDataFrame)
   cols = {}
   col_names = Array(ASCIIString,0)
   for i in 1:ncol(a)
       for j in 1:ncol(b)
          push!(cols, DataArray(a[:,i] .* b[:,j]))
          push!(col_names, strcat(colnames(a)[i],"&",colnames(b)[j]))
       end
   end
   DataFrame(cols, col_names)
end

function interaction_design_matrix(a::AbstractDataFrame, b::AbstractDataFrame, c::AbstractDataFrame)
   cols = {}
   col_names = Array(ASCIIString,0)
   for i in 1:ncol(a)
       for j in 1:ncol(b)
           for k in 1:ncol(b)
              push!(cols, DataArray(a[:,i] .* b[:,j] .* c[:,k]))
              push!(col_names, strcat(colnames(a)[i],"&",colnames(b)[j],"&",colnames(c)[k]))
           end
       end
   end
   DataFrame(cols, col_names)
end

# Temporary: Manually describe the interactions needed for DataFrame Array.
function all_interactions(dfs::Array{Any,1})
    d = DataFrame()
    if length(dfs) == 2
      combos = ([1,2],)
    elseif length(dfs) == 3
      combos = ([1,2], [1,3], [2,3], [1,2,3])
    else
      error("interactions with more than 3 terms not implemented (yet)")
    end
    for combo in combos
       if length(combo) == 2
         a = interaction_design_matrix(dfs[combo[1]],dfs[combo[2]])
       elseif length(combo) == 3
         a = interaction_design_matrix(dfs[combo[1]],dfs[combo[2]],dfs[combo[3]])
       end
       d = insert!(d, a)
    end
    return d
end

# string(Expr) now quotes, which we don't want. This hacks around that, stealing
# from print_to_string
function formula_string(ex::Expr)
    s = memio(0, false)
    Base.show_unquoted(s, ex)
    takebuf_string(s)
end

#
# The main expression to DataFrame expansion function.
# Returns a DataFrame.
#

function expand(ex::Expr, df::AbstractDataFrame)
    f = eval(ex.args[1])
    if method_exists(f, (FormulaExpander, Vector{Any}, DataFrame))
        # These are specialized expander functions (+, *, &, etc.)
        f(FormulaExpander(), ex.args[2:end], df)
    else
        # Everything else is called recursively:
        expand(with(df, ex), formula_string(ex), df)
    end
end

function expand(s::Symbol, df::AbstractDataFrame)
    expand(with(df, s), string(s), df)
end

# TODO: make this array{symbol}?
function expand(args::Array{Any}, df::AbstractDataFrame)
    [expand(x, df) for x in args]
end

function expand(x, name::ByteString, df::AbstractDataFrame)
    # If this happens to be a column group, then expand each and concatenate
    if is_group(df, name)
      preds = get_groups(df)[name]
      dfs = [expand(symbol(x), df) for x in preds]
      return cbind(dfs...)
    end
    # This is the default for expansion: put it right in to a DataFrame.
    DataFrame({x}, [name])
end

#
# Methods for expansion of specific data types
#

# Expand a PooledDataVector into a matrix of indicators for each dummy variable
# TODO: account for NAs?
function expand(poolcol::PooledDataVector, colname::ByteString, df::AbstractDataFrame)
    newcol = {DataArray([convert(Float64,x)::Float64 for x in (poolcol.refs .== i)]) for i in 2:length(poolcol.pool)}
    newcolname = [strcat(colname, ":", x) for x in poolcol.pool[2:length(poolcol.pool)]]
    DataFrame(newcol, convert(Vector{ByteString}, newcolname))
end

#
# Methods for Formula expansion
#
type FormulaExpander; end # This is an indictor type.

function +(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
    d = DataFrame()
    for a in args
        d = insert!(d, expand(a, df))
    end
    d
end
function (&)(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
    interaction_design_matrix(expand(args[1], df), expand(args[2], df))
end
function *(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
    d = +(FormulaExpander(), args, df)
    d = insert!(d, all_interactions(expand(args, df)))
    d
end


text:
# Autogenerated wrapper script for sgtsnepi_jll for i686-linux-musl
export libsgtsnepi

using FFTW_jll
using cilkrts_jll
JLLWrappers.@generate_wrapper_header("sgtsnepi")
JLLWrappers.@declare_library_product(libsgtsnepi, "libsgtsnepi.so.0")
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, cilkrts_jll)
    JLLWrappers.@init_library_product(
        libsgtsnepi,
        "lib/libsgtsnepi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
############################################################################################
#Import
include("tuning/tuning.jl")
include("algorithm.jl")
include("postprocessing.jl")

############################################################################################
#export


text:
#SignalProcessing: Discrete-time-domain dataset generators
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#Integer representation of polynomial x^p1 + x^p2 + x^p3 + x^p4 + 1
_poly(p1::Int) = one(UInt64)<<(p1-1)
_poly(p1::Int, p2::Int) = _poly(p1) + _poly(p2)
_poly(p1::Int, p2::Int, p3::Int, p4::Int) = _poly(p1) + _poly(p2) + _poly(p3) + _poly(p4)

#==Maximum-length LFSR polynomials/taps (XNOR form)
Ref: Alfke, Efficient Shift Registers, LFSR Counters, and Long Pseudo-Random
     Sequence Generators, Xilinx, XAPP 052, v1.1, 1996.==#
const MAXLFSR_POLYNOMIAL = [
	_poly(64,64) #1: not supported
	_poly(64,64) #2: not supported
	_poly(3,2) #3
	_poly(4,3) #4
	_poly(5,3) #5
	_poly(6,5) #6
	_poly(7,6) #7
	_poly(8,6,5,4) #8
	_poly(9,5) #9
	_poly(10,7) #10
	_poly(11,9) #11
	_poly(12,6,4,1) #12
	_poly(13,4,3,1) #13
	_poly(14,5,3,1) #14
	_poly(15,14) #15
	_poly(16,15,13,4) #16
	_poly(17,14) #17
	_poly(18,11) #18
	_poly(19,6,2,1) #19
	_poly(20,17) #20
	_poly(21,19) #21
	_poly(22,21) #22
	_poly(23,18) #23
	_poly(24,23,22,17) #24
	_poly(25,22) #25
	_poly(26,6,2,1) #26
	_poly(27,5,2,1) #27
	_poly(28,25) #28
	_poly(29,27) #29
	_poly(30,6,4,1) #30
	_poly(31,28) #31
	_poly(32,22,2,1) #32
]

#Creates an ideal step response, using time from ref:
#-------------------------------------------------------------------------------
function Base.step(ref::DataTime; ndel::Index=Index(0), amp=1)
	ndel = value(ndel)
	#ASSERT INCREASING
	t = ref.data.t
	xt = zero(t)
	istart = ceil(Int, 1-(t[1]/step(t))) #succeptible to round-offs?
	istart += ndel
	istart = max(istart, 1)

	for i in istart:length(xt)
		xt[i] = amp
	end
	return DataTime(t, xt)
end

#Creates an ideal pulse response, using time from ref:
#-------------------------------------------------------------------------------
function pulse(ref::DataTime; ndel::Index=Index(0), amp=1, npw::Index=Index(0))
	ensure(value(npw) > 0, ArgumentError("npw must be positive"))
	t = ref.data.t
	xt = (step(ref, ndel=ndel, amp=amp).data.xt .-
		step(ref, ndel=ndel+npw, amp=amp).data.xt)
	return DataTime(t, xt)
end

#Creates pulse response of a single-pole system, using time from ref.x:
# tpw: pulse width
#-------------------------------------------------------------------------------
function pulse(ref::DataTime, p::Pole; ndel::Index=Index(0), amp=1, npw::Index=Index(0))
	ndel = value(ndel); npw = value(npw)
	ensure(npw > 0, ArgumentError("npw must be positive"))
	t = ref.data.t
	#Simple wrapper using continuous-time algoritm, for now:
	tdel = step(t)*ndel
	tpw = step(t)*npw
	Π = pulse(DataF1(t), p, tdel=tdel, amp=amp, tpw=tpw)
	return DataTime(ref.data, Π.y)
end

#Creates PRBS sequence, using time from ref.x:
#TODO: should the vector by of type Bool?
#-------------------------------------------------------------------------------
function prbs(; reglen::Int=32, seed::Integer=11, nsamples::Int=0)
	ensure(nsamples > 0, ArgumentError("nsamples must be positive"))
	ensure(in(reglen, 3:32), ArgumentError("Invalid value: 3 <= reglen <= 32"))

	lfsr = UInt64(seed)
	msb = UInt64(1)<<(reglen-1)
	seq = zeros(Bool, nsamples)
	poly = UInt64(MAXLFSR_POLYNOMIAL[reglen])
	mask = ~poly
	#==Since 1 XNOR A => A, we can ignore all taps that are not part of the
	polynomial, simply by forcing all non-tap bits to 1 (OR-ing with ~poly)
	==#

	for i in 1:length(seq)
		reg = lfsr | mask
		bit = msb
		for j in 1:reglen
			bit = ~xor(reg, bit)
			reg <<= 1
		end
		bit = Int((bit & msb) > 0) #Convert resultant MSB to an integer
		seq[i] = bit
		lfsr = (lfsr << 1) | bit
	end
	return seq
end

#Creates a signal pattern from a sequence of data points, and a pulse response:
# seq: a binary data sequence
# p:Pulse response
# nbit: Number of sampling steps in a bit period.
#TODO: shift x instead of y???
#TODO: Add DataTime vectors directly once implemented
function pattern(seq::Vector{T}, p::DataTime; nbit::Index=Index(0)) where {T<:Number}
	nbit = value(nbit)
	ensure(nbit > 0, ArgumentError("nbit must be positive"))
	result = zero(p.data.xt)
	for i in 1:length(seq)
		result += seq[i]*yshift(p, Index((i-1)*nbit)).data.xt
	end
	return DataTime(p.data, result)
end

#Last line


text:
function fmwrite(a, filename::ASCIIString)

    stream = open(filename,"w");
    dummy = Any[]; finalizer(dummy, (x)->close(stream))

    j2m, m2j = typedictionaries()
    fmwrite_internal(stream, a, j2m)
    close(stream)

end

function fmread(filename::ASCIIString)

    stream = open(filename,"r");
    dummy = Any[]; finalizer(dummy, (x)->close(stream))

    j2m, m2j = typedictionaries()
    r = fmread_internal(stream, m2j)
    close(stream)
    return r
end

function typedictionaries()

    j2m = Dict(Int => "int64", Int8 => "int8", Int16 => "int16", Int32 => "int32", Int64 => "int64",
        UInt => "uint64", UInt8 => "uint8", UInt16 => "uint16", UInt32 => "uint32", UInt64 => "uint64",
        Float32 => "single", Float64 => "double", ASCIIString => "char*1", Bool => "logical", Any => "cell",
        Dict => "struct")

    m2j = Dict{ASCIIString,Type}()

    for (k,v) in j2m
        newv = v * repeat("-", 7-length(v))
        j2m[k] = newv
        m2j[newv] = k
    end

    return j2m, m2j

end

function fmwrite_internal(stream, a, j2m)

    ###################
    #  write type

    if typeof(a)==Char
        a = string(a)
    end
    if typeof(a)<:Array
        t = eltype(a)
    else
        t = typeof(a)
    end
    if typeof(a)<:Dict
        t = Dict
    end

    write(stream, j2m[t]);

    ###################
    #  write size

    s = zeros(UInt64, 1, 10)

    if typeof(a)==ASCIIString
        s[1:2] = [1 length(a)]
    elseif typeof(a)<:Dict || ndims(a)==0
        s[1:2] = [1 1]
    else
        if ndims(a)==1
            s[1:2] = [1, size(a,1)];
        else
            for i = 1:ndims(a)
                s[i] = size(a,i)
            end
        end
    end
    if s[1]==0
        s[1] = 1
    end
    if s[2]==0
        s[2] = 1
    end

    #println("s: ",s)

    assert(sum(s)>0);
    write(stream, s);

    ###################
    #  write data

    if typeof(a)<:Dict
        list = Any[]
        for (k,v) in a
            push!(list, [k,v])
        end
        fmwrite_internal(stream, list, j2m)

    elseif t<:Number #|| (typeof(a)<:Array && eltype(a)<:Number)
        write(stream, a)

    elseif t==Any
        for i = 1:length(a)
            fmwrite_internal(stream, a[i], j2m)
        end

    else
        write(stream, a)
    #else
    #    error("how did we get here? typeof(a)=",typeof(a))
    end
end

function fmread_internal(stream, m2j)

    ###################
    #  read type

    t = m2j[ascii(read(stream, UInt8, 7))];

    ###################
    #  read size

    s = read(stream, UInt64, 1, 10)
    s = s[s.!=0]
    println("\n---\ngot here: s ",s," sum(s): ",prod(s),"\n")
    println("type is: $t")
    @show nb_available(stream)

    ###################
    #  read data

    if t<:Dict
        list = fmread_internal(stream, m2j)
        r = Dict()
        for i = 1:length(list)
            r[list[i][1]] = list[i][2]
        end

    elseif t==Any
        r = Array(t, s...)
        for i = 1:length(r)
            r[i] = fmread_internal(stream, m2j)
        end

    elseif t==ASCIIString
        r = ascii(read(stream, UInt8, prod(s)))
    else
        r = Array(t, s...)
        read!(stream, r)
        read(stream, r)
    end

    return r
end


text:
"""
function print_ascii(root::N) where N<:AbstractNode

Prints the ASCII representation of the tree.
Adapted for Julia from the ETE toolkit implementation.

* `root`: Root node of tree to print.
"""
function print_ascii(root::N) where N<:AbstractNode
lines,_ = ascii(root,"-")
println("\n"*join(lines,"\n"))
end #function

function ascii(root::N,char1="-") where N<:AbstractNode
    node_name = root.name
    LEN = max(3,length(node_name))
    PAD = " "^LEN
    PA = " "^(LEN-1)
    if root.nchild != 0
        mids = []
        result = []
        for c in root.children
            if length(root.children) == 1
                char2 = "-"
            elseif c.num == root.children[1].num
                char2 = "/"
            elseif c.num == root.children[end].num
                char2 = "\\"
            else
                char2 = "-"
            end #ifelse
            clines,mid = ascii(c,char2)

            push!(mids,mid+length(result))
            for x in clines
                push!(result,x)
            end #for
            push!(result,"")
        end #for
        pop!(result)
        (lo, hi, final) = (mids[1],mids[end],length(result))
        prefixes = []
        if lo+1 > 0
            for x in 1:(lo+1)
                push!(prefixes,PAD)
            end #for
        end #if
        if hi-lo-1 > 0
            for  x in 1:(hi-lo-1)
                push!(prefixes,PA*"|")
            end #for
        end #if
        if final-hi >0
            for x in 1:(final-hi)
                push!(prefixes,PAD)
            end #for
        end #if
        mid = (Integer(round((lo+hi)/2,RoundDown)))+1

        prefixes[mid] = char1 * ("-"^(LEN-2)) * (prefixes[mid][end])
        if length(prefixes) == length(result)
            result = [p*l for (p,l) in zip(prefixes,result)]
        else
            nuresult = []
            for x in 1:min(length(prefixes),length(result))
                push!(nuresult,prefixes[x]*result[x])
            end #for
            result = nuresult
        end #ifelse
        stem = result[mid]
        #added to handle spacing between a parent and child, each with a long name
        optional = ""
        if !(length(stem) < length(node_name)+2) && isletter(stem[(length(node_name)+2)])
            optional = "-"
        end #if

        result[mid] = stem[1] * node_name * optional * stem[(length(node_name)+2):end]

        return result,mid-1
    else
        return [char1 * "-" * node_name], 0
    end #if
end #function


text:
using ReinforcementLearningEnvironmentViZDoom, Flux, GR, JLD2
using Compat: @info
const withgpu = true
const gpudevice = 0

if withgpu
    using CUDAnative
    device!(gpudevice)
    using CuArrays
    const inputdtype = Float32
else
    const inputdtype = Float64
end
env = ViZDoomEnvironment(living_reward = .1, available_buttons =
[:ATTACK, :MOVE_LEFT, :MOVE_RIGHT, :MOVE_FORWARD, :TURN_LEFT, :TURN_RIGHT])
model = Chain(x -> x./inputdtype(255), Conv((8, 8), 2 => 16, relu, stride = (4, 4)),
                         Conv((4, 4), 16 => 32, relu, stride = (2, 2)),
                         x -> reshape(x, :, size(x, 4)),
                         Dense(2592, 256, relu),
                         Dense(256, length(env.actions)));
learner = DQN(model, opttype = x -> Flux.ADAM(x, .0001),
              loss = huberloss, doubledqn = true,
              updatetargetevery = 2500, nsteps = 10,
              updateevery = 4, replaysize = 10^6, nmarkov = 2,
              startlearningat = 200000);
pchain = [ImageResizeNearestNeighbour((84, 84)),
          x -> reshape(x, (84, 84, 1))]
if withgpu push!(pchain, togpu) end
preprocessor = ImagePreprocessor((160, 120), pchain)
x = RLSetup(learner,
            env,
            ConstantNumberSteps(4 * 10^6),
            preprocessor = preprocessor,
            callbacks = [Progress(5*10^2), EvaluationPerEpisode(TotalReward()),
                         LinearDecreaseEpsilon(5 * 10^4, 10^6, 1, .01)]);
@info "start learning."
@time learn!(x)
@save "modeldeadlycorridor.jld2" model

max_stars_repo_name: jbrea/RLEnvVizDoom

max_issues_repo_name: jbrea/RLEnvVizDoom

max_forks_repo_name: jbrea/RLEnvVizDoom


text:
using Test, OrdinaryDiffEq, HyperbolicDiffEq, Statistics

function calc_error(balance_law, uₐₙₐ, tmin, tmax, numflux, N, usethreads)
    u₀ = x -> uₐₙₐ(tspan[1], x)

    tspan = (tmin, tmax)
    mesh = UniformPeriodicMesh1D(-2., 2., N)
    fv = FirstOrderFV(balance_law, mesh, numflux, usethreads)
    ode = semidiscretise(fv, u₀, tspan)
    sol = solve(ode, OrdinaryDiffEq.Euler(), dt=max_dt(tspan[1], ode.u0, fv),
                save_everystep=false, dense=false)
    uana = compute_coefficients(x->uₐₙₐ(tspan[end],x), mesh)

    N, integrate(u->sum(abs.(u)), sol[end]-uana, mesh)
end

function calc_order_estimate(Ns_and_errors)
    orders = zeros(length(Ns_and_errors)-1)
    for i in 1:length(orders)
        orders[i] = -log(Ns_and_errors[i+1][2]/Ns_and_errors[i][2]) /
                     log(Ns_and_errors[i+1][1]/Ns_and_errors[i][1])
    end
    mean(orders)
end

function calc_order_estimate(balance_law, uₐₙₐ, tspan, numflux, Ns, usethreads=false)
    calc_order_estimate(calc_error.(balance_law, uₐₙₐ, tspan..., numflux, Ns, usethreads))
end

# test convergence orders
Ns = 100 .* 2 .^ (1:6)

balance_law = ConstantLinearAdvection()
tspan = (0., 0.5)
uₐₙₐ = solve(RiemannProblem(balance_law, 0., 1., -0.5) *
                RiemannProblem(balance_law, 1., 0., 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, GodunovFlux(), Ns) > 0.45
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.45
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.45

balance_law = Burgers()
tspan = (0., 0.5)
uₐₙₐ = solve(RiemannProblem(balance_law, 0., 1., -0.5) *
                RiemannProblem(balance_law, 1., 0., 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, GodunovFlux(), Ns) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.8

balance_law = BuckleyLeverette()
tspan = (0., 0.5)
uₐₙₐ = solve(RiemannProblem(balance_law, 0., 1., -0.5) *
                RiemannProblem(balance_law, 1., 0., 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, GodunovFlux(), Ns) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.8

balance_law = ShallowWater()
tspan = (0., 0.1)
u1 = variables(balance_law)(1., 2.)
u2 = variables(balance_law)(2., 0.)
uₐₙₐ = solve(RiemannProblem(balance_law, u1, u2, -1.) *
                RiemannProblem(balance_law, u2, u1, 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        LocalLaxFriedrichsDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        ScalarDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        MatrixDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux2Param(-1., 1.),
                                                        LocalLaxFriedrichsDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, SuliciuFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, KineticFlux(), Ns) > 0.7

balance_law = ShallowWater(10.)
tspan = (0., 0.1)
u1 = variables(balance_law)(1., 2.)
u2 = variables(balance_law)(2., 0.)
uₐₙₐ = solve(RiemannProblem(balance_law, u1, u2, -1.) *
                RiemannProblem(balance_law, u2, u1, 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        LocalLaxFriedrichsDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        ScalarDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux1Param(-1.),
                                                        MatrixDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, FluxPlusDissipation(
                                                        EnergyConservativeFlux2Param(-1., 1.),
                                                        LocalLaxFriedrichsDissipation()
                                                    ), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, SuliciuFlux(), Ns) > 0.7
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, KineticFlux(), Ns) > 0.7

balance_law = HyperbolicDiffEq.Euler()
tspan = (0., 0.25)
u1 = conserved_variables(1.0, 0.0, 1.0, balance_law)
u2 = conserved_variables(0.125, 0.0, 0.1, balance_law)
uₐₙₐ = solve(RiemannProblem(balance_law, u1, u2, -1.) *
                RiemannProblem(balance_law, u2, u1, 0.5))
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, LocalLaxFriedrichsFlux(), Ns) > 0.55
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, HartenLaxVanLeerFlux(), Ns) > 0.55
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, SuliciuFlux(), Ns) > 0.55


text:
include((@__DIR__)*"/../preamble.jl")

include("absorbing_model/model_def.jl")
include("default_params.jl")

tvec = CSV.read("hitting_times_model/hitting_times/data/point_mass/order_1_model_dg1.csv",DataFrame).t

pth = mkpath("hitting_times_model/hitting_times/data/errors")
for ic_string in ["point_mass","exp"]
    ks_errors = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    l1_errors = DataFrame(DG=[],DG_limit=[],Order_1=[],QBDRAP=[])
    ks_error_row = zeros(4)
    l1_error_row = zeros(4)

    sim_cdf = CSV.read("hitting_times_model/hitting_times/data/"*ic_string*"/sims_cdf_evaluated.csv",DataFrame)
    for order in 1:2:21
        for (c,model_string) in enumerate(["dg1","dg2","order13","qbdrap4"])
            approx_data = CSV.read("hitting_times_model/hitting_times/data/"*
                ic_string*"/order_"*string(order)*"_model_"*model_string*".csv",DataFrame)
            ks_error = max(
                maximum(abs.(approx_data.phase_1-sim_cdf.phase_1)),
                maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
            )
            l1_cdf = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                    sim_cdf.phase_1[1:end-1]-sim_cdf.phase_1[2:end]).*diff(approx_data.t))/2 +
                sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                    sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.t))/2

            ks_error_row[c] = ks_error
            l1_error_row[c] = l1_cdf
        end
        push!(ks_errors,ks_error_row)
        push!(l1_errors,l1_error_row)
    end
    CSV.write(pth*"/ks_"*ic_string*".csv",ks_errors)
    CSV.write(pth*"/l1_"*ic_string*".csv",l1_errors)
end


text:
struct SLInvalidPath <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::SLInvalidPath)
    print(io, ex.msg)
end


text:
using Multilane
using StatsBase
using MCTS
using JLD

solvers = Dict{String, Any}()

dpws = DPWSolver(depth=20,
                 n_iterations=100,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

solvers["dpw_100"] = deepcopy(dpws)
dpws.n_iterations = 300
solvers["dpw_300"] = deepcopy(dpws)
dpws.n_iterations = 500
solvers["dpw_500"] = deepcopy(dpws)
dpws.n_iterations = 1000
solvers["dpw_1000"] = deepcopy(dpws)

curve = TestSet(lambda=[1., 10., 17.78, 31.62, 56.23, 100.], N=500)

tests = []
for (s,_) in solvers
    push!(tests, TestSet(curve, solver_key=s))
end

objects = gen_initials(tests)

@show objects["param_table"]
objects["solvers"] = solvers

results = evaluate(tests, objects)

filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results["histories"] = nothing
save(filename, results)
println("results saved to $filename")


text:
#
# Demo program for timestepping unsteady flows
#
# Currently there's a memory leak in HDF5.jl, so
# this runner will shut down after a certain number
# of steps. To restart automatically, run it as
# while true; do julia output_scripts/demo_unsteady.jl; done

using Revise
using Compat.@info
using JLD
using FileIO
using HDF5

include("timestep_base.jl")

# Bounding box
boxW, boxH = 5, 3

# flow_case = :obstacles
# Re = 30

flow_case = :vortex
Re = 25
offset = true

if flow_case == :obstacles
    basedir = "data/obstacles"
elseif flow_case == :vortex
    if offset
        basedir = "data/vortex_Re$Re"
    else
        basedir = "data/vortex_Re$(Re)_center"
    end
else
    error("Unknown: $(flow_case)")
end
mkpath(basedir)

# Cache for multiple runs
if !isdefined(:solvers)
    solvers = Dict()
end

# Resume timestepping
geoname = "$basedir/geo.jld"
dataname = i -> @sprintf("%s/vstep_%06d.jld", basedir, i)
resume = 0
while isfile(dataname(resume+1))
    resume += 1
end
@show resume

function hook_quiver(solver, sol, ti, i)
    # Unpack
    U1 = sol.U1
    U2 = sol.U2
    X = solver.X
    Y = solver.Y
    dcurve = solver.dcurve
    interior = solver.interior

    figure(1)
    clf()
    k = 10 # Point skip in quiver plot
    quiver(X[1:k:end, 1:k:end], Y[1:k:end, 1:k:end], U1[1:k:end, 1:k:end], U2[1:k:end, 1:k:end])
    axis("equal")
end

function hook_save(solver, sol, ti, i)
    @time begin
        if i==1
            filename = geoname
            println("Writing ", filename)
            save(filename,
                 "X", solver.X, "Y", solver.Y,
                 "interior", solver.interior,
                 "dcurve", solver.dcurve
                 )
        end
        filename = dataname(i)
        println("Writing ", filename)
        U1    = sol.U1
        U2    = sol.U2
        dU1dx = sol.dU1dx
        dU1dy = sol.dU1dy
        dU2dx = sol.dU2dx
        dU2dy = sol.dU2dy
        save(File(format"JLD", filename),
             "ti", ti,
             "U1"    , U1   ,
             "U2"    , U2   ,
             "dU1dx" , dU1dx,
             "dU1dy" , dU1dy,
             "dU2dx" , dU2dx,
             "dU2dy" , dU2dy,
             "RHS1"  , sol.RHS1,
             "RHS2"  , sol.RHS2
             )
    end
end

function hook_vortex(solver, sol, ti, i)
    hook_save(solver, sol, ti, i)
    if mod(i, 1)==0
        println("Vorticity plot:")
        @time hook_vorticity(solver, sol, ti, i)
    end
    if mod(i, 20)==0
        hook_streamlines(solver, sol, ti, i)
    end
end

function hook_vorticity(solver, sol, ti, i)
    figure(1)
    clf()
    X = solver.X
    Y = solver.Y
    vorticity = sol.dU2dx - sol.dU1dy
    s = 1
    pcolormesh(X[1:s:end,1:s:end], Y[1:s:end,1:s:end], vorticity[1:s:end,1:s:end],
           cmap=PyPlot.cm[:coolwarm])
    clim(-3,3)
    plot_boundary(solver.dcurve)
    axis("image")
    axis("off")
end

function hook_streamlines(solver, sol, ti, i)
    figure(2)
    clf()
    U1 = sol.U1
    U2 = sol.U2
    X = solver.X
    Y = solver.Y
    dcurve = solver.dcurve
    interior = solver.interior
    ystart = collect(linspace(-boxH+0.05,boxH-0.05,40))
    xstart = -boxW+0.1 + zeros(ystart)
    start_points = [xstart ystart]
    streamplot(y=Y[1,:],
               x=X[:,1],
               u=U1',
               v=U2',
               density=20,
               start_points=start_points,
               integration_direction="both",
               maxlength=6)

    plot_boundary(dcurve)
    axis("image")
    axis("off")
end

function demo_unsteady(flow_case, Re, resume)
    # Defaults
    fastdirect_TOL = 1e-8

    # Geometry setup
    if flow_case == :obstacles
        numpanels1 = 500
        numpanels2 = 50
        Ngrid = 500
        dt = 0.01
        sdc_order = 2
        pux_R = 0.4
        Nstep = 10000

        curve1 = AnalyticDomains.rectangle(width=boxW, height=boxH)
        # Setup a bunch of starfish
        srand(1)
        curve2 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = (-3.0, -1.0), rotation=2*pi*rand(), interior=true)
        curve3 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = (-1.0, -0.1), rotation=2*pi*rand(), interior=true)
        curve4 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = (-0.9,  1.5), rotation=2*pi*rand(), interior=true)
        curve5 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = ( 1.0,  0.0), rotation=2*pi*rand(), interior=true)
        curve6 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = (-3.5,  0.7), rotation=2*pi*rand(), interior=true)
        curve7 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = ( 1.7, -2.0), rotation=2*pi*rand(), interior=true)
        curve8 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = (-0.7, -1.9), rotation=2*pi*rand(), interior=true)
        curve9 = AnalyticDomains.starfish(;amplitude = 0.2, radius = 0.5, center = ( 2.9,  1.7), rotation=2*pi*rand(), interior=true)

        curve = [curve1, curve2, curve3, curve4, curve5, curve6, curve7, curve8, curve9]
        numpanels = [numpanels1, numpanels2, numpanels2, numpanels2, numpanels2, numpanels2, numpanels2, numpanels2, numpanels2]
    elseif flow_case == :vortex
        if Re==200
            # Stable, but slow
            # Oscillations around 2500
            Re = 200
            dt = 0.01/2
            Ngrid = 500*2
            pux_R = 0.5/2
        elseif Re==100
            # Separates after ~2000 steps with offset
            Re = 100
            dt = 0.01
            Ngrid = 500
            pux_R = 0.5
        elseif Re==50
            # offset and SDC=1: 3500 steps oscillations visible
            Re = 50
            dt = 0.02
            Ngrid = 500
            pux_R = 0.5
        elseif Re==25
            # No oscillations after 10000 steps
            Re = 25
            dt = 1/Re
            Ngrid = 500
            pux_R = 0.5
        end

        sdc_order = 2
        Nstep = 20000

        numpanels1 = 500
        numpanels2 = 100
        if offset
            cylinder = (-boxW+1.5, -0.1) # Offset triggers faster separation
        else
            cylinder = (-boxW+1.5, 0.0) # No offset
        end

        curve1 = AnalyticDomains.rectangle(width=boxW, height=boxH)
        curve2 = AnalyticDomains.starfish(;amplitude = 0.0, radius = 0.5, center = cylinder, interior=true)
        curve = [curve1, curve2]
        numpanels = [numpanels1, numpanels2]
    else
        error()
    end

    # Restart every NN steps (HDF5.jl memory leak!)
    restart = Ngrid > 500 ? 300 : 1000
    Nstep = min(Nstep, resume+restart)
    @info("Will run steps $resume:$Nstep on flow case $flow_case")

    # Solver precomputation
    dt_substep = sdc_substep(dt, sdc_order)
    alpha = sqrt(Re/dt_substep)
    @info("alpha = $alpha")
    if haskey(solvers, alpha)
        println(" * USING CACHE")
        solver = solvers[alpha]
    else
        solver = precompute_solver(alpha=alpha,
                                   pux_R=pux_R,
                                   Ngrid=Ngrid,
                                   curve = curve,
                                   numpanels = numpanels,
                                   fastdirect_TOL = fastdirect_TOL=1e-8
                                   )
        solvers[alpha] = solver
    end

    # Set BCs
    dcurve = solver.dcurve
    ubc1 = zeros(dcurve.numpoints)
    ubc2 = zeros(dcurve.numpoints)
    if flow_case == :obstacles
        # Flow on outer boundary
        @. ubc1[dcurve.curvenum==1] = 1.0
    elseif flow_case == :vortex
        # # Set flow rate on left/right box boundaries
        # idx1 = 1:dcurve.panelorder*numpanels1
        # xbdry = dcurve.points[1,idx1]
        # ybdry = dcurve.points[2,idx1]
        # sigmax = boxW/10
        # u1bdry(x,y) = (exp(-((x-boxW)/sigmax).^2) + exp(-((x+boxW)/sigmax).^2))
        # ubc1[idx1] = u1bdry.(xbdry,ybdry)
        # @. ubc1[idx1] = 1.0

        # Flow on outer boundary
        @. ubc1[dcurve.curvenum==1] = 1.0
    else
        error("Don't know what BCs to set")
    end

    # Estimate BC resolution
    panelorder = dcurve.panelorder
    L = mss.legendre_matrix(panelorder)
    bc_resolution = 0.0
    for i=1:dcurve.numpanels
        xcoeff = L*ubc1[(1:panelorder) + panelorder*(i-1)]
        ycoeff = L*ubc2[(1:panelorder) + panelorder*(i-1)]
        xres = max(abs(xcoeff[end]), abs(xcoeff[end-1]))
        yres = max(abs(ycoeff[end]), abs(ycoeff[end-1]))
        res = max(xres, yres)
        bc_resolution = max(bc_resolution, res)
    end
    info("BC resolution: ", bc_resolution)

    # Start from rest
    X = solver.X
    U1, U2 = zeros(X), zeros(X)
    RHS1, RHS2 = zeros(X), zeros(X)
    dU1dx, dU1dy = zeros(X), zeros(X)
    dU2dx, dU2dy = zeros(X), zeros(X)
    sol0 = Solution(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy, RHS1, RHS2, [], [])

    t0 = 0.0
    istart = 0

    hook = hook_vortex

    if resume > 0
        info("Resuming from ", resume)
        istart = resume
        filename = dataname(istart)
        h5open(filename, "r") do io
            U1 = read(io, "U1")
            U2 = read(io, "U2")
            dU1dx = read(io, "dU1dx")
            dU1dy = read(io, "dU1dy")
            dU2dx = read(io, "dU2dx")
            dU2dy = read(io, "dU2dy")
            RHS1 = read(io, "RHS1")
            RHS2 = read(io, "RHS2")
            sol0 = Solution(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy, RHS1, RHS2, [], [])
        end
    end

    t, sol = timestep(Nstep, dt, sdc_order, hook, solver, t0, sol0, ubc1, ubc2; istart=istart)

end

demo_unsteady(flow_case, Re, resume)


text:
# Collection of test programs to help with testing and debugging of the Simulator.
module TestPrograms
    # Wouldn't make much sense to have a collection of test programs if we
    # didn't include the simulator module now would it?
    using AsapSim

    export test_core, testfifo_core, runfor, writeto!, readfrom!

    #  --- Test reading from all sources  --- #
    include("Sources.jl")
    # --- Test writing to all destinations --- #
    include("Destinations.jl")
    # --- Test branch behavior --- #
    include("Branch.jl")
    # --- Test RPT Instruction --- #
    include("Repeat.jl")
    # --- Test Conditional Execution --- #
    #include("CondExec.jl") (TODO)

    # Simple programs
    include("Simple.jl")
    include("Sort.jl")
    include("Temp.jl")

    ####################
    # Helper Functions #
    ####################
    function runfor(core, cycles)
        for _ in 1:cycles
            AsapSim.update!(core)
        end
    end

    function writeto!(fifo, payload, index)
        # Need to clock fifo to get it to update.
        for i in 1:AsapSim.buffersize(fifo)
            writeupdate!(fifo)
        end
        while AsapSim.iswriteready(fifo) && index <= length(payload)
            write(fifo, payload[index])
            writeupdate!(fifo)
            index += 1
        end
        return index
    end

    function readfrom!(fifo, buffer)
        for i in 1:AsapSim.buffersize(fifo)
            readupdate!(fifo)
        end
        while AsapSim.isreadready(fifo)
            push!(buffer, read(fifo))
            # Read from the fifo and mark it for incrementing.
            AsapSim.increment!(fifo)
            readupdate!(fifo)
        end
    end

    # Create a test-core with some fields initialized for easier testing.
    const RETURN_ADDRESS = 10
    # Set the accumulator to all 1's. Make sure we only grab the lower
    # 16 bits.
    const ACCUMULATOR    = 2^40 - 1

    # Initial DMEM contents
    # Address (index 1) => value
    const DMEM = Dict(
        128 + 1 => 128,
        129 + 1 => 129,
        130 + 1 => 130,
        131 + 1 => 131,

        # Initialize some values for the address generator test.
        12 + 1 => 10,
        14 + 1 => 20,
        16 + 1 => 30
    )

    function test_core(program)
        core = AsapCore(;
            # Initialize return address to test reading from it.
            return_address = RETURN_ADDRESS,
            accumulator = ACCUMULATOR,

            # Initialize pointers to DMEM addresses. Will initialize DMEM
            # addresses to verify these were dereferenced correctly.
            pointers = [128, 129, 130, 131],
        )

        # Initialize DMEM.
        for (k,v) in DMEM
            core.dmem[k] = v
        end

        # Set address generator 1 (0)
        core.address_generators[1] = AsapSim.AddressGenerator(
            start   = 10,
            stop    = 30,
            current = 12,
            stride  = 2,
        )

        # Attach the function to the core.
        core.program = program

        return core
    end

    function testfifo_core(program)
        core = AsapCore(;
            program = program,
            fifos = [AsapSim.TestFifo{Int16}(),AsapSim.TestFifo{Int16}()],
            outputs = AsapSim.SmallAssociative{Int,AsapSim.TestFifo{Int16}}(),
            #outputs = Dict{Int,AsapSim.TestFifo{Int16}}(),
        )

        return core
    end
end


text:
module Milk

using Reexport

@reexport using XGrad
@reexport using NNlib
@reexport using GradDescent

export
    # utils
    xavier_init,
    class_to_index,
    one_hot,
    # gradient descent
    ModelOptimizer,
    update_params!,
    # losses
    cross_entropy_loss,
    # other
    accuracy,
    to_cuda

include("core.jl")

end # module