text:
#### for module multiobjectiveGraphs ####
export genericWeightCategory, weightCategory, objTypes # abstract type and list of concrete
# 6 predefined types
export weightMinSum, weightMinProduct, weightMinMax, weightMaxSum, weightMaxProduct, weightMaxMin
# Programming Interface
export typeCode, defaultValue, mimic, compare, combine , descr

##############################################################
#################### generic type for one objective value ####
##############################################################
abstract type genericWeightCategory end
abstract type weightCategory{Tobjval} <: genericWeightCategory end

###################  weight types supervisor #####################
struct weightCategorySupervisor

end

############################################################
# each concrete weightCategory must implement
# 	... typeCode ... wich returns a small positive integer
# 		specific to this weightCategory and so is unique
# 			see function objTypes for existing values
# ... descr ... wich return the description of the objective type and value
# 	... defaultValue ... wich returns the base value for combining.  ###### really needed  ?? -> change initializations !!! , a arbitary base value may suffice
# 		e.g. 0 for sum, 1 for product
# 	... mimic ... wich yields a new element similar to its first parameter with or without value (several functions)
# 	... compare(a,b) ... whose result is a couple of booleans (a is better than b or equals b, b is better than a or equals a)
#      compare must be a total order (especially  : forall (a,b) if a<=b & b<=a then a=b)
# 	... combine ... to create a new weight from two such weights
# 		that can be a sum, a product, a min, a max or anything
# 		that is associative ( combine(x,combine(y,z)) == combine(combine(x,y),z) )
#      and monotonous ( i.e. compare(x, combine(x,y)) must always yield (true, _) ) ### this constraint can be relaxed with appropriate detection in algorithms
#      in order to ensure there is no "negative weight" circuit in the graph
# don't forget to add each new type to the array (by redefining the function objTypes ???). ### must be tested
############################################################

#############################################################
################# objective types coding ##########################
function objTypes( a::weightCategory{Tobjval} )  where  Tobjval
  return [
   				weightMinSum{Tobjval}(),   			# --> 1
  				weightMinProduct{Tobjval}(),   	# --> 2
  				weightMinMax{Tobjval}(),   			# --> 3
  				weightMaxSum{Tobjval}(),   		# --> 4
  				weightMaxProduct{Tobjval}(),   	# --> 5
  				weightMaxMin{Tobjval}()   			# --> 6
  			]
end

##############################################################
##########        weights that sum up   ( min = better  )  ###
##############################################################
struct weightMinSum{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMinSum{Tobjval}() where Tobjval
	   new(0)                                       ## default weight                    					 ### 0 must be convertible to Tobjval
	   				### setting it to 1 can prevent from zero weighted closed path !!!!
	 end
	 function weightMinSum{Tobjval}(nb) where Tobjval
	   new( convert(Tobjval,nb) )																				### nb must be convertible to Tobjval
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMinSum{Tobjval} ) where  Tobjval
   		return UInt8(1)
   end
############## objective type name and value (for trace)
   function descr(a::weightMinSum{Tobjval} ) where  Tobjval
   		return "min∑("*string(a.value)*")"
   end
   ############## initialisation value
	function defaultValue(a::weightMinSum{Tobjval} ) where  Tobjval
		w::weightMinSum{Tobjval} = weightMinSum{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMinSum{Tobjval}, value ) where  Tobjval
		return weightMinSum{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMinSum{Tobjval},b::weightMinSum{Tobjval}) where  Tobjval
	 	return ( a.value <= b.value, b.value <= a.value )  ####### Tobjval must implement <= #######
	 end
############## combination
	 function combine(a::weightMinSum{Tobjval},b::weightMinSum{Tobjval}) where  Tobjval
	 	return weightMinSum{Tobjval}(a.value + b.value) ####### Tobjval must implement + #######
	 end

##############################################################
###########        weights that multiply ( min = better  )    ################
##############################################################
struct weightMinProduct{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMinProduct{Tobjval}() where Tobjval
	   new(1)                                       ## defaults to one, all values assumed to be >= 1
	 end
	 function weightMinProduct{Tobjval}(nb::Number) where Tobjval
	   new( convert(Tobjval,nb) )
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMinProduct{Tobjval} ) where  Tobjval
   		return UInt8(2)
   end
############## objective type name and value (for trace)
   function descr(a::weightMinProduct{Tobjval} ) where  Tobjval
   		return "min∏("*string(a.value)*")"
   end
############## initialisation value
	function defaultValue(a::weightMinProduct{Tobjval} ) where  Tobjval
		w::weightMinProduct{Tobjval} = weightMinProduct{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMinProduct{Tobjval}, value ) where  Tobjval
		return weightMinProduct{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMinProduct{Tobjval},b::weightMinProduct{Tobjval}) where  Tobjval
	 	return ( a.value <= b.value, b.value <= a.value )
	 end
 ############## combination
	 function combine(a::weightMinProduct{Tobjval},b::weightMinProduct{Tobjval}) where  Tobjval
	 	return weightMinProduct{Tobjval}(a.value * b.value)
	 end

##############################################################
##########           bottleneck objectives  ( min = better  )       #############
##############################################################
struct weightMinMax{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMinMax{Tobjval}() where Tobjval
	   new(0)                                       ## begin with 0 (all values > 0)
	 end
	 function weightMinMax{Tobjval}(nb) where Tobjval
	   new( convert(Tobjval,nb) )
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMinMax{Tobjval} ) where  Tobjval
   		return UInt8(3)
   end
############## objective type name and value (for trace)
   function descr(a::weightMinMax{Tobjval} ) where  Tobjval
   		return "minMAX("*string(a.value)*")"
   end
############## initialisation value
	function defaultValue(a::weightMinMax{Tobjval} ) where  Tobjval
		w::weightMinMax{Tobjval} = weightMinMax{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMinMax{Tobjval}, value ) where  Tobjval
	 	return weightMinMax{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMinMax{Tobjval},b::weightMinMax{Tobjval}) where  Tobjval
	 	return  ( a.value <= b.value, b.value <= a.value )
	 end
############## combination
	 function combine(a::weightMinMax{Tobjval},b::weightMinMax{Tobjval}) where  Tobjval
	 	return weightMinMax{Tobjval}(max(a.value, b.value))
	 end

##############################################################
##########       weights that sum up   ( max = better  )  #################
##############################################################
struct weightMaxSum{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMaxSum{Tobjval}() where Tobjval
	   new(0)                                       ## default weight is zero
	 end
	 function weightMaxSum{Tobjval}(nb) where Tobjval
	   new( convert(Tobjval,nb) )
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMaxSum{Tobjval} ) where  Tobjval
   		return UInt8(4)
   end
############## objective type name and value (for trace)
   function descr(a::weightMaxSum{Tobjval} ) where  Tobjval
   		return "max∑("*string(a.value)*")"
   end
############## initialisation value
	function defaultValue(a::weightMaxSum{Tobjval} ) where  Tobjval
		w::weightMaxSum{Tobjval} = weightMaxSum{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMaxSum{Tobjval}, value ) where  Tobjval
		return weightMaxSum{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMaxSum{Tobjval},b::weightMaxSum{Tobjval}) where  Tobjval
	 	return ( b.value <= a.value, a.value <= b.value )
	 end
############## combination
	 function combine(a::weightMaxSum{Tobjval},b::weightMaxSum{Tobjval}) where  Tobjval
	 	return weightMaxSum{Tobjval}(a.value + b.value)
	 end

##############################################################
###########       weights that multiply ( max = better  )    ################
##############################################################
struct weightMaxProduct{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMaxProduct{Tobjval}() where Tobjval
	   new(1)                                       ## initialisation = 1 ( all values <= 1 -> probabilities)
	 end
	 function weightMaxProduct{Tobjval}(nb) where Tobjval
	   new( convert(Tobjval,nb) )
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMaxProduct{Tobjval} ) where  Tobjval
   		return UInt8(5)
   end
############## objective type name and value (for trace)
   function descr(a::weightMaxProduct{Tobjval} ) where  Tobjval
   		return "max∏("*string(a.value)*")"
   end
############## initialisation value
	function defaultValue(a::weightMaxProduct{Tobjval} ) where  Tobjval
		w::weightMaxProduct{Tobjval} = weightMaxProduct{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMaxProduct{Tobjval}, value ) where  Tobjval
		return weightMaxProduct{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMaxProduct{Tobjval},b::weightMaxProduct{Tobjval}) where  Tobjval
	 	return ( b.value <= a.value, a.value <= b.value )
	 end
 ############## combination
	 function combine(a::weightMaxProduct{Tobjval},b::weightMaxProduct{Tobjval}) where  Tobjval
	 	return weightMaxProduct{Tobjval}(a.value * b.value)
	 end

##############################################################
##########           bottleneck objectives  ( max = better )       #############
##############################################################
struct weightMaxMin{Tobjval}  <: weightCategory{Tobjval}
				# immutable  -> perhaps make it mutable to gain execution  time and memory
   value::Tobjval
   ############## constructors
	function weightMaxMin{Tobjval}() where Tobjval
	   new(typemax(Tobjval))                                       ## initialisation = infinity
	 end
	 function weightMaxMin{Tobjval}(nb::Number) where Tobjval
	   new( convert(Tobjval,nb) )
	 end
end
############## objective type code (for graph text file format)
   function typeCode(a::weightMaxMin{Tobjval} ) where  Tobjval
   		return UInt8(6)
   end
 ############## objective type name and value (for trace)
   function descr(a::weightMaxMin{Tobjval} ) where  Tobjval
   		return "maxMIN("*string(a.value)*")"
   end
 ############## initialisation value
	function defaultValue(a::weightMaxMin{Tobjval} ) where  Tobjval
		w::weightMaxMin{Tobjval} = weightMaxMin{Tobjval}()
		return w.value
	end
############## same type with another value
	function mimic( a::weightMaxMin{Tobjval}, value::Number ) where  Tobjval
	 	return weightMaxMin{Tobjval}( value )
	end
############## comparaison
	 function compare(a::weightMaxMin{Tobjval},b::weightMaxMin{Tobjval}) where  Tobjval
	 	return  ( b.value <= a.value, a.value <= b.value )
	 end
############## combination
	 function combine(a::weightMaxMin{Tobjval},b::weightMaxMin{Tobjval}) where  Tobjval
	 	return weightMaxMin{Tobjval}(min(a.value, b.value))
	 end

############ unit tests ###################

println("objectives loaded")
### TO DO #############

max_stars_repo_name: DRobbes/multiobjectiveGraphs

max_issues_repo_name: DRobbes/multiobjectiveGraphs

max_forks_repo_name: DRobbes/multiobjectiveGraphs


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

export BrickletIndustrialDualAnalogInV2Calibration
struct BrickletIndustrialDualAnalogInV2Calibration
    offset::Vector{Integer}
    gain::Vector{Integer}
end

export BrickletIndustrialDualAnalogInV2ChannelLEDStatusConfig
struct BrickletIndustrialDualAnalogInV2ChannelLEDStatusConfig
    min::Integer
    max::Integer
    config::Integer
end

export BrickletIndustrialDualAnalogInV2AllVoltagesCallbackConfiguration
struct BrickletIndustrialDualAnalogInV2AllVoltagesCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

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

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

export BrickletIndustrialDualAnalogInV2
"""
Measures two DC voltages between -35V and +35V with 24bit resolution each
"""
mutable struct BrickletIndustrialDualAnalogInV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage for the given channel.

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

export set_voltage_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Voltage` 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:`Voltage` 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_voltage_callback_configuration(device::BrickletIndustrialDualAnalogInV2, channel, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_voltage_callback_configuration(channel, period, value_has_to_change, option, min, max)
end

export get_voltage_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Voltage Callback Configuration`.
"""
function get_voltage_callback_configuration(device::BrickletIndustrialDualAnalogInV2, channel)
    return device.deviceInternal.get_voltage_callback_configuration(channel)
end

export set_sample_rate
"""
    $(SIGNATURES)

Sets the sample rate. The sample rate can be between 1 sample per second
and 976 samples per second. Decreasing the sample rate will also decrease the
noise on the data.
"""
function set_sample_rate(device::BrickletIndustrialDualAnalogInV2, rate)
    device.deviceInternal.set_sample_rate(rate)
end

export get_sample_rate
"""
    $(SIGNATURES)

Returns the sample rate as set by :func:`Set Sample Rate`.
"""
function get_sample_rate(device::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_sample_rate()
end

export set_calibration
"""
    $(SIGNATURES)

Sets offset and gain of MCP3911 internal calibration registers.

See MCP3911 datasheet 7.7 and 7.8. The Industrial Dual Analog In Bricklet 2.0
is already factory calibrated by Tinkerforge. It should not be necessary
for you to use this function
"""
function set_calibration(device::BrickletIndustrialDualAnalogInV2, offset, gain)
    device.deviceInternal.set_calibration(offset, gain)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_calibration()
end

export get_adc_values
"""
    $(SIGNATURES)

Returns the ADC values as given by the MCP3911 IC. This function
is needed for proper calibration, see :func:`Set Calibration`.
"""
function get_adc_values(device::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_adc_values()
end

export set_channel_led_config
"""
    $(SIGNATURES)

Each channel has a corresponding LED. You can turn the LED off, on or show a
heartbeat. You can also set the LED to "Channel Status". In this mode the
LED can either be turned on with a pre-defined threshold or the intensity
of the LED can change with the measured value.

You can configure the channel status behavior with :func:`Set Channel LED Status Config`.

By default all channel LEDs are configured as "Channel Status".
"""
function set_channel_led_config(device::BrickletIndustrialDualAnalogInV2, channel, config)
    device.deviceInternal.set_channel_led_config(channel, config)
end

export get_channel_led_config
"""
    $(SIGNATURES)

Returns the channel LED configuration as set by :func:`Set Channel LED Config`
"""
function get_channel_led_config(device::BrickletIndustrialDualAnalogInV2, channel)
    return device.deviceInternal.get_channel_led_config(channel)
end

export set_channel_led_status_config
"""
    $(SIGNATURES)

Sets the channel LED status config. This config is used if the channel LED is
configured as "Channel Status", see :func:`Set Channel LED Config`.

For each channel you can choose between threshold and intensity mode.

In threshold mode you can define a positive or a negative threshold.
For a positive threshold set the "min" parameter to the threshold value in mV
above which the LED should turn on and set the "max" parameter to 0. Example:
If you set a positive threshold of 10V, the LED will turn on as soon as the
voltage exceeds 10V and turn off again if it goes below 10V.
For a negative threshold set the "max" parameter to the threshold value in mV
below which the LED should turn on and set the "min" parameter to 0. Example:
If you set a negative threshold of 10V, the LED will turn on as soon as the
voltage goes below 10V and the LED will turn off when the voltage exceeds 10V.

In intensity mode you can define a range in mV that is used to scale the brightness
of the LED. Example with min=4V, max=20V: The LED is off at 4V, on at 20V
and the brightness is linearly scaled between the values 4V and 20V. If the
min value is greater than the max value, the LED brightness is scaled the other
way around.
"""
function set_channel_led_status_config(device::BrickletIndustrialDualAnalogInV2, channel, min, max, config)
    device.deviceInternal.set_channel_led_status_config(channel, min, max, config)
end

export get_channel_led_status_config
"""
    $(SIGNATURES)

Returns the channel LED status configuration as set by
:func:`Set Channel LED Status Config`.
"""
function get_channel_led_status_config(device::BrickletIndustrialDualAnalogInV2, channel)
    return device.deviceInternal.get_channel_led_status_config(channel)
end

export get_all_voltages
"""
    $(SIGNATURES)

Returns the voltages for all channels.

If you want to get the value periodically, it is recommended to use the
:cb:`All Voltages` callback. You can set the callback configuration
with :func:`Set All Voltages Callback Configuration`.

.. versionadded:: 2.0.6\$nbsp;(Plugin)
"""
function get_all_voltages(device::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_all_voltages()
end

export set_all_voltages_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Voltages`
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 at least one of the values has changed. If the values 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.

.. versionadded:: 2.0.6\$nbsp;(Plugin)
"""
function set_all_voltages_callback_configuration(device::BrickletIndustrialDualAnalogInV2, period, value_has_to_change)
    device.deviceInternal.set_all_voltages_callback_configuration(period, value_has_to_change)
end

export get_all_voltages_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set All Voltages Callback Configuration`.

.. versionadded:: 2.0.6\$nbsp;(Plugin)
"""
function get_all_voltages_callback_configuration(device::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_all_voltages_callback_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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2, 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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2, 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::BrickletIndustrialDualAnalogInV2, 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::BrickletIndustrialDualAnalogInV2, 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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2, 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::BrickletIndustrialDualAnalogInV2)
    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::BrickletIndustrialDualAnalogInV2)
    return device.deviceInternal.get_identity()
end

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

ext: jl


text:
# This file was generated, do not modify it. # hide
DataFrames.aggregate(iris, :Species, std)


text:
include("Oving13/Ov13.jl")

function main()

	println(certifysubsetsum_long([3, 5], [14, 23, 18, 2, 4], 22))
end

main()


text:
using Test

@testset "JLDArchives" begin
    include("v0.2.jl")
    include("v0.4.jl")
end


text:
module SimpleModel

using NLOptControl
using VehicleModels
using DataFrames
using Parameters

include("CaseModule.jl")
using .CaseModule

export
      initializeSimpleModel,
      updateSimpleModel,
      runSimpleModel,
      compareSimpleModels,
      resetDesignVariables

type OptData
      SA
      AX
      OBJ
      XO
      t
end

"""
initializeSimpleModels(n,c,pa);
# might have to call a single function and just have different kwargs
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function initializeSimpleModel(n,c,pa;x_min::Float64=0.0)  # can intially pass different x_min
error("this needs to be updated")
  # initialize
  s1=Settings();

  # define
  @unpack_Vpara pa # other vehicle model's inital parameters
  X0 = [x0_,y0_,psi0_,u0_];
  XF = [NaN,NaN,NaN,NaN];

  XL = [x_min,y_min,psi_min,u_min];
  XU = [x_max,y_max,psi_max,u_max];
  CL = [sa_min,-2.0];
  CU = [sa_max,2.5];
  define(KinematicBicycle;numStates=4,numControls=2,X0=X0,XF=XF,XL=XL,XU=XU,CL=CL,CU=CU)

  # build
  configure!(n1,Nck=n.Nck;(:integrationScheme => n.integrationScheme),(:finalTimeDV => true))
  defineSolver!(n1,solver=:IPOPT)
  error("update this!")
  mdl=defineSolver!(n,c);

  # setup OCP
  params1 = [VparaKB()];   # simple vehicle parameters TODO--> are they the same?
  r1=OCPdef!(mdl1,n1,s1,params1);
  @NLobjective(mdl1, Min,  n1.tf + (r1.x[end,1]-c.x_ref)^2 + (r1.x[end,2]-c.y_ref)^2);

  # obstacles
  Q = size(c.A)[1]; # number of obstacles
  @NLparameter(mdl1, a[i=1:Q] == c.A[i]);
  @NLparameter(mdl1, b[i=1:Q] == c.B[i]);
  @NLparameter(mdl1, X_0[i=1:Q] == c.X0_obs[i]);
  @NLparameter(mdl1, Y_0[i=1:Q] == c.Y0_obs[i]);
  @NLparameter(mdl1, speed_x[i=1:Q] == c.s_x[i]);
  @NLparameter(mdl1, speed_y[i=1:Q] == c.s_y[i]);

  # obstacle postion after the intial postion
  X_obs=@NLexpression(mdl1, [j=1:Q,i=1:n1.numStatePoints], X_0[j] + speed_x[j]*n1.tV[i]);
  Y_obs=@NLexpression(mdl1, [j=1:Q,i=1:n1.numStatePoints], Y_0[j] + speed_y[j]*n1.tV[i]);

  # constraint on position
  obs_con=@NLconstraint(mdl1, [j=1:Q,i=1:n1.numStatePoints-1], 1 <= ((r1.x[(i+1),1]-X_obs[j,i])^2)/((a[j]+c.sm)^2) + ((r1.x[(i+1),2]-Y_obs[j,i])^2)/((b[j]+c.sm)^2));
  newConstraint!(r1,obs_con,:obs_con);

  # solve
  optimize!(mdl1,n1,r1,s1)

  return mdl1, n1, r1
end

# TODO  http://docs.julialang.org/en/stable/manual/parallel-computing/#clustermanagers
#=
immutable LocalManager <: ClusterManager
    np::Integer
end

function launch(manager::LocalManager, params::Dict, launched::Array, c::Condition)
    ...
end

function manage(manager::LocalManager, id::Integer, config::WorkerConfig, op::Symbol)
    ...
end
=#
"""
# create several SimpleModels to be run in parallel
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function createParallelModels(n,c,pa)
#  http://docs.julialang.org/en/stable/manual/parallel-computing/#clustermanagers
  #addprocs(3);
  Xmin=[100.,180.,100.];
  Xmax=[300.,300.,220.];
  #sm1=remotecall(initializeSimpleModel,2,(n,c,pa;x_min=Xmin[1]));
  #sm2=remotecall(initializeSimpleModel,3,(n,c,pa;x_min=Xmin[2]));
  #sm3=remotecall(initializeSimpleModel,4,(n,c,pa;x_min=Xmin[3]));
  sm1=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[1],x_max=Xmax[1])
  sm2=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[2],x_max=Xmax[2])
  sm3=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[3],x_max=Xmax[3])

  sm=[sm1,sm2,sm3];
  return sm
end

"""
# based off of predicted vehicle the states of a SimpleModel
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function updateSimpleModel(r,sm)

end

"""
# based off of predicted vehicle the states of a SimpleModel
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function runSimpleModel(r)

end

"""
# based off of predicted vehicle the states of a SimpleModel
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function compareSimpleModels(r)

end

"""
# use the best output from the SimpleModels update certain DesignVariables
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""

function resetDesignVariables()

end

end # module


text:
#  C:\dev\GitHub\autodiff\logsumexp.jl
#  AUTOGENERATED FROM C:\dev\GitHub\autodiff\logsumexp.ipynb on 2015-08-03T13:16:36
using Base.Test

# logsumexp()
function logsumexp(x::Array{Float64,1})
  A = maximum(x);
  ema = exp(x-A);
  sema = sum(ema);
  log(sema) + A;
end

r = rand(5);
@test_approx_eq_eps logsumexp(r) log(sum(exp(r))) 1.0e-8
@printf("logsumexp: %f = %f\n", logsumexp(r), log(sum(exp(r))));

# logsumexp()
function logsumexp_both(x::Array{Float64,1})
  A = maximum(x);
  ema = exp(x-A);
  sema = sum(ema);
  l = log(sema) + A;
  Jacobian = ema/sema;
  return (l, Jacobian);
end

#logsumexp(x::Array{Float64,1}) = logsumexp_both(x)[1]

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
"""

FWD (Glaescher et al. 2010)
---------------------------

Forward Learner

See learner/forwardlearner.jl for the core implementation.

"""
function getLL_Forward_Learner(data::Array{Int,2}, eta::Float64,
                                gamma::Float64, temp::Float64;
                                initTsas = nothing, initR = nothing)

    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_Forward_Learner supports only one subject")
    end

    LL_unbiased = 0.
    action_selection_probs = zeros(Float64, nr_samples)
    # signals:
    SPE_record = zeros(Float64, nr_samples)
    RiskPE_record = zeros(Float64, nr_samples)
    Qsa_fwd_record = zeros(Float64, nr_samples)

    final_TSAS_record = Dict{Int,Any}()
    final_Reward_record = Dict{Int,Any}()
    current_puzzle_id = -1
    # Puzzle = independent block of the experiment.
    # The code was made for general purposes, in this experiment we
    # have only ONE puzzle.
    learner = []
    for i = 1:nr_samples-1
        if current_puzzle_id != data[i,Int(c_puzzle)]
            current_puzzle_id = data[i,Int(c_puzzle)]
            # new puzzle. reset learning
            if initTsas == nothing
                learner = ForwardLearner(nr_states = nr_states,
                                    nr_actions = nr_actions,
                                    η = eta, γ = gamma,
                                    init_R = 0.,
                                    )
            else
                learner = ForwardLearner(nr_states = nr_states,
                                    nr_actions = nr_actions,
                                    η = eta, γ = gamma,
                                    TSAS = copy(initTsas[current_puzzle_id]), # first argument, i-th puzzle
                                    R = copy(initR[current_puzzle_id])
                                    )
            end
        end
        state = data[i, Int(c_state)]
        action = data[i, Int(c_action)]
        reward = data[i, Int(c_reward)]
        next_state = data[i, Int(c_nextState)]
        Qsa_fwd_record[i] = copy(learner.QSA[state,action])

        p_actions = getSoftmaxProbs(learner.QSA[state,:], temp, [0., 0.])
        action_selection_probs[i] = copy(p_actions[action])
        LL_unbiased += log(p_actions[action])

        # ForwardLearner
        update!(learner, state, action, next_state, reward)
        SPE_record[i] = copy(learner.spe[1])

        final_TSAS_record[current_puzzle_id] = copy(learner.TSAS)
        final_Reward_record[current_puzzle_id] = copy(learner.R)
    end

    results_dict = Dict{String,Any}(
        "QSA_record"=> Qsa_fwd_record,
        "final_Tsas" => final_TSAS_record,
        "final_Qsa" => learner.QSA,
        "final_Reward" => final_Reward_record
        )

    LL = LL_unbiased
    trial_by_trial_p_sa = (action_selection_probs)
    algo = "Forward Learner"
    param_v = [eta, gamma, temp]
    parma_str = "eta=$(round(eta,digits=4)), gamma=$(round(gamma,digits=4)),
                temp=$(round(temp,digits=4)),
                initTsas=$(repr(initTsas)), initR=$(repr(initR))"
    signals = Step_by_step_signals(
        data_stats.all_subject_ids[1],
        [], #RPE
        SPE_record,
        trial_by_trial_p_sa,
        LL,
        data_stats,
        algo,
        param_v,
        parma_str,
        99, #biased_action_idx, dummy value
        results_dict
        )
    return LL, signals
end


text:
import LinearAlgebra: norm, normalize, dot
import Random: seed!

@test invlerp(0.0,1.0,0.5) ≈ 0.5
@test invlerp(0.0,1.0,0.4) ≈ 0.4
@test invlerp(0.0,1.0,0.7) ≈ 0.7
@test invlerp(10.0,20.0,12.0) ≈ 0.2

@testset "VecE2" begin
    include("test_vecE2.jl")
end

@testset "VecE3" begin
    include("test_vecE3.jl")
end

@testset "VecSE2" begin
    include("test_vecSE2.jl")
end

@testset "quaternions" begin
    include("test_quat.jl")
end

@testset "coordinate transforms" begin
    include("test_coordinate_transforms.jl")
end

@testset "geom" begin
    include("test_geomE2.jl")
end

@testset "diff" begin
    include("test_diff.jl")
end


text:
## PART 7

include("autocorrelation.jl");

println("\nPart 7");
println("\n\tPick subjets v1 and v2 of each dataset,");
println("\tsizes N1 = 1,000 & N2 = 10,000.");

v1 = v[[1:N1],:];
v2 = v[[1:N2],:];

println("\n\tEstimate autocorrelations");
println("\n\t...for N1...");
xcorr = [0:10:400];
ycorr1 = zeros( length( xcorr ), 5 );
for a in [1:5]
  # Calculate autocorrelations
  ycorr1[:,a] = autocorrelation(N1, v1, a, xcorr);
  # Normalize by dividing through by Ĉ[1,a]
  ycorr1[:,a] /= ycorr1[1,a];
end
println("\n\t...for N2 (See plots below)");
ycorr2 = zeros( length( xcorr ), 5 );
for a in [1:5]
  # Calculate autocorrelations
  ycorr2[:,a] = autocorrelation(N2, v2, a, xcorr);
  # Normalize by dividing through by Ĉ[1,a]
  ycorr2[:,a] /= ycorr2[1,a];
end

# Make plots for each alpha
println("\n\tPlacing estimated autocorrelation plots 'plots1' and 'plots2':");
using Gadfly
plots1 = Plot[];
plots2 = Plot[];
for i in [1:5]
  # Plot
  push!(plots1, Gadfly.plot ( x = xcorr, y = ycorr1[:,i], Geom.line, Geom.point,
      Guide.xlabel("Separation"), Guide.ylabel("Normalized Autocorrelation"),
      Guide.title(join(["Estimated Autocorrelation for Dataset",string(i),
                        "with Sample Size N1 = $N1"]," "," "))))
  push!(plots2, Gadfly.plot ( x = xcorr, y = ycorr2[:,i], Geom.line, Geom.point,
      Guide.xlabel("Separation"), Guide.ylabel("Normalized Autocorrelation"),
      Guide.title(join(["Estimated Autocorrelation for Dataset",string(i),
                        "with Sample Size N2 = $N2"]," "," "))))
end

# Estimate autocorrelation times
println("\n\tEstimating autocorrelation times τ1 and τ2");
τ1 = Float64[];
τ2 = Float64[];
for a in [1:5]
  corra = autocorrelation(N1, v1, a, ns);
  corra /= corra[1];                          # Normalize
  push!(τ1,sum(corra) - 0.5);                 # Get autocorrelation time
end
for a in [1:5]
  corra = autocorrelation(N2, v2, a, ns);
  corra /= corra[1];                          # Normalize
  push!(τ2,sum(corra) - 0.5);                 # Get autocorrelation time
end

println("\n\tEstimated Autocorrelation Times for N1 = $N1:\n");
for i in [1:5]
  println("\t\tτ̂1v$i:\t$(τ1[i])");
end
println("\n\tEstimated Autocorrelation Times for N2 = $N2:\n");
for i in [1:5]
  println("\t\tτ̂2v$i:\t$(τ2[i])");
end
println("\n\tThe estimated autocorrelation times are clearly far better for the larger sample.");

# Estimate standard deviation
σ1 = Float64[];
σ2 = Float64[];
println("\n\tCreating estimators σ1 for the standard deviation for N1 = $N1");
for a in [1:5]
  diff1 = v1[:,a] .- mean(v1[:,a]);
  σ1a = sqrt( dot( diff1, diff1 ) / (N1-1) );
  println("\t\tσ̂1[$a] = $σ1a");
  push!(σ1, σ1a);
end
println("\n\tCreating estimators σ2 for the standard deviation for N2 = $N2");
for a in [1:5]
  diff2 = v2[:,a] .- mean(v2[:,a]);
  σ2a = sqrt( dot( diff2, diff2 ) / (N2-1) );
  println("\t\tσ̂2[$a] = $σ2a");
  push!(σ2, σ2a);
end

# Estimate standard deviation of the sample mean
σN1 = sqrt( 2 / N1 .* τ1 ) .* σ1;
σN2 = sqrt( 2 / N2 .* τ2 ) .* σ2;
println("\n\tEstimating standard deviations of the sample mean for N1 = $N1:");
for i in [1:5]
  println("\n\t\tσ̂N1[$i] = $(σN1[i])");
end
println("\n\tEstimating standard deviations of the sample mean for N2 = $N2:");
for i in [1:5]
  println("\n\t\tσ̂N2[$i] = $(σN2[i])");
end

# Compare to the actual measured values of the standard deviations of the sample mean
sigmean1 = σ̂1./σN1;
sigmean2 = σ̂2./σN2;
println("\n\tTake ratio of measured standard deviation of sample means over the estimated");
println("\tstandard deviation just found in order to measure the accuracy of the estimate.");
println("\n\tFind the ratio for N1 = $N1:")
for i in [1:5]
  println("\n\t\tσ̂1[$i] / σ̂N1[$i] = $(sigmean1[i])");
end
println("\n\tFind the ratio for N2 = $N2:")
for i in [1:5]
  println("\n\t\tσ̂2[$i] / σ̂N2[$i] = $(sigmean2[i])");
end

println("\n\tFor both sample sizes, these estimates are close to the true values of");
println("\tthe standard deviation of the sample means, though the agreement is far");
println("\tbetter for the larger sample size.");

# Find sample means
v̄1 = zeros(5);
v̄2 = zeros(5);
for i in [1:5]
  v̄1[i] = mean(v1[:,i]);
  v̄2[i] = mean(v2[:,i]);
end

# Finally, compare the normalized covariance matrices ρ1 and ρ2 to the true value ρ̂.
println("\n\tFinally, estimate the normalized covariance matrices ρ1 and ρ2.\n");
println("\tEstimating the covariance matrix c1 for sample size N1 = $N1");
D1 = v1 - ones(N1) * transpose(v̄1);
c1 = transpose(D1) * D1 / N1;
println("\tEstimating the covariance matrix c2 for sample size N2 = $N2");
D2 = v2 - ones(N2) * transpose(v̄2);
c2 = transpose(D2) * D2 / N2;

println("\tEstimating the normalized covariant matrix ρ1 for N1 = $N1");
ρ1 = c1 ./ ( σ1 * transpose(σ1) );
println("\tEstimating the normalized covariant matrix ρ2 for N2 = $N2");
ρ2 = c2 ./ ( σ2 * transpose(σ2) );


text:
# Tests for our auditing infrastructure

@testset "Auditor - ISA tests" begin
    begin
        build_path = tempname()
        mkpath(build_path)
        isa_platform = Linux(:x86_64; compiler_abi=CompilerABI(:gcc8, :cxx_any))
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], isa_platform)

        main_sse = ExecutableProduct(prefix, "main_sse", :main_sse)
        main_avx = ExecutableProduct(prefix, "main_avx", :main_avx)
        main_avx2 = ExecutableProduct(prefix, "main_avx2", :main_avx2)
        products = [main_sse, main_avx, main_avx2]

        cd(joinpath(dirname(@__FILE__),"build_tests","isa_tests")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # Build isa tests
            script="""
            /usr/bin/make clean
            /usr/bin/make install
            """

            # Build it
            @test build(ur, "isa_tests", products, script, isa_platform, prefix; verbose=true)

            # Ensure it's satisfied
            @test all(satisfied(r; verbose=true) for r in products)
        end

        # Next, test isa of these files
        readmeta(locate(main_sse)) do oh
            isa_sse = BinaryBuilder.analyze_instruction_set(oh, isa_platform; verbose=true)
            @test isa_sse == :core2
        end

        readmeta(locate(main_avx)) do oh
            isa_avx = BinaryBuilder.analyze_instruction_set(oh, isa_platform; verbose=true)
            @test isa_avx == :sandybridge
        end

        readmeta(locate(main_avx2)) do oh
            isa_avx2 = BinaryBuilder.analyze_instruction_set(oh, isa_platform; verbose=true)
            @test isa_avx2 == :haswell
        end

        # Delete the build path
        rm(build_path, recursive = true)
    end
end

@testset "Auditor - .dll moving" begin
    begin
        build_path = tempname()
        mkpath(build_path)
        dll_platform = Windows(:x86_64)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], dll_platform)
        cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # First, build libfoo, but with a dumb script that doesn't know to put .dll files in bin
            dumb_script = """
            /usr/bin/make clean
            /usr/bin/make install libdir=\$prefix/lib
            """

            @test build(ur, "foo", libfoo_products(prefix), libfoo_make_script, dll_platform, prefix; autofix=false)
        end

        # Test that libfoo puts its .dll's into lib, even on windows:
        @test !isfile(joinpath(prefix, "bin", "libfoo.dll"))
        @test isfile(joinpath(prefix, "lib", "libfoo.dll"))

        # Test that `audit()` moves it to `bin`.
        BinaryBuilder.audit(prefix; platform=dll_platform, verbose=true, autofix=true)
        @test isfile(joinpath(prefix, "bin", "libfoo.dll"))
        @test !isfile(joinpath(prefix, "lib", "libfoo.dll"))
    end
end

@testset "Auditor - .dylib identity mismatch" begin
    begin
        build_path = tempname()
        mkpath(build_path)
        dylib_platform = MacOS()
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], dylib_platform)
        cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # First, build libfoo, but with a dumb script that doesn't know to put .dll files in bin
            dumb_script = """
            /usr/bin/make clean
            /usr/bin/make install libdir=\$prefix/lib
            """

            @test build(ur, "foo", libfoo_products(prefix), libfoo_make_script, dylib_platform, prefix; autofix=false)
        end

        function get_dylib_id(path)
            return readmeta(path) do oh
                return dylib_name(first([lc for lc in MachOLoadCmds(oh) if typeof(lc) <: MachOIdDylibCmd]))
            end
        end
        # Test that libfoo gives an incorrect dylib ID:
        @test !startswith(get_dylib_id(joinpath(prefix, "lib", "libfoo.dylib")), "@rpath")

        # Test that `audit()` changes the
        BinaryBuilder.audit(prefix; platform=dylib_platform, verbose=true, autofix=true)
        @test startswith(get_dylib_id(joinpath(prefix, "lib", "libfoo.dylib")), "@rpath")
    end
end

@testset "Auditor - absolute paths" begin
    prefix = Prefix(tempname())
    try
        sharedir = joinpath(prefix.path, "share")
        mkpath(sharedir)
        open(joinpath(sharedir, "foo.conf"), "w") do f
            write(f, "share_dir = \"$sharedir\"")
        end

        # Test that `audit()` warns about an absolute path within the prefix
        @test_warn "share/foo.conf" BinaryBuilder.audit(prefix)
    finally
        rm(prefix.path; recursive=true)
    end
end

@testset "Auditor - gcc version" begin
    begin
        build_path = tempname()

        # These tests assume our gcc version is concrete
        our_gcc_version = BinaryProvider.compiler_abi(platform_key_abi()).gcc_version
        @test our_gcc_version != :gcc_any

        # Get one that isn't us.
        other_gcc_version = :gcc4
        if our_gcc_version == other_gcc_version
            other_gcc_version = :gcc7
        end

        our_platform = platform_key_abi()
        other_platform = BinaryBuilder.replace_gcc_version(our_platform, other_gcc_version)
        fail_platform = BinaryBuilder.replace_gcc_version(our_platform, :gcc_any)

        # Build `hello` for our own platform, ensure it builds and runs:
        mkpath(build_path)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], our_platform)
        hello = ExecutableProduct(prefix, "hello", :hello)

        cd(joinpath(dirname(@__FILE__),"build_tests","fortran_hello")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # Build it
            @test build(ur, "hello", [hello], "make install", our_platform, prefix; verbose=true)

            # Ensure it's satisfied
            @test satisfied(hello; verbose=true)
        end

        # Ensure that we can actually run this `hello`:
        withenv(prefix) do
            @test success(`hello`)
        end
        rm(build_path, recursive = true)

        # Now, build one with not-our-GCC:
        mkpath(build_path)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], other_platform)
        hello = ExecutableProduct(prefix, "hello", :hello)
        cd(joinpath(dirname(@__FILE__),"build_tests","fortran_hello")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # Build it
            @test build(ur, "hello", [hello], "make install", other_platform, prefix; verbose=true)

            # Ensure it's satisfied (because `hello` exists, even if it can't run)
            @test satisfied(hello; verbose=true)
        end

        # Ensure that we can't actually run this `hello`:
        withenv(prefix) do
            @test !success(`hello`)
        end
        rm(build_path, recursive = true)

        # Finally, build with `:gcc_any`, which should throw because it needs to be gfortran-specific
        mkpath(build_path)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], fail_platform)
        hello = ExecutableProduct(prefix, "hello", :hello)
        cd(joinpath(dirname(@__FILE__),"build_tests","fortran_hello")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # Build it
            @test_warn "links to libgfortran!" !build(ur, "hello", [hello], "make install", fail_platform, prefix; verbose=true)
        end
        rm(build_path, recursive = true)
    end
end


text:
struct Params_S2{ComplexArray,Total}
    # boundary conditions
    periodic::Bool
    # normalization
    total::Total

    # algorithm-specific

    # fft buffers
    complex_img::ComplexArray
end

function Params_S2(img; periodic::Bool=true)
    box = size(img)
    complex_box = periodic ? box : box .* 2

    total = cnt_total(img, periodic)

    p = Params_S2(
        periodic,
        total,
        similar(img, ComplexF32, complex_box),
    )
    cf_type = periodic ? :periodic_point_point : :central_symmetry
    p, cf_type
end

"""
    correllation_function!(res, img, params::Params_S2)

Compute S₂ correlation function in positive directions
"""
function correllation_function!(res, img, params::Params_S2)
    ix = CartesianIndices(img)

    f = params.complex_img .= 0

    v_f = view(f, ix)
    v_f .= img

    self_correlation!(f)

    res .= real.(v_f) ./ params.total
end

"""
    s2(image; periodic = false)

Calculate `S₂` (two point) correlation function map
for the N-dimensional image and return a `CFMap` object.

The `image` contains the probability of the voxel being in the correct phase.

# Examples
```jldoctest
julia> s2([1 0; 0 1]; periodic=true).result
2×2 Matrix{Float32}:
 0.5  0.0
 0.0  0.5
```
"""
function s2(image; periodic::Bool=false)
    corr_function_map(image, Params_S2; periodic)
end

function s2(image, phase::Int; periodic::Bool=false)
    one_phase_img = image .== phase
    corr_function_map(one_phase_img, Params_S2; periodic)
end


text:
"""

Union of types to define an `isaprox` for testing.

"""
const CMTypes = Union{Result,Density,Volume,SolSummary,Options}

"""

```
Base.isapprox(r1::T, r2::T; debug=false) where T <: CMTypes
```

Function to test if two runs offered similar results. Mostly used in the
package testing routines

"""
function Base.isapprox(r1::T, r2::T; debug = false) where {T<:CMTypes}
    check = true
    diff_list = Symbol[]
    for field in fieldnames(T)
        if field == :files
            continue
        end
        x = getfield(r1, field)
        y = getfield(r2, field)
        if !isapprox(x, y)
            check = false
            if debug
                push!(diff_list, field)
            end
        end
    end
    if debug
        for field in diff_list
            println(" Data in $field field differ. ")
        end
    end
    return check
end


text:
# Autogenerated wrapper script for ColPack_jll for x86_64-apple-darwin
export libcolpack

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ColPack")
JLLWrappers.@declare_library_product(libcolpack, "@rpath/libColPack.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcolpack,
        "lib/libColPack.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Oceananigans.Architectures
using Oceananigans.BoundaryConditions
using Oceananigans.TurbulenceClosures: calculate_diffusivities!

import Oceananigans.TimeSteppers: update_state!

"""
    update_state!(model::IncompressibleModel)

Update peripheral aspects of the model (halo regions, diffusivities, hydrostatic pressure) to the current model state.
"""
function update_state!(model::IncompressibleModel)

    # Fill halos for velocities and tracers
    fill_halo_regions!(merge(model.velocities, model.tracers), model.architecture,  model.clock, fields(model))

    # Calculate diffusivities
    calculate_diffusivities!(model.diffusivities, model.architecture, model.grid, model.closure,
                             model.buoyancy, model.velocities, model.tracers)

    fill_halo_regions!(model.diffusivities, model.architecture, model.clock, fields(model))

    # Calculate hydrostatic pressure
    pressure_calculation = launch!(model.architecture, model.grid, :xy, update_hydrostatic_pressure!,
                                   model.pressures.pHY′, model.grid, model.buoyancy, model.tracers,
                                   dependencies=Event(device(model.architecture)))

    # Fill halo regions for pressure
    wait(device(model.architecture), pressure_calculation)

    fill_halo_regions!(model.pressures.pHY′, model.architecture)

    return nothing
end


text:
using IndirectTrajOpt
using Heuristics
using IndirectCoStateInit
using IndirectShooting
using StaticArrays

# Initialize BVP function w/o STM
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
tspan = (0.0, 8.6404*24*3600/ps.crp.TU)
bvpFunc(y0) = cr3bpOptIntegrate(y0, tspan, ps, copyParams = true, termCallbacks = true)

# Initialize BVP function w/ STM
psSTM = initCR3BPIndirectWithSTMParams("Low Thrust 10 CR3BP")
bvpFuncNoSTM(z0) = cr3bpOptIntegrate(z0, tspan, ps, SingleOutput(); copyParams = true)
bvpFuncWSTM(z0, ϵ) = cr3bpOptWithSTMIntegrate(z0, tspan, ϵ, psSTM; copyParams = true)

# Set initial and final conditions
ics = [-0.019488511458668, -0.016033479812051, 0.0,
        8.918881923678198, -4.081793688818725, 0.0,
        1.0]
fcs = [ 0.823385182067467, 0.0, -0.022277556273235,
        0.0, 0.134184170262437, 0.0, 0.0]

# Initialize FFS Initializer
csInitializer = FSSCoStateInitializer(bvpFunc, ics, fcs;
                                    optimizer = :MS_PSO,
                                    numParticles = 1000,
                                    initMethod = :Uniform,
                                    displayInterval = 5,
                                    maxStallIters = 50,
                                    #UBs = [40, 40, 40, 2, 2, 2, 2],
                                    #LBs = [-40, -40, -40, -2, -2, -2, -2],
                                    UBs = [100, 100, 100, 10, 10, 10, 10],
                                    LBs = [-100, -100, -100, -10, -10, -10, -10],
                                    iUBs = [40, 40, 1, 2, 2, 2, 2],
                                    iLBs = [-40, -40, -1, -2, -2, -2, -2]
                                    );

# Initialize co-states
initialize!(csInitializer)

# Initialize Shooting Solver
solver = FSSSolver(GetInitializedCostates(csInitializer), ics, fcs, bvpFuncNoSTM, bvpFuncWSTM;
    homotopy = true, homotopyParamVec = [(j^2 - 1)/(25^2 - 1) for j in 25:-1:1])

solve!(solver)


text:
import AlphaZero.GI

using Crayons
using StaticArrays

const NUM_COLS = 7
const NUM_ROWS = 6
const NUM_CELLS = NUM_COLS * NUM_ROWS
const TO_CONNECT = 4

const Player = UInt8
const WHITE = 0x01
const BLACK = 0x02

other(p::Player) = 0x03 - p

const Cell = UInt8
const EMPTY = 0x00
const Board = SMatrix{NUM_COLS, NUM_ROWS, Cell, NUM_CELLS}

const INITIAL_BOARD = @SMatrix zeros(Cell, NUM_COLS, NUM_ROWS)

mutable struct Game <: GI.AbstractGame
  board :: Board
  curplayer :: Player
  finished :: Bool
  winner :: Player
  amask :: Vector{Bool} # actions mask
  # Actions history, which uniquely identifies the current board position
  # Used by external solvers
  history :: Union{Nothing, Vector{Int}}
end

function Game()
  board = INITIAL_BOARD
  curplayer = WHITE
  finished = false
  winner = 0x00
  amask = trues(NUM_COLS)
  history = Int[]
  Game(board, curplayer, finished, winner, amask, history)
end

GI.Board(::Type{Game}) = Board
GI.Action(::Type{Game}) = Int

const ACTIONS = collect(1:NUM_COLS)

GI.actions(::Type{Game}) = ACTIONS

function Base.copy(g::Game)
  history = isnothing(g.history) ? nothing : copy(g.history)
  Game(g.board, g.curplayer, g.finished, g.winner, copy(g.amask), history)
end

history(g::Game) = g.history

#####
##### Defining game rules
#####

function first_free(board, col)
  row = 1
  while row <= NUM_ROWS && board[col, row] != EMPTY
    row += 1
  end
  return row
end

function update_actions_mask!(g::Game)
  g.amask = map(ACTIONS) do col
    first_free(g.board, col) <= NUM_ROWS
  end
end

GI.actions_mask(g::Game) = g.amask

valid_pos((col, row)) = 1 <= col <= NUM_COLS && 1 <= row <= NUM_ROWS

function num_connected_dir(board, player, pos, dir)
  @assert board[pos...] == player
  p = pos .+ dir
  n = 0
  while valid_pos(p) && board[p...] == player
    n += 1
    p = p .+ dir
  end
  return n
end

function num_connected_axis(board, player, pos, axis)
  @assert board[pos...] == player
  num_after = num_connected_dir(board, player, pos, axis)
  num_before = num_connected_dir(board, player, pos, (0, 0) .- axis)
  return 1 + num_before + num_after
end

function winning_pattern_at(board, player, pos)
  return any(((1, 1), (1, -1), (1, 0), (0, 1))) do axis
    num_connected_axis(board, player, pos, axis) >= TO_CONNECT
  end
end

# Update the game status assuming g.curplayer just played at pos=(col, row)
function update_status!(g::Game, pos)
  update_actions_mask!(g)
  if winning_pattern_at(g.board, g.curplayer, pos)
    g.winner = g.curplayer
    g.finished = true
  else
    g.finished = !any(g.amask)
  end
end

function GI.play!(g::Game, col)
  isnothing(g.history) || push!(g.history, col)
  row = first_free(g.board, col)
  g.board = setindex(g.board, g.curplayer, col, row)
  update_status!(g, (col, row))
  g.curplayer = other(g.curplayer)
end

function Game(board::Board, white_playing=true)
  g = Game()
  g.history = nothing
  g.board = board
  g.curplayer = white_playing ? WHITE : BLACK
  update_actions_mask!(g)
  any(g.amask) || (g.finished = true)
  for col in 1:NUM_COLS
    top = first_free(g.board, col)
    top == 1 && continue
    row = top - 1
    c = board[col, row]
    if c != EMPTY && winning_pattern_at(board, c, (col, row))
      g.winner = c
      g.finished = true
      break
    end
  end
  return g
end

GI.board(g::Game) = g.board

cell_symmetric(c::Cell) = c == EMPTY ? EMPTY : other(c)

function GI.board_symmetric(g::Game) :: Board
  return @SMatrix [
    cell_symmetric(g.board[col, row])
    for col in 1:NUM_COLS, row in 1:NUM_ROWS]
end

GI.white_playing(g::Game) = g.curplayer == WHITE

#####
##### Reward shaping
#####

function GI.white_reward(g::Game)
  if g.finished
    g.winner == WHITE && (return  1.)
    g.winner == BLACK && (return -1.)
    return 0.
  else
    return nothing
  end
end

#####
##### Simple heuristic for minmax
#####

const Pos = Tuple{Int, Int}
const Alignment = Vector{Pos}

function alignment_from(pos, dir) :: Union{Alignment, Nothing}
  al = Alignment()
  for i in 1:TO_CONNECT
    valid_pos(pos) || (return nothing)
    push!(al, pos)
    pos = pos .+ dir
  end
  return al
end

function alignments_with(dir) :: Vector{Alignment}
  als = [alignment_from((x, y), dir) for x in 1:NUM_COLS for y in 1:NUM_ROWS]
  return filter(al -> !isnothing(al), als)
end

const ALIGNMENTS = [
  alignments_with((1,  1));
  alignments_with((1, -1));
  alignments_with((0,  1));
  alignments_with((1,  0))]

function alignment_value_for(g::Game, player, alignment)
  γ = 0.1
  N = 0
  for pos in alignment
    cell = g.board[pos...]
    if cell == player
      N += 1
    elseif cell == other(player)
      return 0.
    end
  end
  return γ ^ (TO_CONNECT - 1 - N)
end

function heuristic_value_for(g::Game, player)
  return sum(alignment_value_for(g, player, al) for al in ALIGNMENTS)
end

function GI.heuristic_value(g::Game)
  mine = heuristic_value_for(g, g.curplayer)
  yours = heuristic_value_for(g, other(g.curplayer))
  return mine - yours
end

#####
##### ML interface
#####

function GI.vectorize_board(::Type{Game}, board)
  return Float32[
    board[col, row] == c
    for col in 1:NUM_COLS,
        row in 1:NUM_ROWS,
        c in [EMPTY, WHITE, BLACK]]
end

#####
##### Symmetries
#####

function flipped_board(board)
  return @SMatrix[board[col, row]
    for col in reverse(1:NUM_COLS), row in 1:NUM_ROWS]
end

function GI.symmetries(::Type{Game}, board)
  symb = flipped_board(board)
  σ = reverse(collect(1:NUM_COLS))
  return [(symb, σ)]
end

#####
##### User interface
#####

GI.action_string(::Type{Game}, a) = string(a)

function GI.parse_action(g::Game, str)
  try
    p = parse(Int, str)
    1 <= p <= NUM_COLS ? p : nothing
  catch
    nothing
  end
end

# 1 2 3 4 5 6 7
# . . . . . . .
# . . . . . . .
# . . . . . . .
# . . o x . . .
# . o o o . . .
# o x x x . x .

player_color(p) = p == WHITE ? crayon"light_red" : crayon"light_blue"
player_name(p)  = p == WHITE ? "Red" : "Blue"
player_mark(p)  = p == WHITE ? "o" : "x"
cell_mark(c)    = c == EMPTY ? "." : player_mark(c)
cell_color(c)   = c == EMPTY ? crayon"" : player_color(c)

function GI.print_state(g::Game; with_position_names=true, botmargin=true)
  pname = player_name(g.curplayer)
  pcol = player_color(g.curplayer)
  print(pcol, pname, " plays:", crayon"reset", "\n\n")
  # Print legend
  for col in 1:NUM_COLS
    print(GI.action_string(Game, col), " ")
  end
  print("\n")
  # Print board
  for row in NUM_ROWS:-1:1
    for col in 1:NUM_COLS
      c = g.board[col, row]
      print(cell_color(c), cell_mark(c), crayon"reset", " ")
    end
    print("\n")
  end
  botmargin && print("\n")
end

function GI.read_state(::Type{Game})
  board = Array(INITIAL_BOARD)
  try
    for col in 1:NUM_COLS
      input = readline()
      for (row, c) in enumerate(input)
        c = lowercase(c)
        if c ∈ ['o', 'w', '1']
          board[col, row] = WHITE
        elseif c ∈ ['x', 'b', '2']
          board[col, row] = BLACK
        end
      end
    end
    nw = count(==(WHITE), board)
    nb = count(==(BLACK), board)
    if nw == nb
      wp = true
    elseif nw == nb + 1
      wp = false
    else
      return nothing
    end
    return Game(Board(board), wp)
  catch e
    return nothing
  end
end


text:
module CovidSurvey

using DrWatson
using MCMCChains
using Distributions
using Plots
using StatsBase
using Dates
using DistributionsAD
using Random: AbstractRNG
import Bijectors
using CSV
using Underscores
using DataFrames
using Base:@kwdef

abstract type ObservationInit end
abstract type ObservationsModel end

export
    #National,
    Regional,
    International,
    LogNormalMeanStd,
    NegativeBinomial2,
    NegativeBinomial3,
    KLogistic,
    KLogit,
    InverseGamma2,
    GammaMeanCv,
    GammaMeanStd,
    vectup2tupvec,
    arrarrarr2arr,
    plot_confidence_timeseries,
    plot_confidence_timeseries!,
    hpdi

include("io.jl")
include("utils.jl")
include("model_utils.jl")           # <= stuff that might also be included by sub-modules
include("visualizations.jl")   # <= visualization stuff
include("data.jl")

# Simulations on a national level
#include("national/National.jl")
include("regional/Regional.jl")
include("international/International.jl")

__precompile__(true)
end # module

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# This file was generated, do not modify it. # hide
model1 = linear
model2 = knn


text:
using PyPlot
using ODE
using Dierckx
using Roots

setprecision(128);

include("cosmology.jl")
include("rays.jl")
include("mass.jl")
include("harmonic.jl")
include("numeric.jl");
include("thin.jl")
include("time.jl")
include("precision.jl")
include("defaultmodels.jl")
include("fermat.jl")

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
#
#--------------------------------------#
export mass
#--------------------------------------#
"""
 (v,u) = (Q*R'*v)' * B * (Q*R'*u)

 implemented as

 (QQ' * R'R * B_loc) * u_loc
"""
function mass(u::Array
             ,msh::Mesh)

    @unpack B = msh

    Bu = B .* u
    #Bu = gatherScatter(Bu,QQtx,QQty)
    #Bu = mask(Bu,M)

return Bu
end
#--------------------------------------#
# Dealiased implementation
function mass(u::Array
             ,msh1::Mesh
             ,msh2::Mesh)

    return mass(u,msh1)
end
#--------------------------------------#
function mass(u,M,B,Jr,Js,QQtx,QQty
             ,mult);

Ju = ABu(Js,Jr,u);

if(length(B)==0); BJu =      Ju;
else              BJu = @. B*Ju;
end

Bu = ABu(Js',Jr',BJu);

#Bu = Zygote.hook(d->hmp(d),Bu);
Bu = Zygote.hook(d->d .* mult,Bu);

Bu = gatherScatter(Bu,QQtx,QQty);
Bu = mask(Bu,M);

return Bu
end
#--------------------------------------#
function mass(f::Field,msh::Mesh)

    @unpack B   = msh
    @unpack u,M = f

    Bu = B .* u

    Bu = Field(Bu,M)
#   Bu = gatherScatter(Bu,QQtx,QQty)
#   Bu = mask(Bu,M)

return Bu
end
#--------------------------------------#


text:
######### Stan program example  ###########

using Stan, Mamba

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

const binorm = "
transformed data {
    matrix[2,2] Sigma;
    vector[2] mu;

    mu[1] <- 0.0;
    mu[2] <- 0.0;
    Sigma[1,1] <- 1.0;
    Sigma[2,2] <- 1.0;
    Sigma[1,2] <- 0.10;
    Sigma[2,1] <- 0.10;
}
parameters {
    vector[2] y;
}
model {
      y ~ multi_normal(mu,Sigma);
}
"

stanmodel = Stanmodel(name="binormal", model=binorm, Sample(save_warmup=true));

sim1 = stan(stanmodel, CmdStanDir=CMDSTAN_HOME)

## Subset Sampler Output
sim = sim1[1:size(sim1, 1), ["lp__", "y.1", "y.2"], 1:size(sim1, 3)]

describe(sim)
println()

## Brooks, Gelman and Rubin Convergence Diagnostic
try
  gelmandiag(sim1, mpsrf=true, transform=true) |> display
catch e
  #println(e)
  gelmandiag(sim, mpsrf=false, transform=true) |> display
end

## Geweke Convergence Diagnostic
gewekediag(sim) |> display

## Highest Posterior Density Intervals
hpd(sim) |> display

## Cross-Correlations
cor(sim) |> display

## Lag-Autocorrelations
autocor(sim) |> display

## Deviance Information Criterion
#dic(sim) |> display

## Plotting

p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-summaryplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-summaryplot-$(i).svg"`)
        end
      end : println()

end # cd


text:
export pointwise_log_likelihoods

const ARRAY_DIMS_WARNING = "The supplied array of mcmc samples indicates you have more
parameters than mcmc samples.This is possible, but highly unusual. Please check that your
array of mcmc samples has the following dimensions: [n_samples,n_parms,n_chains]."

"""
    pointwise_log_likelihoods(
        ll_fun::Function,
        samples::AbstractArray{<:Real,3},
        data;
        splat::Bool=true
    )

Compute the pointwise log likelihood.

# Arguments
  - $LIKELIHOOD_FUNCTION_ARG
  - `samples::AbstractArray`: A three dimensional array of MCMC samples. Here, the first
    dimension should indicate the iteration of the MCMC ; the second dimension should
    indicate the parameter ; and the third dimension represents the chains.
  - `data`: A vector of data used to estimate the parameters of the model.
  - `splat`: If `true` (default), `f` must be a function of `n` different parameters.
    Otherwise, `f` is assumed to be a function of a single parameter vector.

# Returns
  - `Array`: A three dimensional array of pointwise log-likelihoods.
"""
function pointwise_log_likelihoods(
    ll_fun::Function,
    samples::AbstractArray{<:Union{Real, Missing}, 3},
    data;
    splat::Bool=true,
)
    n_posterior, n_parms, n_chains = size(samples)
    if n_parms > n_posterior
        @info ARRAY_DIMS_WARNING
    end
    if splat
        fun = (p, d) -> ll_fun(p..., d)
    else
        fun = (p, d) -> ll_fun(p, d)
    end
    n_posterior, _, n_chains = size(samples)
    n_data = length(data)
    pointwise_lls = similar(samples, n_data, n_posterior, n_chains)
    for index in CartesianIndices(pointwise_lls)
        datum, iteration, chain = Tuple(index)
        pointwise_lls[datum, iteration, chain] = fun(
            samples[iteration, :, chain], data[datum]
        )
    end
    return pointwise_lls
end

function pointwise_log_likelihoods(
    ll_fun::Function,
    samples::AbstractMatrix{<:Union{Real, Missing}},
    data;
    chain_index::AbstractVector{<:Integer}=_assume_one_chain(samples),
    kwargs...,
)
    samples = _convert_to_array(samples, chain_index)
    return pointwise_log_likelihoods(ll_fun, samples, data)
end


text:
module test_0_ResultDict

using ModiaResult
using Unitful
ModiaResult.@usingModiaPlot

time0 = [0.0, 15.0]
t     = ([time0], [time0], ModiaResult.Independent)

time1 = 0.0  : 0.1 : 3.0
time2 = 5.0  : 0.1 : 9.5
time3 = 11.0 : 0.1 : 15.0
time4 = 0.0  : 0.1 : 15.0
sigA1 = 0.9*sin.(time1)u"m"
sigA2 = cos.(time2)u"m"
sigA3 = 1.1*sin.(time3)u"m"
sigA  = ([time1, time2, time3],
         [sigA1, sigA2, sigA3 ],
         ModiaResult.Continuous)
sigB  = ([time4], [0.7*sin.(time4)u"m/s"], ModiaResult.Continuous)
sigC  = ([time2], [sin.(time2)u"N*m"]    , ModiaResult.Clocked)

result = ModiaResult.ResultDict("time" => t,
                                "sigA" => sigA,
                                "sigB" => sigB,
                                "sigC" => sigC,
                                defaultHeading = "Signals from test_70_ResultDict.jl")

println("\n... test_70_ResultDict.jl:\n")
ModiaResult.printResultInfo(result)

plot(result, ("sigA", "sigB", "sigC"))

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
using Encounter
using Sensor
using CollisionAvoidanceSystem
using PilotResponse
using DynamicModel
using WorldModel

import Base.convert

#The fact that these are distributed between the simulators (i.e., TCASSimulatorImpl.jl),
#is a bit messy.
function convert(::Type{StochasticLinearPRCommand}, command::Union{CorrAEMCommand, LLAEMCommand})
    StochasticLinearPRCommand(command.t, command.v_d, command.h_d, command.psi_d, 0.0)
end

function convert(::Type{LLDetPRCommand}, command::Union{CorrAEMCommand, LLAEMCommand})
    LLDetPRCommand(command.t, command.v_d, command.h_d, command.psi_d)
end
function convert(::Type{LLDetPRAircraft}, state::ASWMState)
    LLDetPRAircraft(state.vh)
end

function convert(::Type{SimpleADMCommand}, command::StochasticLinearPRCommand)
    SimpleADMCommand(command.t, command.v_d, command.h_d, command.psi_d)
end

function convert(::Type{SimpleADMCommand}, command::LLDetPRCommand)
    SimpleADMCommand(command.t, command.v_d, command.h_d, command.psi_d)
end

function convert(::Type{ACASXSensorInput}, states::Vector{ASWMState})
    sr_states = ACASXSensorState[ ACASXSensorState(state.x,state.y,state.h,
        state.vx,state.vy,state.vh) for state in states]
  ACASXSensorInput(sr_states)
end

function convert(::Type{StochasticLinearPRRA},RA::ACASXOutput)
    ra_active = (RA.dh_min > -9999.0 || RA.dh_max < 9999.0)
    StochasticLinearPRRA(ra_active,RA.target_rate,RA.dh_min,RA.dh_max)
end

function convert(::Type{StochasticLinearPRRA},RA::Union{SimpleTCASResolutionAdvisory,Void})
    ra_active = RA != nothing
    StochasticLinearPRRA(ra_active,ra_active ? RA.h_d : 0.0,-9999.0,9999.0)
end

function convert(::Type{LLDetPRRA},RA::ACASXOutput)
    LLDetPRRA(RA.dh_min,RA.dh_max,RA.target_rate,RA.ddh)
end

function convert(::Type{LLDetPRRA},RA::Union{SimpleTCASResolutionAdvisory,Void})
    ra_active = RA != nothing
    if !ra_active
        #not active
        return LLDetPRRA(-9999.0,9999.0,0.0)  #TODO: put these hardcodings in a central location
    else
        #RA
        if RA.h_d > 0
            #climb
            return LLDetPRRA(RA.h_d, 9999.0, RA.h_d)
        elseif RA.h_d < 0
            #descend
            return LLDetPRRA(-9999.0, RA.h_d, RA.h_d)
        else
            #level
            return LLDetPRRA(RA.h_d, RA.h_d, RA.h_d)
        end
    end
    error("Shouldn't have gotten here")
end

function convert(::Type{LLADMCommand}, out::LLDetPROutput)
    LLADMCommand(out.t, out.h_d, deg2rad(out.psi_d), out.v_d, out.dh_min,
        out.dh_max, out.target_rate, out.ddh)
end

function convert(::Type{LLADMState}, s::CorrAEMInitialState)
    LLADMState(s.t, s.v, s.x, s.y, s.h, deg2rad(s.psi), asin(s.h_d / s.v), 0.0, 0.0)
end

convert(::Type{ASWMState}, s::LLADMOutputState) = ASWMState(s.t, s.x, s.y, s.h, s.vx, s.vy, s.vh)


text:
import SparseArrays.sparse

function lar2plasm(EV, FE)
    EV = sparse(EV)
    faces = Array{Array{Int64, 1}, 1}()

    for f in 1:size(FE, 1)

        verts = Array{Int64, 1}()
        done = Array{Int64, 1}()
        face = FE[f, :]
        edges = EV[face.nzind, :]

        push!(done, 1)
        vs = edges[1, :].nzind
        if face.nzval[1] < 0
            vs = vs[end:-1:1]
        end
        startv, nextv = vs

        # println(startv)
        push!(verts, startv)
        while nextv != startv
            e = setdiff(edges[:, nextv].nzind, done)[1]
            push!(done, e)
            vs = edges[e, :].nzind
            if face.nzval[e] < 0
                vs = vs[end:-1:1]
            end
            curv, nextv = vs
            push!(verts, curv)
        end

        push!(faces, verts)
    end

    "FV = " * string(faces)[15:end]
end

function lar2py(V, EV)
    text = "V = [[0,0]"
    for i in 1:size(V, 1)
        text = text * "," * string(vec(V[i, :]))
    end
    text = text * "]\n"

    text = text * "EV = ["
    for i in 1:size(EV, 1)
        if i > 1 text = text * "," end
        text = text * string(EV[i, :].nzind)
    end
    text = text * "]"
    return text
end

function visualize_numbers(set, font_size)

    f = open("tmp.py", "w")

    write(f, "from larlib import *\n" *
        lar2py(set[1], set[2]) * "\n" *
        lar2plasm(set[2], set[3]) * "\n" *
        "VV = AA(LIST)(range(len(V)))\n" *
        "submodel = STRUCT(MKPOLS((V,EV)))\n" *
        "VIEW(larModelNumbering(1,1,1)(V,[VV,EV,FV],submodel,$font_size))\n")

    close(f)
    run(`python tmp.py`)
end

function visualize_boundary(set, boundary, font_size)
    f = open("tmp.py", "w")

    write(f, "from larlib import *\n" *
        lar2py(set[1], set[2]) * "\n" *
        lar2plasm(set[2], set[3]) * "\n" *
        "VV = AA(LIST)(range(len(V)))\n" *
        "submodel = STRUCT(MKPOLS((V,EV)))\n" *
        "model = larModelNumbering(1,1,1)(V,[VV,EV,FV],submodel,$font_size)\n" *
        lar2py(set[1], set[2][boundary.nzind, :]) * "\n" *
        "overmodel = STRUCT(MKPOLS((V,EV)))\n" *
        "VIEW(STRUCT([model, COLOR(RED)(overmodel)]))")

    close(f)
    run(`python tmp.py`)
end


text:
let
    trajdata = get_test_trajdata()

    scene = get!(Scene(), trajdata, 1)
    @test length(scene) == 2
    for (i,veh) in enumerate(scene)
        @test scene[i].state == get_vehiclestate(trajdata, i, 1)
        @test scene[i].def == get_vehicledef(trajdata, i)
    end

    scene2 = Scene(deepcopy(scene.vehicles), 2)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_vehiclestate(trajdata, i, 1)
        @test scene2[i].def == get_vehicledef(trajdata, i)
    end

    @test get_by_id(scene, 1) == scene[1]

    empty!(scene2)
    @test length(scene2) == 0

    copy!(scene2, scene)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_vehiclestate(trajdata, i, 1)
        @test scene2[i].def == get_vehicledef(trajdata, i)
    end

    delete!(scene2, scene2[1])
    @test length(scene2) == 1
    @test scene2[1].state == get_vehiclestate(trajdata, 2, 1)
    @test scene2[1].def == get_vehicledef(trajdata, 2)
    scene2[1] = deepcopy(scene[1])
    @test scene2[1].state == get_vehiclestate(trajdata, 1, 1)
    @test scene2[1].def == get_vehicledef(trajdata, 1)

    @test get_index_of_first_vehicle_with_id(scene, 1) == 1
    @test get_index_of_first_vehicle_with_id(scene, 2) == 2
    @test get_first_available_id(scene) == 3

    @test iscarinframe(scene, 1)
    @test iscarinframe(scene, 2)
    @test !iscarinframe(scene, 3)

    veh = get_vehicle(scene, 2)
    @test veh.state == get_vehiclestate(trajdata, 2, 1)
    @test veh.def == get_vehicledef(trajdata, 2)

    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 1, trajdata.roadway) == NeighborLongitudinalResult(2, 3.0078125)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 2, trajdata.roadway) == NeighborLongitudinalResult(0, 250.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 1, trajdata.roadway) == NeighborLongitudinalResult(2, 4.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 2, trajdata.roadway) == NeighborLongitudinalResult(0, 250.0)

    push!(scene, get_vehiclestate(trajdata, 1, 1))
end

let
    roadway = gen_stadium_roadway(1)
    scene = Scene(2)

    function place_at(s)
        roadproj = proj(VecSE2(0.0,0.0,0.0), roadway)
        roadind = RoadIndex(roadproj.curveproj.ind, roadproj.tag)
        roadind = move_along(roadind, roadway, s)
        frenet = Frenet(roadind, roadway[roadind].s, 0.0, 0.0)
        VehicleState(frenet, roadway, 0.0)
    end

    push!(scene, Vehicle(place_at(0.0), VehicleDef(1, AgentClass.CAR, 2.0, 1.0)))
    push!(scene, Vehicle(place_at(0.0), VehicleDef(2, AgentClass.CAR, 2.0, 1.0)))
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test isapprox(foreinfo.Δs, 0.0)

    scene[2].state = place_at(100.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 100.0)

    foreinfo = get_neighbor_fore_along_lane(scene, 2, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 1
    @test isapprox(foreinfo.Δs, 277.07, atol=1e-2)

    scene[2].state = place_at(145.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 145.0, atol=1e-5)

    scene[2].state = place_at(240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 240.0, atol=1e-5)

    scene[1].state = place_at(240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 0.0, atol=1e-5)

    ################################################
    # get_frenet_relative_position

    scene[1].state = place_at(0.0)
    scene[2].state = place_at(10.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 0
    @test isapprox(frp.target.ind.t, 0.1, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    scene[1].state = place_at(10.0)
    scene[2].state = place_at(20.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 0
    @test isapprox(frp.target.ind.t, 0.2, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    scene[1].state = place_at(0.0)
    scene[2].state = place_at(120.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 120.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    scene[1].state = place_at(0.0)
    scene[2].state = place_at(250.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag != scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 250.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)
end


text:
include("../examples/lorenz63.jl")
include("../src/lss.jl")
using LinearAlgebra
function refine_orbit(n)
	#create epsilon orbit.
	n_runup = 5000
    s = [10., 28., 8/3]
    u0 = rand(3,1)
	u0 = lorenz63(u0, s, n_runup)[end,:,:]
	u_trj = lorenz63(u0, s, n)[:,:,1]
	n, d = size(u_trj)
	eps = 0.1
	noise = eps*randn(n)
	u_obs = copy(u_trj)
	u_obs[:,3] .+= noise
	# run gradient descent to minimize observation error.
	n_gd_steps = 4
	gamma = 0.1
	dJds = zeros(n_gd_steps)
	d_u = 2
	u = zeros(n, d, n_gd_steps)
	# compute_sens assumes many objective functions
	dJ = zeros(1,d,n)
	error = zeros(n_gd_steps)
	for i = 1:n_gd_steps
		# Set up LSS
		# J is now sum of squared observation error.
		dJ[1,3,:] = -2.0*noise
    	du_trj = dlorenz63(u_trj, s)
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	#f = vectorField(u_trj,s)
    	f = zeros(d,n)
    	y, xi = lss(du_trj, X, f, s, d_u)
		dJds[i] = compute_sens(y, xi, dJ, f)[1]
    	println(dJds[i])
		# Make u_trj the computed shadowing trajectory
		ds = gamma*dJds[i]
		s[2] = s[2] - ds
		u_trj = u_trj + ds*y'
		error[i] =  sum((u_trj[:,3] .- u_obs[:,3]).^2)/n
		u[:,:,i] = u_trj
	end
	return u_obs[:,3], u[:,3,:], error
end
function assimilate()
	n_repeat = 1000
	n = 2000
	n_gd_steps = 4
	z_obs = zeros(n, n_repeat)
	z_trj = zeros(n, n_gd_steps, n_repeat)
	errors = zeros(n_gd_steps, n_repeat)
	for i = 1:n_repeat
		res1, res2, res3 = refine_parameter(n-1)
		z_obs[:,i] = res1
		z_trj[:,:,i] = res2
		errors[:,i] = res3
	end
	return z_obs, z_trj, errors
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
export fparse, fstring

function fparse(::Type{T}, s::AbstractString) where {T<:Integer}
    return parse(T, s)
end
function fparse(::Type{T}, s::AbstractString) where {T<:Float32}
    return parse(T, replace(lowercase(s), r"(?<=[^e])(?=[+-])" => "f"))
end
function fparse(::Type{T}, s::AbstractString) where {T<:Float64}
    return parse(T, replace(lowercase(s), r"d"i => "e"))
end
function fparse(::Type{Complex{T}}, str::AbstractString) where {T<:AbstractFloat}
    if first(str) == '(' && last(str) == ')' && length(split(str, ',')) == 2
        re, im = split(str[2:end-1], ',', limit = 2)
        return Complex(parse(T, re), parse(T, im))
    else
        throw(Meta.ParseError("$str must be in complex number form (x, y)."))
    end
end
function fparse(::Type{Bool}, s::AbstractString)
    str = lowercase(s)
    if str in (".true.", ".t.", "true", 't')
        return true
    elseif str in (".false.", ".f.", "false", 'f')
        return false
    else
        throw(Meta.ParseError("$str is not a valid logical constant."))
    end
end
function fparse(::Type{String}, str::AbstractString)
    m = match(r"([\"'])((?:\\\1|.)*?)\1", str)
    if m === nothing
        throw(Meta.ParseError("$str is not a valid string!"))
    else
        quotation_mark, content = m.captures
        # Replace escaped strings
        return string(replace(content, repeat(quotation_mark, 2) => quotation_mark))
    end
end

fstring(v::Integer) = string(v)
function fstring(v::Float32; scientific::Bool = false)
    str = string(v)
    return scientific ? replace(str, r"f"i => "e") : str
end
function fstring(v::Float64, scientific::Bool = false)
    str = string(v)
    return scientific ? replace(str, r"e"i => "d") : str
end
fstring(v::Bool) = v ? ".true." : ".false."
fstring(v::Union{AbstractString,AbstractChar}) = "'" * v * "'"
# fstring(::Namelist)  # TODO:


text:
# This file includes stepping of abstract types.

##### Input-Output reading and writing.
"""
    $SIGNATURES

Returns current time of `comp` read from its `trigger` link.

!!! note
    To read time of `comp`, `comp` must be launched. See also: [`launch(comp::AbstractComponent)`](@ref).
"""
readtime!(comp::AbstractComponent) = take!(comp.trigger)

"""
    $SIGNATURES

Returns the state of `comp` if `comp` is `AbstractDynamicSystem`. Otherwise, returns `nothing`.
"""
readstate(comp::AbstractComponent) = typeof(comp) <: AbstractDynamicSystem ? comp.state : nothing

"""
    $SIGNATURES

Returns the input value of `comp` if the `input` of `comp` is `Inport`. Otherwise, returns `nothing`.

!!! note
    To read input value of `comp`, `comp` must be launched. See also: [`launch(comp::AbstractComponent)`](@ref)
"""
function readinput!(comp::AbstractComponent)
    typeof(comp) <: AbstractSource && return nothing
    typeof(comp.input) <: Inport ? take!(comp.input) : nothing
end

"""
    $SIGNATURES

Writes `out` to the output of `comp` if the `output` of `comp` is `Outport`. Otherwise, does `nothing`.
"""
function writeoutput!(comp::AbstractComponent, out)
    typeof(comp) <: AbstractSink && return nothing
    typeof(comp.output) <: Outport ? put!(comp.output, out) : nothing
end

"""
    $SIGNATURES

Computes the output of `comp` according to its `readout` if `readout` is not `nothing`. Otherwise, `nothing` is done. `x` is
the state, `u` is the value of input, `t` is the time.
"""
function computeoutput end
computeoutput(comp::AbstractSource, x, u, t) = comp.readout(t)
computeoutput(comp::AbstractStaticSystem, x, u, t) =
    typeof(comp.readout) <: Nothing ? nothing : comp.readout(u, t)
function computeoutput(comp::AbstractDynamicSystem, x, u, t)
    typeof(comp.readout) <: Nothing && return nothing
    typeof(u) <: Nothing ? comp.readout(x, u, t) : comp.readout(x, map(uu -> t -> uu, u), t)
end
    # typeof(comp.readout) <: Nothing ? nothing : comp.readout(x, constructinput(comp, u, t), t)
computeoutput(comp::AbstractSink, x, u, t) = nothing

"""
    evolve!(comp::AbstractSource, u, t)

Does nothing. `u` is the value of `input` and `t` is time.

    evolve!(comp::AbstractSink, u, t)

Writes `t` to time buffer `timebuf` and `u` to `databuf` of `comp`. `u` is the value of `input` and `t` is time.

    evolve!(comp::AbstractStaticSystem, u, t)

Writes `u` to `buffer` of `comp` if `comp` is an `AbstractMemory`. Otherwise, `nothing` is done. `u` is the value of `input`
and `t` is time.

    evolve!(comp::AbstractDynamicSystem, u, t)

Solves the differential equation of the system of `comp` for the time interval `(comp.t, t)` for the inital condition `x`
where `x` is the current state of `comp` . `u` is the input function defined for `(comp.t, t)`. The `comp` is updated with
the computed state and time `t`.
"""
function evolve! end
evolve!(comp::AbstractSource, u, t) = nothing
evolve!(comp::AbstractSink, u, t) = (write!(comp.timebuf, t); write!(comp.databuf, u); comp.sinkcallback(comp); nothing)
function evolve!(comp::AbstractStaticSystem, u, t)
    if typeof(comp) <: AbstractMemory
        timebuf = comp.timebuf
        databuf = comp.databuf
        write!(timebuf, t)
        write!(databuf, u)
    end
end
function evolve!(comp::AbstractDynamicSystem, u, t)
    # For DDESystems, the problem for a time span of (t, t) cannot be solved. Thus, there will be no evolution in such a
    # case.
    integrator = comp.integrator
    interpolator = integrator.sol.prob.p
    update_interpolator!(interpolator, u, t)
    comp.t == t && return comp.state

    # Advance the system and update the system.
    step!(integrator, t - comp.t, true)
    comp.t = integrator.t
    comp.state = integrator.u

    # Return comp state
    return comp.state
end
update_interpolator!(interp::Nothing) = nothing
update_interpolator!(interp::Nothing, u, t) = nothing
function update_interpolator!(interp::Interpolant, u, t)
    write!(interp.timebuf, t)
    write!(interp.databuf, u)
    update!(interp)
end

##### Task management
"""
    $SIGNATURES

Reads the time `t` from the `trigger` link of `comp`. If `comp` is an `AbstractMemory`, a backward step is taken. Otherwise,
a forward step is taken. See also: [`forwardstep`](@ref), [`backwardstep`](@ref).
"""
function takestep!(comp::AbstractComponent)
    t = readtime!(comp)
    t === NaN && return t
    typeof(comp) <: AbstractMemory ? backwardstep(comp, t) : forwardstep(comp, t)
end

"""
    $SIGNATURES

Makes `comp` takes a forward step.  The input value `u` and state `x` of `comp` are read. Using `x`, `u` and time `t`,
`comp` is evolved. The output `y` of `comp` is computed and written into the output bus of `comp`.
"""
function forwardstep(comp, t)
    u = readinput!(comp)
    x = evolve!(comp, u, t)
    y = computeoutput(comp, x, u, t)
    writeoutput!(comp, y)
    applycallbacks(comp)
    return t
end

"""
    $SIGNATURES

Reads the state `x`. Using the time `t` and `x`, computes and writes the ouput value `y` of `comp`. Then, the input value `u`
is read and `comp` is evolved.
"""
function backwardstep(comp, t)
    x = readstate(comp)
    y = computeoutput(comp, x, nothing, t)
    writeoutput!(comp, y)
    u = readinput!(comp)
    xn = evolve!(comp, u, t)
    applycallbacks(comp)
    return t
end

"""
    $SIGNATURES

Returns a tuple of tasks so that `trigger` link and `output` bus of `comp` is drivable. When launched, `comp` is ready to be
driven from its `trigger` link. See also: [`drive!(comp::AbstractComponent, t)`](@ref)
"""
function launch(comp::AbstractComponent)
    @async begin
        while true
            takestep!(comp) === NaN && break
            put!(comp.handshake, true)
        end
        typeof(comp) <: AbstractSink && close(comp)
    end
end

"""
    $SIGNATURES

Writes `t` to the `trigger` link of `comp`. When driven, `comp` takes a step. See also:
[`takestep!(comp::AbstractComponent)`](@ref)
"""
drive!(comp::AbstractComponent, t) = put!(comp.trigger, t)

"""
    $SIGNATURES

Read `handshake` link of `comp`. When not approved or `false` is read from the `handshake` link, the task launched for the
`trigger` link of `comp` gets stuck during `comp` is taking step.
"""
approve!(comp::AbstractComponent) = take!(comp.handshake)

"""
    $SIGNATURES

Closes the `trigger` link and `output` bus of `comp`.
"""
function terminate!(comp::AbstractComponent)
    typeof(comp) <: AbstractSink || typeof(comp.output) <: Nothing || close(comp.output)
    close(comp.trigger)
    return
end

##### SubSystem interface
"""
    $SIGNATURES

Launches all subcomponents of `comp`. See also: [`launch(comp::AbstractComponent)`](@ref)
"""
function launch(comp::AbstractSubSystem)
    comptask = @async begin
        while true
            if takestep!(comp) === NaN
                put!(comp.triggerport, fill(NaN, length(comp.components)))
                break
            end
            put!(comp.handshake, true)
        end
    end
    [launch.(comp.components)..., comptask]
end

"""
    $SIGNATURES

Makes `comp` to take a step by making each subcomponent of `comp` take a step. See also:
[`takestep!(comp::AbstractComponent)`](@ref)
"""
function takestep!(comp::AbstractSubSystem)
    t = readtime!(comp)
    t === NaN && return t
    put!(comp.triggerport, fill(t, length(comp.components)))
    all(take!(comp.handshakeport)) || @warn "Could not be approved in the subsystem"
    # foreach(takestep!, comp.components) approve!(comp) ||  @warn "Could not be approved in the subsystem"
    # put!(comp.handshake, true)
end

"""
    $SIGNATURES

Drives `comp` by driving each subcomponent of `comp`. See also: [`drive!(comp::AbstractComponent, t)`](@ref)
"""
drive!(comp::AbstractSubSystem, t) = foreach(component -> drive!(component, t), comp.components)

"""
    $SIGNATURES

Approves `comp` by approving each subcomponent of `comp`. See also: [`approve!(comp::AbstractComponent)`](@ref)
"""
approve!(comp::AbstractSubSystem) = all(approve!.(comp.components))

"""
    $SIGNATURES

Terminates `comp` by terminating each subcomponent of `comp`. See also: [`terminate!(comp::AbstractComponent)`](@ref)
"""
terminate!(comp::AbstractSubSystem) = foreach(terminate!, comp.components)


text:
using Turing, ContinuousBenchmarks
using JLD2, FileIO, Profile, ProfileView

setadbackend(:reverse_diff)
turnprogress(false)

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "lda-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda-stan.model.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda.model.jl"))

data = map(ind -> ldastandata[1][ind], ["K", "V", "M", "N", "w", "doc", "beta", "alpha"])
sample(ldamodel(data...), HMC(2, 0.025, 10))
Profile.clear()
@profile sample(ldamodel(data...), HMC(20, 0.025, 10))

ProfileView.svgwrite("ldamodel.svg")


text:
mutable struct RI
    water
    cornea
    lens
    lens2distal_retina
    distal_retina2proximal_retina
    gap
    behind_mirror
end

mutable struct Thick
    water
    cornea
    lens
    lens2distal_retina
    distal_retina2proximal_retina
    gap
    behind_mirror
end

mutable struct Radii
    cornea_distal
    lens_distal
    lens_proximal
    distal_retina
    proximal_retina
    mirror
end

ri = RI(zeros(7)...)
thick = Thick(μm*zeros(7)...)
radii = Radii(μm*zeros(6)...)


text:
using Images
using ImageSegmentation
using BenchmarkTools
using ProgressMeter

run_complement(img) = complement.(img)

setup_mean(img) = repeat(img, ntuple(_->1, ndims(img))..., 5)
function run_mean(imgr)
    return mean(imgr; dims=ndims(imgr))
end

run_gradient(img) = imgradients(img, KernelFactors.sobel)

run_blur(img) = imfilter(img, KernelFactors.gaussian(ntuple(_->5.0f0, ndims(img)), ntuple(_->21, ndims(img))))
run_blur_iir(img) = imfilter(img, KernelFactors.IIRGaussian(ntuple(_->5, ndims(img)))) # a bit faster

run_histeq(img) = adjust_histogram(img, Equalization(nbins = 256, minval = 0, maxval = 1))

function time_generics(workdir)
    tdata = Dict(task=>Dict{String,Float64}() for task in ["complement", "mean", "gradient", "blur", "histeq"])
    time_generics!(tdata, workdir, keys(tdata))
end

function time_generics!(tdata, workdir, taskitems)
    @showprogress "Timing generic operations: " for fn in readdir(workdir)
        img = load(joinpath(workdir, fn))
        fn, _ = splitext(fn)
        "complement" ∈ taskitems && (tdata["complement"][fn] = @belapsed run_complement($img) samples=16 evals=1)
        "mean" ∈ taskitems       && (tdata["mean"][fn] = @belapsed run_mean(imgr) samples=16 evals=1 setup=(imgr = setup_mean($img)))
        "gradient" ∈ taskitems   && (tdata["gradient"][fn] = @belapsed run_gradient($img) samples=16 evals=1)
        "blur" ∈ taskitems       && (tdata["blur"][fn] = @belapsed run_blur($img) samples=16 evals=1)
        "histeq" ∈ taskitems     && (tdata["histeq"][fn] = @belapsed run_histeq($img) samples=16 evals=1)
    end
    return tdata
end

# tdata = time_generics("/tmp/imgs")

run_label(img) = label_components(img, trues(ntuple(_->3, ndims(img))))

run_disttform(img) = distance_transform(feature_transform(img))

function run_flood(img)
    middle = map(axes(img)) do ax
        (first(ax) + last(ax)) ÷ 2
    end
    flood(>=(0.5), img, CartesianIndex(middle))
end

function time_special(workdir)
    imgs = ["components2d" => run_label, "components3d" => run_label,
            "dblcone2d" => run_disttform, "dblcone3d" => run_disttform,
            "spiral2d" => run_flood, "spiral3d" => run_flood]
    fls = readdir(workdir)
    tdata = Dict{String,Float64}()
    @showprogress "Timing special operations: " for (name, f) in imgs
        idx = findfirst(str->startswith(str, name), fls)
        img = load(joinpath(workdir, fls[idx]))
        if f === run_disttform
            img = convert(Array{Bool}, img)
        end
        tdata[name] = @belapsed $f($img) evals=1
    end
    return tdata
end

# tdata = time_special("/tmp/imgspecial")


text:
using Test

@testset "adt" begin
include("adt.jl")
end


text:
using DocStringExtensions
using Base.Test

VERSION < v"0.5.0-dev" ? warn("Untestable on Julia 0.4.") : include("tests.jl")