text:
#
# MikrotronMC408x.jl --
#
# Configuration and control of CoaXPress camera Mikrotron MC408x.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

module MikrotronMC408x

using ScientificCameras, ScientificCameras.PixelFormats
for sym in names(ScientificCameras)
    if sym != :ScientificCameras
        @eval begin
            import ScientificCameras: $sym
        end
    end
end
import ScientificCameras: ScientificCamera, ROI, setspeed!

using Phoenix
for sym in Phoenix._DEVEL_SYMBOLS
    @eval begin
        import Phoenix: $sym
    end
end

macro exportmethods()
    :(export
      setfixedpatternnoisereduction!,
      getfixedpatternnoisereduction,
      setinfofieldframecounter!,
      getinfofieldframecounter,
      setinfofieldtimestamp!,
      getinfofieldtimestamp,
      setinfofieldroi!,
      getinfofieldroi,
      setfiltermode!,
      getfiltermode)
end
@exportmethods

#
# CoaXPress camera constants for a Mikrotron MC408x camera.  These values have
# been extracted from the XML configuration file of the camera.
#
# FIXME: The XML file says that register HsUpConnection is at address 0x4030,
#         while CoaXPress standard says that it is at address 0x403C.
#

const ACQUISITION_MODE  = RegisterValue{UInt32,ReadWrite}(0x8200)
const ACQUISITION_START = RegisterCommand{UInt32}(0x8204, 1)
const ACQUISITION_STOP  = RegisterCommand{UInt32}(0x8208, 1)

const ACQUISITION_BURST_FRAME_COUNT = RegisterValue{UInt32,ReadWrite}(0x8914)
const TRIGGER_SELECTOR              = RegisterEnum{ReadWrite}(0x8900)
const TRIGGER_MODE                  = RegisterEnum{ReadWrite}(0x8904)
const TRIGGER_SOURCE                = RegisterEnum{ReadWrite}(0x8908)
const TRIGGER_ACTIVATION            = RegisterEnum{ReadWrite}(0x890C)
const TRIGGER_COUNT                 = RegisterValue{Nothing,Unreachable}(0x891C)
const TRIGGER_DEBOUNCER             = RegisterValue{Nothing,Unreachable}(0x8918)
const TRIGGER_SOFTWARE              = RegisterEnum{WriteOnly}(0x8910)
const TEST_IMAGE_SELECTOR           = RegisterEnum{ReadWrite}(0x9000)
const EXPOSURE_MODE                 = RegisterEnum{ReadWrite}(0x8944)
const EXPOSURE_TIME                 = RegisterValue{UInt32,ReadWrite}(0x8840)
const EXPOSURE_TIME_MAX             = RegisterValue{UInt32,ReadOnly}(0x8818)
const ACQUISITION_FRAME_RATE        = RegisterValue{UInt32,ReadWrite}(0x8814)
const ACQUISITION_FRAME_RATE_MIN    = 10 # FIXME: the XML file says 16
const ACQUISITION_FRAME_RATE_MAX    = RegisterValue{UInt32,ReadOnly}(0x881C)
const SEQUENCER_CONFIGURATION_MODE  = RegisterValue{UInt32,Unreachable}(0x8874)
const SEQUENCER_MODE                = RegisterValue{UInt32,Unreachable}(0x8870)
const SEQUENCER_SET_SELECTOR        = RegisterValue{Nothing,Unreachable}(0x8878)
const SEQUENCER_SET_SAVE            = RegisterValue{Nothing,Unreachable}(0x887C)
const SEQUENCER_SET_NEXT            = RegisterValue{Nothing,Unreachable}(0x8888)
const USER_SET_SELECTOR             = RegisterEnum{ReadWrite}(0x8820)
const USER_SET_LOAD                 = RegisterCommand{UInt32}(0x8824, 1)
const USER_SET_SAVE                 = RegisterCommand{UInt32}(0x8828, 1)
const USER_SET_DEFAULT_SELECTOR     = RegisterEnum{ReadWrite}(0x882C)
const DEVICE_RESET                  = RegisterCommand{UInt32}(0x8300, 1)
const CONNECTION_RESET              = RegisterCommand{UInt32}(0x4000, 1)
const CONNECTION_CONFIG             = RegisterValue{UInt32,ReadWrite}(0x4014)
const CONNECTION_DEFAULT            = RegisterValue{UInt32,ReadOnly}(0x4018)

# Image Format Control.
const HORIZONTAL_INCREMENT          = 16
const VERTICAL_INCREMENT            =  2
const REGION_SELECTOR               = RegisterEnum{ReadWrite}(0x8180)
const REGION_MODE                   = RegisterEnum{ReadWrite}(0x8184)
const REGION_DESTINATION            = RegisterEnum{ReadWrite}(0x8188) # always 0x0, i.e. "Stream1"
const WIDTH                         = RegisterValue{UInt32,ReadWrite}(0x8118)
const HEIGHT                        = RegisterValue{UInt32,ReadWrite}(0x811c)
const OFFSET_X                      = RegisterValue{UInt32,ReadWrite}(0x8800)
const OFFSET_Y                      = RegisterValue{UInt32,ReadWrite}(0x8804)
const DECIMATION_HORIZONTAL         = RegisterValue{UInt32,ReadWrite}(0x8190)
const DECIMATION_VERTICAL           = RegisterValue{UInt32,ReadWrite}(0x818c)
const SENSOR_WIDTH                  = RegisterValue{UInt32,ReadOnly}(0x8808)
const SENSOR_HEIGHT                 = RegisterValue{UInt32,ReadOnly}(0x880c)
# WidthMax = ((div(SENSOR_WIDTH - OFFSET_X), 16)*16) # use a bit mask
# HeightMax = ((div(SENSOR_HEIGHT - OFFSET_Y), 2)*2) # use a bit mask
const TAP_GEOMETRY = RegisterEnum{ReadOnly}(0x8160) # always 0x0, i.e. "Geometry_1X_1Y"
const PIXEL_FORMAT = RegisterEnum{ReadWrite}(0x8144)
const IMAGE1_STREAM_ID = RegisterValue{UInt32,ReadOnly}(0x8164)
# const DEVICE_SCANTYPE = RegisterEnum{ReadOnly}(????) always 0x0, i.e. "Areascan"

# Possible values for CXP_DEVICE_INFORMATION_SELECTOR.
const DEVICE_INFORMATION_SELECTOR_SERIAL_NUMBER     =  0
const DEVICE_INFORMATION_SELECTOR_DEVICE_TYPE       =  1
const DEVICE_INFORMATION_SELECTOR_DEVICE_SUBTYPE    =  2
const DEVICE_INFORMATION_SELECTOR_HARDWARE_REVISION =  3
const DEVICE_INFORMATION_SELECTOR_FPGA_VERSION      =  4
const DEVICE_INFORMATION_SELECTOR_SOFTWARE_VERSION  =  5
const DEVICE_INFORMATION_SELECTOR_POWER_SOURCE      = 20
const DEVICE_INFORMATION_SELECTOR_POWER_CONSUMPTION = 21
const DEVICE_INFORMATION_SELECTOR_POWER_VOLTAGE     = 22
const DEVICE_INFORMATION_SELECTOR_TEMPERATURE       = 23

const REGION_MODE_OFF = UInt32(0x00)
const REGION_MODE_ON  = UInt32(0x01)

const PIXEL_FORMAT_MONO8     = UInt32(0x0101)
const PIXEL_FORMAT_MONO10    = UInt32(0x0102)
const PIXEL_FORMAT_BAYERGR8  = UInt32(0x0311)
const PIXEL_FORMAT_BAYERGR10 = UInt32(0x0312)

const REGION_SELECTOR_REGION0 = UInt32(0x0)
const REGION_SELECTOR_REGION1 = UInt32(0x1)
const REGION_SELECTOR_REGION2 = UInt32(0x2)

# Bits for CONNECTION_CONFIG register (the value is a conbination of speed and
# number of connections):
const CONNECTION_CONFIG_SPEED1250 = UInt32(0x00028)
const CONNECTION_CONFIG_SPEED2500 = UInt32(0x00030)
const CONNECTION_CONFIG_SPEED3125 = UInt32(0x00038)
const CONNECTION_CONFIG_SPEED5000 = UInt32(0x00040)
const CONNECTION_CONFIG_SPEED6250 = UInt32(0x00048)
const CONNECTION_CONFIG_CONNECTION1 = UInt32(0x10000)
const CONNECTION_CONFIG_CONNECTION2 = UInt32(0x20000)
const CONNECTION_CONFIG_CONNECTION3 = UInt32(0x30000)
const CONNECTION_CONFIG_CONNECTION4 = UInt32(0x40000)

const GAIN                            = RegisterValue{UInt32,ReadWrite}(0x8850)
const GAIN_MIN                        =   50
const GAIN_MAX                        = 1000
const BLACK_LEVEL                     = RegisterValue{UInt32,ReadWrite}(0x8854)
const BLACK_LEVEL_MIN                 =    0
const BLACK_LEVEL_MAX                 =  500
const GAMMA                           = RegisterValue{Float32,ReadWrite}(0x8858)
const GAMMA_MIN                       =  0.1
const GAMMA_MAX                       =  3.0
const GAMMA_INCREMENT                 =  0.1
#const LINE_SOURCE                    = RegisterEnum{ReadWrite}(0x????)
#const LINE_SELECTOR                  = RegisterEnum{ReadWrite}(0x????)
const LINE_INVERTER                   = RegisterEnum{ReadWrite}(0x8A20)
const TX_LOGICAL_CONNECTION_RESET     = RegisterValue{Nothing,Unreachable}(0x9010)
const PRST_ENABLE                     = RegisterValue{UInt32,ReadWrite}(0x9200)
const PULSE_DRAIN_ENABLE              = RegisterValue{UInt32,ReadWrite}(0x9204)
const CUSTOM_SENSOR_CLK_ENABLE        = RegisterValue{Nothing,Unreachable}(0x9300)
const CUSTOM_SENSOR_CLK               = RegisterValue{Nothing,Unreachable}(0x9304)
const DEVICE_INFORMATION              = RegisterValue{UInt32,ReadOnly}(0x8A04)
const DEVICE_INFORMATION_SELECTOR     = RegisterValue{UInt32,ReadWrite}(0x8A00)
const ANALOG_REGISTER_SET_SELECTOR    = RegisterValue{Nothing,Unreachable}(0x20000)
const ANALOG_REGISTER_SELECTOR        = RegisterValue{Nothing,Unreachable}(0x20004)
const ANALOG_VALUE                    = RegisterValue{Nothing,Unreachable}(0x20008)
const INFO_FIELD_FRAME_COUNTER_ENABLE = RegisterValue{UInt32,ReadWrite}(0x9310)
const INFO_FIELD_TIME_STAMP_ENABLE    = RegisterValue{UInt32,ReadWrite}(0x9314)
const INFO_FIELD_ROI_ENABLE           = RegisterValue{UInt32,ReadWrite}(0x9318)
const FIXED_PATTERN_NOISE_REDUCTION   = RegisterValue{UInt32,ReadWrite}(0x8A10)
const FILTER_MODE                     = RegisterEnum{ReadWrite}(0x10014)
const PIXEL_TYPE_F                    = RegisterValue{Nothing,Unreachable}(0x51004)
const DIN1_CONNECTOR_TYPE             = RegisterValue{Nothing,Unreachable}(0x8A30)
const IS_IMPLEMENTED_MULTI_ROI        = RegisterValue{Nothing,Unreachable}(0x50004)
const IS_IMPLEMENTED_SEQUENCER        = RegisterValue{Nothing,Unreachable}(0x50008)
const CAMERA_TYPE_HEX                 = RegisterValue{Nothing,Unreachable}(0x51000)
const CAMERA_STATUS                   = RegisterValue{Nothing,Unreachable}(0x10002200)
const IS_STOPPED                      = RegisterValue{Nothing,Unreachable}(0x10002204)

const FILTER_MODE_RAW                 = UInt32(0)
const FILTER_MODE_MONO                = UInt32(1)
const FILTER_MODE_COLOR               = UInt32(2)

# Singleton to uniquely identify this camera model.
struct MikrotronMC408xModel <: CameraModel; end

# Initialize the camera after board is open.
function _openhook(cam::Camera{MikrotronMC408xModel})

    # Sanity checks.
    assert_coaxpress(cam)
    vendorname = cam[CXP_DEVICE_VENDOR_NAME]
    if vendorname != "Mikrotron GmbH"
        error("bad device vendor name (got \"$vendorname\", expecting \"Mikrotron GmbH\")")
    end
    modelname = cam[CXP_DEVICE_MODEL_NAME]
    if length(modelname) != 6 || modelname[1:5] != "MC408" ||
        (modelname[6] != '2' && modelname[6] != '3' &&
         modelname[6] != '6' && modelname[6] != '7')
        error("bad device model name (got \"$modelname\", expecting \"MC408[2367]\")")
    end

    # Get current hardware settings.
    xsub   = Int(cam[DECIMATION_HORIZONTAL]) # in pixels
    ysub   = Int(cam[DECIMATION_VERTICAL])   # in pixels
    width  = Int(cam[WIDTH])                 # in macro-pixels
    height = Int(cam[HEIGHT])                # in macro-pixels

    # The following settings are the same as the contents of the configuration
    # file "Mikrotron_MC4080_CXP.pcf".
    cam[PHX_BOARD_VARIANT]      = PHX_DIGITAL
    cam[PHX_CAM_TYPE]           = PHX_CAM_AREASCAN_ROI
    cam[PHX_CAM_FORMAT]         = PHX_CAM_NON_INTERLACED
    cam[PHX_CAM_CLOCK_POLARITY] = PHX_CAM_CLOCK_POS
    cam[PHX_CAM_SRC_COL]        = PHX_CAM_SRC_MONO
    cam[PHX_CAM_DATA_VALID]     = PHX_DISABLE
    cam[PHX_CAM_HTAP_NUM]       = 1
    cam[PHX_CAM_HTAP_DIR]       = PHX_CAM_HTAP_LEFT
    cam[PHX_CAM_HTAP_TYPE]      = PHX_CAM_HTAP_LINEAR
    cam[PHX_CAM_HTAP_ORDER]     = PHX_CAM_HTAP_ASCENDING
    cam[PHX_CAM_VTAP_NUM]       = 1
    cam[PHX_CAM_VTAP_DIR]       = PHX_CAM_VTAP_TOP
    cam[PHX_CAM_VTAP_TYPE]      = PHX_CAM_VTAP_LINEAR
    cam[PHX_CAM_VTAP_ORDER]     = PHX_CAM_VTAP_ASCENDING
    cam[PHX_COMMS_DATA]         = PHX_COMMS_DATA_8
    cam[PHX_COMMS_STOP]         = PHX_COMMS_STOP_1
    cam[PHX_COMMS_PARITY]       = PHX_COMMS_PARITY_NONE
    cam[PHX_COMMS_SPEED]        = 9600
    cam[PHX_COMMS_FLOW]         = PHX_COMMS_FLOW_NONE

    # Set the format of the image sent by the camera.
    srcformat, srcdepth, dstformat = equivalentformat(cam)
    cam[PHX_CAM_SRC_COL]   = srcformat
    cam[PHX_CAM_SRC_DEPTH] = srcdepth
    cam[PHX_DST_FORMAT]    = dstformat
    setactiveregion!(cam, width, height)

    # Acquisition settings.
    cam[PHX_ACQ_BLOCKING]          = PHX_ENABLE
    cam[PHX_ACQ_CONTINUOUS]        = PHX_ENABLE
    cam[PHX_ACQ_XSUB]              = PHX_ACQ_X1
    cam[PHX_ACQ_YSUB]              = PHX_ACQ_X1
    cam[PHX_ACQ_NUM_BUFFERS]       = 1
    cam[PHX_ACQ_IMAGES_PER_BUFFER] = 1
    cam[PHX_ACQ_BUFFER_START]      = 1
    cam[PHX_DATASTREAM_VALID]      = PHX_DATASTREAM_ALWAYS
    cam[PHX_TIMEOUT_DMA]           = 1_000 # milliseconds

    # Set source ROI to match the size of the image sent by the camera.
    setsourceregion!(cam, 0, 0, width, height)

    # Setup destination buffer parameters.  The value of `PHX_BUF_DST_XLENGTH`
    # is the number of bytes per line of the destination buffer (it must be
    # larger of equal the width of the ROI times the number of bits per pixel
    # rounded up to a number of bytes), the value of `PHX_BUF_DST_YLENGTH` is
    # the number of lines in the destination buffer (it must be larger or equal
    # `PHX_ROI_DST_YOFFSET` plus `PHX_ROI_YLENGTH`.
    bits = capture_format_bits(dstformat)
    cam[PHX_ROI_DST_XOFFSET] = 0
    cam[PHX_ROI_DST_YOFFSET] = 0
    cam[PHX_BUF_DST_XLENGTH] = div(width*bits + 7, 8)
    cam[PHX_BUF_DST_YLENGTH] = height
    cam[PHX_BIT_SHIFT]       = 0 # FIXME: PHX_BIT_SHIFT_ALIGN_LSB not defined

    # Use native byte order for the destination buffer.
    if ENDIAN_BOM == 0x04030201
        # Little-endian byte order.
        cam[PHX_DST_ENDIAN] = PHX_DST_LITTLE_ENDIAN
    elseif ENDIAN_BOM == 0x01020304
        # Big-endian byte order.
        cam[PHX_DST_ENDIAN] = PHX_DST_BIG_ENDIAN
    end

    return nothing # FIXME: return PXH_OK?
end

"""
    equivalentformat(pixfmt) -> (srcfmt, srcdepth, dstfmt)

yields values of parameters `PHX_CAM_SRC_COL`, `PHX_CAM_SRC_DEPTH`
and `PHX_DST_FORMAT` corresponding to the camera pixel format `pixfmt`
which is one of `PIXEL_FORMAT_MONO8`, `PIXEL_FORMAT_MONO10`,
`PIXEL_FORMAT_BAYERGR8` or `PIXEL_FORMAT_BAYERGR10`.

"""
function equivalentformat(pixfmt::Integer)
    local srcfmt::ParamValue, srcdepth::ParamValue
    local dstfmt::ParamValue
    if pixfmt == PIXEL_FORMAT_MONO8
        srcfmt   = PHX_CAM_SRC_MONO
        srcdepth = 8
        dstfmt   = PHX_DST_FORMAT_Y8
    elseif pixfmt == PIXEL_FORMAT_MONO10
        srcfmt   = PHX_CAM_SRC_MONO
        srcdepth = 10
        dstfmt   = PHX_DST_FORMAT_Y16
    elseif pixfmt == PIXEL_FORMAT_BAYERGR8
        srcfmt   = PHX_CAM_SRC_BAY_RGGB
        srcdepth = 8
        dstfmt   = PHX_DST_FORMAT_BAY8
    elseif pixfmt == PIXEL_FORMAT_BAYERGR10
        srcfmt   = PHX_CAM_SRC_BAY_RGGB
        srcdepth = 10
        dstfmt   = PHX_DST_FORMAT_BAY16
    else
        error("unknown pixel format 0x", string(pixfmt, base=16))
    end
    return (srcfmt, srcdepth, dstfmt)
end

equivalentformat(cam::Camera{MikrotronMC408xModel}) =
    equivalentformat(cam[PIXEL_FORMAT])

"""
    setactiveregion!(cam, width, height)

set the dimensions of the active region for camera `cam`.  This is a low-level
method which does not check its arguments.  It is equivalent to:

    cam[PHX_CAM_ACTIVE_XOFFSET] = 0
    cam[PHX_CAM_ACTIVE_YOFFSET] = 0
    cam[PHX_CAM_ACTIVE_XLENGTH] = width
    cam[PHX_CAM_ACTIVE_YLENGTH] = height
    cam[PHX_CAM_XBINNING]       = 1
    cam[PHX_CAM_YBINNING]       = 1

See also: [`setsourceregion!`](@ref).
"""
function setactiveregion!(cam::Camera{MikrotronMC408xModel},
                          width::Integer, height::Integer)
    cam[PHX_CAM_ACTIVE_XOFFSET] = 0
    cam[PHX_CAM_ACTIVE_YOFFSET] = 0
    cam[PHX_CAM_ACTIVE_XLENGTH] = width
    cam[PHX_CAM_ACTIVE_YLENGTH] = height
    cam[PHX_CAM_XBINNING]       = 1
    cam[PHX_CAM_YBINNING]       = 1
    nothing
end

"""
    setsourceregion!(cam, xoff, yoff, width, height)

set the offsets and dimensions of the *source* region for camera `cam`.  This
is a low-level method which does not check its arguments.  It is equivalent to:

    cam[PHX_ROI_SRC_XOFFSET] = xoff
    cam[PHX_ROI_SRC_YOFFSET] = yoff
    cam[PHX_ROI_XLENGTH]     = width
    cam[PHX_ROI_YLENGTH]     = height

The *source* region is the *capture* region, a.k.a. *ROI*, relative to the
*active* region.

See also: [`setcaptureregion!`](@ref), [`setroi!`](@ref).
"""
function setsourceregion!(cam::Camera{MikrotronMC408xModel},
                          xoff::Integer, yoff::Integer,
                          width::Integer, height::Integer)
    cam[PHX_ROI_SRC_XOFFSET] = xoff
    cam[PHX_ROI_SRC_YOFFSET] = yoff
    cam[PHX_ROI_XLENGTH]     = width
    cam[PHX_ROI_YLENGTH]     = height
    nothing
end

_starthook(cam::Camera{MikrotronMC408xModel}) =
    exec(cam, ACQUISITION_START)

_stophook(cam::Camera{MikrotronMC408xModel}) =
    exec(cam, ACQUISITION_STOP)

# Extend method.
getfullwidth(cam::Camera{MikrotronMC408xModel}) =
    Int(cam[SENSOR_WIDTH])

# Extend method.
getfullheight(cam::Camera{MikrotronMC408xModel}) =
    Int(cam[SENSOR_HEIGHT])

# Extend method.
function getroi(cam::Camera{MikrotronMC408xModel};
                quiet::Bool = false)

    @assert cam.state > 0

    # The frame grabber receives an image of `camwidth` by `camheight`
    # macro-pixels and all ROI settings on the frame grabber are with
    # macro-pixels.  On the camera, only the width and height of the ROI are
    # expressed in macro-pixels, the offsets and decimation factors are in
    # pixels.

    # Retrieve current settings for the camera active region in macro-pixels.
    xsub      = Int(cam[DECIMATION_HORIZONTAL]) # in pixels
    ysub      = Int(cam[DECIMATION_VERTICAL])   # in pixels
    camxoff   = Int(cam[OFFSET_X])              # in pixels
    camyoff   = Int(cam[OFFSET_Y])              # in pixels
    camwidth  = Int(cam[WIDTH])                 # in macro-pixels
    camheight = Int(cam[HEIGHT])                # in macro-pixels

    # Make sure the active region is correct.
    setactiveregion!(cam, camwidth, camheight)

    # Retrieve current settings for the frame grabber ROI.
    srcxoff   = Int(cam[PHX_ROI_SRC_XOFFSET])   # in macro-pixels
    srcyoff   = Int(cam[PHX_ROI_SRC_YOFFSET])   # in macro-pixels
    width  = Int(cam[PHX_ROI_XLENGTH])       # in macro-pixels
    height = Int(cam[PHX_ROI_YLENGTH])       # in macro-pixels

    # Check settings of the source region and fix them (try to clip first and
    # reset in last resort).
    clip = false
    reset = false
    if srcxoff < 0
        clip = true
        width += srcxoff
        srcxoff = 0
    end
    if srcxoff + width > camwidth
        clip = true
        width = camwidth - srcxoff
    end
    if srcyoff < 0
        clip = true
        height += srcyoff
        srcyoff = 0
    end
    if srcyoff + height > camheight
        clip = true
        height = camheight - srcyoff
    end
    if width < 1 || height < 1
        reset = true
        srcxoff   = 0
        srcyoff   = 0
        width  = camwidth
        height = camheight
    end
    if reset || clip
        quiet || @warn (reset ?
                        "non-overlapping ROI has been reset to active region" :
                        "ROI has been clipped within active region")
        setsourceregion!(cam, srcxoff, srcyoff, width, height)
    end

    # Compute actual ROI and return it.
    xoff = camxoff + srcxoff*xsub
    yoff = camyoff + srcyoff*ysub
    return ROI(xsub, ysub, xoff, yoff, width, height)
end

# Extend method.
function setroi!(cam::Camera{MikrotronMC408xModel}, roi::ROI;
                 quiet::Bool = false)
    # Check arguments and retrieve current camera settings.
    @assert cam.state == 1
    fullwidth  = getfullwidth(cam)
    fullheight = getfullheight(cam)
    checkroi(roi, fullwidth, fullheight)

    # Compute new settings for the camera (to use the smallest active region
    # within constraints) and the source region.
    camxoff, camwidth,  srcxoff = fitroi(roi.xsub, roi.xoff, roi.width,
                                         HORIZONTAL_INCREMENT, fullwidth,
                                         "horizontal")
    camyoff, camheight, srcyoff = fitroi(roi.ysub, roi.yoff, roi.height,
                                         VERTICAL_INCREMENT, fullheight,
                                         "vertical")

    # Compute actual size of ROI.
    width  = min(roi.width,  camwidth  - srcxoff)
    height = min(roi.height, camheight - srcyoff)

    if true
        # FIXME: There is a bug in the frame grabber which yields corrupted
        # images when the size of the ROI does not match what is sent by the
        # camera.
        srcxoff = 0
        srcyoff = 0
        width = camwidth
        height = camheight
    end

    # Fix settings in a specific order such that the actual camera settings are
    # always valid.
    fixcamroi!(cam,
               DECIMATION_HORIZONTAL, roi.xsub,
               OFFSET_X,              camxoff,
               WIDTH,                 camwidth)
    fixcamroi!(cam,
               DECIMATION_VERTICAL,   roi.ysub,
               OFFSET_Y,              camyoff,
               HEIGHT,                camheight)

    # Set frame grabber parameters.
    setactiveregion!(cam, camwidth, camheight)
    setsourceregion!(cam, srcxoff, srcyoff, width, height)
    return nothing
end

"""
    rounddown(a, b)

yields largest multiple of `b` which is smaller or equal `a`.  Arguments and
result are integers, `a` must be nonnegative, `b` must be striclty positive and
result is nonnegative.

See also: [`roundup`](@ref)

"""
rounddown(a::Integer, b::Integer) = div(a, b)*b

"""
    roundup(a, b)

yields smallest multiple of `b` which is larger or equal `a`.  Arguments and
result are integers, `a` must be nonnegative, `b` must be striclty positive and
result is nonnegative.

See also: [`rounddown`](@ref)

"""
roundup(a::Integer, b::Integer) = rounddown((b - one(b)) + a, b)

"""
     divrnd(a, b) -> q

yields `a/b` rounded to the nearest integer.  Arguments are assumed to be
positive integers.  The result is the same as `round(Int, a/b)` but (slightly)
faster.

"""
divrnd(a::Integer, b::Integer) = div(2a + b, 2b)

"""
    fitroi(sub, off, len, inc, lim, dir) -> devoff, devlen, srcoff

yields the offset (in pixels) and the length (in macro-pixels) of the ROI for
the camera and the offset (in macro-pixels) of the ROI for the frame grabber to
fit a 1D ROI where `sub` is the subsampling factor, `off` is the offset of the
ROI relative to the sensor, `len` is the length of the ROI (in macro-pixels),
`inc` is the sensor increment (in pixels) and `lim` is the maximum length (in
pixels) of the device.  Last argument is `"horizontal"` or `"vertical"` and is
used for error messages.

"""
function fitroi(sub::Int, off::Int, len::Int, inc::Int, lim::Int, dir::String)
    # Check consistency of arguments.
    @assert inc ≥ 1
    @assert lim ≥ 1 && rem(lim, inc) == 0
    @assert sub ≥ 1
    @assert len ≥ 1
    @assert off ≥ 0
    @assert off + len*sub ≤ lim

    # We have to make a compromise between exactness of the ROI (which is not
    # always possible if `sub > 1`) and smallness of the size of the region
    # sent by the camera.  If the region sent by the camera is large enough to
    # encompass the requested ROI, then the maximal error on each side of the
    # ROI is striclty less than a macro-pixel because it is always possible to
    # reduce the final image by macro-pixel adjustments.  The strategy is then
    # to choose the smallest possible region on the camera which encompasses
    # the requested ROI and, then, to compute macro-pixel adjustments to
    # approximate the requested ROI.

    # Offset in multiple of `inc` pixels which best fits by below the requested
    # offset.
    devoff = rounddown(off, inc)

    # Length of the region sent by the camera contrained to be a multiple of
    # `inc` and `sub` (hence a multiple of their Least Common Multiple), at or
    # after end of ROI (if possible) and whithin camera limits.
    mul = lcm(inc, sub)
    pixlen = min(roundup(off + len*sub - devoff, mul),
                 rounddown(lim - devoff, mul))

    # Size of region sent by the camera in macro-pixels.
    devlen = div(pixlen, sub)

    # Source offset in macro-pixels to best fits by below the requested ROI.
    srcoff = div(off - devoff, sub)

    return (devoff, devlen, srcoff)
end

# Other possible strategy
# =======================
#
# The region of interest is `roi = (a,b]` in pixels with
#
#     a = off
#     b = off + len⋅sub
#
# which is approximated by:
#
#     ap = k⋅inc + l⋅sub
#     bp = k⋅inc + m⋅lcm(inc,sub) - n⋅sub
#
# where: `k ≥ 0`, `l ≥ 0`, `m > 0` and `n ≥ 0` are all nonnegative integers
# and:
#
#     devoff = k⋅inc
#     devlen = m⋅lcm(inc,sub)/sub
#     srcoff = l
#
# so as to minimize the misfit:
#
#     err = abs(ap - a) + abs(bp - b)
#
# and for the same misfit, the objective is to minimize the number of
# transmitted macro-pixels hence `devlen` (or equivalently `m`).
#
# If `k` and `m` are given, the best `l` and `n` are:
#
#     l = argmin_{l ≥ 0} abs(k⋅inc + l⋅sub - a)
#       = max(0, round(Int, (a - k⋅inc)/sub))
#
#     n = argmin_{n ≥ 0} abs(k⋅inc + m⋅lcm(inc,sub) - n⋅sub - b)
#       = max(0, round(Int, (k⋅inc + m⋅lcm(inc,sub) - b)/sub))
#

"""
    fixcamroi!(cam, subkey, offkey, lenkey, sub, off, len)

applies 1D ROI settings to the camera `cam` in a specific order so that the
configuration is valid at any moment.  Arguments `subkey`, `offkey` and
`lenkey` are the CoaXPress registers corresponding to `sub` the subsampling
factor (in pixels), `off` the offset of the ROI relative to the sensor (in
pixels) and `len` the length of the ROI (in macro-pixels).

"""
function fixcamroi!(cam::Camera{MikrotronMC408xModel},
                    subkey::Register, sub::Int,
                    offkey::Register, off::Int,
                    lenkey::Register, len::Int)
    oldsub = Int(cam[subkey])
    oldoff = Int(cam[offkey])
    oldlen = Int(cam[lenkey])
    if oldsub > sub
        cam[subkey] = sub
    end
    if oldoff > off
        cam[offkey] = off
    end
    if oldlen != len
        cam[lenkey] = len
    end
    if oldsub < sub
        cam[subkey] = sub
    end
    if oldoff < off
        cam[offkey] = off
    end
    return nothing
end

# Extend method.
function supportedpixelformats(cam::Camera{MikrotronMC408xModel}, buf::Bool)
    if buf
        return CAPTURE_FORMATS
    else
        format = cam[PIXEL_FORMAT]
        if (format == PIXEL_FORMAT_MONO8 ||
            format == PIXEL_FORMAT_MONO10)
            return Union{Monochrome{8},Monochrome{10}}
        elseif (format == PIXEL_FORMAT_BAYERGR8 ||
                format == PIXEL_FORMAT_BAYERGR10)
            return Union{BayerGRBG{8},BayerGRBG{10}}
        else
            error("unexpected pixel format!")
        end
    end
end

# Extend method.
function getpixelformat(cam::Camera{MikrotronMC408xModel})
    fmt = cam[PIXEL_FORMAT]
    return (fmt == PIXEL_FORMAT_MONO8     ? Monochrome{8}  :
            fmt == PIXEL_FORMAT_MONO10    ? Monochrome{10} :
            fmt == PIXEL_FORMAT_BAYERGR8  ? BayerGRBG{8}   :
            fmt == PIXEL_FORMAT_BAYERGR10 ? BayerGRBG{10}  :
            error("unexpected pixel format!"))
end

# Extend method.
function setpixelformat!(cam::Camera{MikrotronMC408xModel},
                         ::Type{T}) where {T <: PixelFormat}
    # Check state.
    if cam.state != 1
        if cam.state == 0
            error("camera not yet open")
        elseif cam.state == 2
            error("acquisition is running")
        else
            error("camera instance corrupted")
        end
    end

    # Determine best matching formats.
    oldfmt = cam[PIXEL_FORMAT]
    newfmt = guesspixelformat(oldfmt, T)
    srccol, srcdepth, dstfmt = equivalentformat(newfmt)

    # Apply the settings.
    if newfmt != oldfmt
        cam[PIXEL_FORMAT] = newfmt
    end
    cam[PHX_CAM_SRC_COL]   = srccol
    cam[PHX_CAM_SRC_DEPTH] = srcdepth
    cam[PHX_DST_FORMAT]    = dstfmt
    return nothing
end

# Determine camera pixel format.
function guesspixelformat(oldfmt::Integer,
                          ::Type{C}) where {C <:PixelFormat}
    if C <: Monochrome
        if oldfmt != PIXEL_FORMAT_MONO8 && oldfmt != PIXEL_FORMAT_MONO10
            throw(ArgumentError("not a monochrome camera"))
        end
        if C == Monochrome{8}
            return PIXEL_FORMAT_MONO8
        elseif C == Monochrome{10}
            return PIXEL_FORMAT_MONO10
        elseif C != Monochrome
            throw(ArgumentError("unsupported number of bits per pixel"))
        end
    elseif C <: ColorFormat
        if oldfmt != PIXEL_FORMAT_BAYERGR8 && oldfmt != PIXEL_FORMAT_BAYERGR10
            throw(ArgumentError("not a color Bayer camera"))
        end
        if (C == ColorFormat{8} ||
            C == BayerFormat{8} ||
            C == BayerGRBG{8})
            return PIXEL_FORMAT_BAYERGR8
        elseif (C == ColorFormat{10} ||
                C == BayerFormat{10} ||
                C == BayerGRBG{10})
            return PIXEL_FORMAT_BAYERGR10
        elseif (C != ColorFormat &&
                C != BayerFormat &&
                C != BayerGRBG)
            throw(ArgumentError("unsupported number of bits per pixel or color format"))
        end
    end
end

getframespersecond(cam::Camera{MikrotronMC408xModel}) =
    convert(Float64, cam[ACQUISITION_FRAME_RATE])

getexposuretime(cam::Camera{MikrotronMC408xModel}) =
    convert(Float64, cam[EXPOSURE_TIME]/1_000_000)

# Extend method.
getspeed(cam::Camera{MikrotronMC408xModel}) =
    (getframespersecond(cam), getexposuretime(cam))

# Extend method.
function setspeed!(cam::Camera{MikrotronMC408xModel},
                   fps::Float64, exp::Float64)
    checkspeed(cam, fps, exp)
    newfps = round(UInt32, fps)
    newexp = round(UInt32, exp*1_000_000)
    oldfps = cam[ACQUISITION_FRAME_RATE]
    oldexp = cam[EXPOSURE_TIME]
    if newfps < oldfps
        cam[ACQUISITION_FRAME_RATE] = newfps
    end
    if newexp != oldexp
        cam[EXPOSURE_TIME] = newexp
    end
    if newfps > oldfps
        cam[ACQUISITION_FRAME_RATE] = newfps
    end
    nothing
end

function setspeed!(cam::Camera{MikrotronMC408xModel},
                   fps::Real, exp::AbstractString)
    setspeed!(cam, fps, Symbol(exp))
end

function setspeed!(cam::Camera{MikrotronMC408xModel},
                   fps::Real, exp::Symbol)
    if exp == :max
        setspeed!(cam, fps, 1e-6)
        setspeed!(cam, fps, 1e-6*cam[EXPOSURE_TIME_MAX])
    else
        error("exposure can be a value of `:max`")
    end
    nothing
end

# Extend method.
getgain(cam::Camera{MikrotronMC408xModel}) =
    convert(Float64, cam[GAIN])

# Extend method.
setgain!(cam::Camera{MikrotronMC408xModel}, gain::Float64) =
    setifneeded!(cam, GAIN, round(UInt32, gain))

# Extend method.
getbias(cam::Camera{MikrotronMC408xModel}) =
    convert(Float64, cam[BLACK_LEVEL])

# Extend method.
setbias!(cam::Camera{MikrotronMC408xModel}, bias::Float64) =
    setifneeded!(cam, BLACK_LEVEL, round(UInt32, bias))

# Extend method.
getgamma(cam::Camera{MikrotronMC408xModel}) =
    convert(Float64, cam[GAMMA])

# Extend method.
setgamma!(cam::Camera{MikrotronMC408xModel}, gamma::Float64) =
    setifneeded!(cam, GAMMA, gamma)

function setifneeded!(cam::Camera{MikrotronMC408xModel},
                      reg::RegisterValue{T,ReadWrite},
                      val) where {T<:Integer}
    newval = round(T, val)
    if cam[reg] != newval
        cam[reg] = newval
    end
    nothing
end

function setifneeded!(cam::Camera{MikrotronMC408xModel},
                      reg::RegisterValue{T,ReadWrite},
                      val) where {T<:AbstractFloat}
    newval = convert(T, val)
    if cam[reg] != newval
        cam[reg] = newval
    end
    nothing
end

# This overloading of the method is to treat specifically certain problematic
# parameters such as the pixel format.
function setparam!(cam::Camera{MikrotronMC408xModel},
                   key::RegisterValue{T,A},
                   val::T) :: T where {T<:Real,A<:Writable}
    #
    # Unfortunately, setting some parameters (as the pixel format or the gamma
    # correction) returns an error with an absurd code
    # (`PHX_ERROR_MALLOC_FAILED`) even if the value has been correctly set.
    # The error cannot be just ignored as further reads of registers yield
    # wrong values.  The strategy is to close and re-open the camera when such
    # an error occurs which solves the problem in practice to the cost of the
    # time spent to close and re-open (0.4 sec.).  To avoid alarming the user,
    # printing of error messages is disabled during this process.
    #
    # Other (unsuccessful) strategies have been tried:
    # - Calling PHX_ControlReset yields a PHX_ERROR_BAD_HANDLE error.
    # - Re-reading the register until its value is OK (this does not work for
    #   write-only registers).
    #
    errmode = printerror(false) # temporarily switch off reporting of errors
    buf = Ref{T}(val)
    status = _setparam!(cam, key, buf)
    printerror(errmode) # restore previous mode
    if status == PHX_ERROR_MALLOC_FAILED && 1 ≤ cam.state ≤ 2
        if cam.state == 2
            @warn "Acquisition will be aborted due to setting a bogus parameter"
        end
        close(cam)
        open(cam; quiet=true)
        if _getparam(cam, key, buf) != PHX_OK || buf[] != val
            error("failed to change register value")
        end
    elseif status != PHX_OK
        _printlasterror()
        checkstatus(status)
    end
    return val
end

"""
    setfixedpatternnoisereduction!(cam, val)

switches the Fixed Pattern Noise Reduction feature on or off for camera `cam`.

    getfixedpatternnoisereduction(cam)

yields the current setting.

"""
setfixedpatternnoisereduction!(cam::Camera{MikrotronMC408xModel}, val::Bool) =
    cam[FIXED_PATTERN_NOISE_REDUCTION] = (val ? one(UInt32) : zero(UInt32))

getfixedpatternnoisereduction(cam::Camera{MikrotronMC408xModel}) =
    (cam[FIXED_PATTERN_NOISE_REDUCTION] != zero(UInt32))

@doc @doc(setfixedpatternnoisereduction!) getfixedpatternnoisereduction

"""
    setinfofieldframecounter!(cam, val)

enables/disables the *frame counter* info field in the images captured by
camera `cam`.  If enabled, the Frame Counter info field overwrites the 1-st to
4-th pixels of the images as a 32-bit unsigned integer.

    getinfofieldframecounter(cam)

yields the current setting.

    getinfofieldframecounter(cam, i)

yields the value of the frame counter info field in `i`-th captured image
buffer.

"""
setinfofieldframecounter!(cam::Camera{MikrotronMC408xModel}, val::Bool) =
    cam[INFO_FIELD_FRAME_COUNTER_ENABLE] = (val ? one(UInt32) : zero(UInt32))

getinfofieldframecounter(cam::Camera{MikrotronMC408xModel}) =
    (cam[INFO_FIELD_FRAME_COUNTER_ENABLE] != zero(UInt32))

getinfofieldframecounter(cam::Camera{MikrotronMC408xModel}, i::Integer) =
    convert(Int, getinfofieldvalue(UInt32, cam.imgs[i], 0))

@doc @doc(setinfofieldframecounter!) getinfofieldframecounter

"""
    setinfofieldtimestamp!(cam, val)

enables/disables the *time stamp* info field in the images captured by camera
`cam`.  The frequency of the time stamp counter amounts to 25 MHz (period = 40
nanoseconds) and is stored as a 32-bit unsigned integer.  If enabled, the time
stamp info field overwrites the 5-th to 8-th pixels of the images.

    getinfofieldtimestamp(cam)

yields the current setting.

    getinfofieldtimestamp(cam, i)

yields the value (in seconds) of the time stamp info field in `i`-th captured
image buffer.

"""
setinfofieldtimestamp!(cam::Camera{MikrotronMC408xModel}, val::Bool) =
    cam[INFO_FIELD_TIME_STAMP_ENABLE] = (val ? one(UInt32) : zero(UInt32))

getinfofieldtimestamp(cam::Camera{MikrotronMC408xModel}) =
    (cam[INFO_FIELD_TIME_STAMP_ENABLE] != zero(UInt32))

# The frequency of the time stamp counter amounts to 25 MHz (period = 40 nanoseconds).
getinfofieldtimestamp(cam::Camera{MikrotronMC408xModel}, i::Integer) =
    getinfofieldvalue(UInt32, cam.imgs[i], 4)*40e-9

getinfofieldvalue(::Type{UInt16}, img::Array{UInt8,2}, off::Int) =
    ((convert(UInt16, img[off+2]) <<  8) | convert(UInt16, img[off+1]))

getinfofieldvalue(::Type{UInt32}, img::Array{UInt8,2}, off::Int) =
    ((convert(UInt32, img[off+4]) << 24) |
     (convert(UInt32, img[off+3]) << 16) |
     (convert(UInt32, img[off+2]) <<  8) | convert(UInt32, img[off+1]))

getinfofieldvalue(::Type{UInt32}, img::Array{UInt16,2}, off::Int) =
    ((convert(UInt32, img[off+4] & 0x00FF) << 24) |
     (convert(UInt32, img[off+3] & 0x00FF) << 16) |
     (convert(UInt32, img[off+2] & 0x00FF) <<  8) | convert(UInt32, img[off+1] & 0x00FF))

@doc @doc(setinfofieldtimestamp!) getinfofieldtimestamp

"""
    setinfofieldroi!(cam, val)

enables/disables the *ROI* info field in the images captured by camera `cam`.
If enabled, the ROI info field overwrites the 9-th to 16-th pixels of the
images.  The ROI info is only available in 8-bit pixel format mode.

    getinfofieldroi(cam)

yields the current setting.

    getinfofieldroi(cam, i)

yields the value of the ROI info field in `i`-th captured image buffer.

"""
setinfofieldroi!(cam::Camera{MikrotronMC408xModel}, val::Bool) =
    cam[INFO_FIELD_ROI_ENABLE] = (val ? one(UInt32) : zero(UInt32))

getinfofieldroi(cam::Camera{MikrotronMC408xModel}) =
    (cam[INFO_FIELD_ROI_ENABLE] != zero(UInt32))

function getinfofieldroi(cam::Camera{MikrotronMC408xModel}, i::Integer)
    img = cam.imgs[i]
    xoff   = convert(Int, getinfofieldvalue(UInt16, img,  8))
    width  = convert(Int, getinfofieldvalue(UInt16, img, 10))
    yoff   = convert(Int, getinfofieldvalue(UInt16, img, 12))
    height = convert(Int, getinfofieldvalue(UInt16, img, 14))
    return (xoff, yoff, width, height)
end

@doc @doc(setinfofieldroi!) getinfofieldroi

"""
    setfiltermode!(cam, val)

enables/disables the image filter mode of the camera `cam`.

    getfiltermode(cam)

yiels the current setting.

"""
function setfiltermode!(cam::Camera{MikrotronMC408xModel}, val::Bool)
    if val
        format = cam[PIXEL_FORMAT]
        if (format == PIXEL_FORMAT_BAYERGR8 ||
            format == PIXEL_FORMAT_BAYERGR10)
            cam[FILTER_MODE] = FILTER_MODE_COLOR
        elseif (format == PIXEL_FORMAT_MONO8 ||
                format == PIXEL_FORMAT_MONO10)
            cam[FILTER_MODE] = FILTER_MODE_MONO
        else
            error("unknown pixel format (0x$(string(format, base=16)))")
        end
    else
        cam[FILTER_MODE] = FILTER_MODE_RAW
    end
    return nothing
end

getfiltermode(cam::Camera{MikrotronMC408xModel}) =
    (cam[FILTER_MODE] != FILTER_MODE_RAW)

end # module


text:
isdefined(:MNIST) || include(Pkg.dir("Knet/examples/mnist.jl"))
using Knet, BenchmarkTools
using MNIST: weights, minibatch, xtrn, ytrn
using Knet: gpusync, reluback, xentloss, xentback
using AutoGrad: unbroadcast

atype  = KnetArray{Float32}
hidden = [64]
batch  = 100
data   = minibatch(xtrn, ytrn, batch, atype=atype)
model  = weights(64, atype=atype)

forw = quote
    a1 = w[1]*x                         # 01 0.56 0.56
    a2 = w[2].+a1                       # 02 0.60 0.04
    a3 = relu(a2)                       # 03 0.63 0.03  0.64 0.04
    #a3 = max(0,a2)
    a4 = w[3]*a3                        # 04 0.76 0.13  0.77 0.13
    a5 = w[4].+a4                       # 05 0.79 0.03  0.80 0.04
    #a6 = y - a5
    a6 = logp(a5,1)                     # 06 1.02 0.23  1.03 0.23
    #a6 = a5 .- log(sum(exp(a5),1))
    #a7 = a6.*a6
    #a7 = a6.^2
    #a7 = abs2(a6)
    a7 = y .* a6                        # 07 1.06 0.04  1.06 0.03
    #a7 = sumabs2(a6)
    a8 = sum(a7)                        # 08 1.32 0.26  1.35 0.28
    #a6 = xentloss(a5,y,1)              # xx 1.32 0.52
end                                     # 20 1.43 0.11  1.42 0.10 (record)

back = quote
    d7 = ones(a7)                       # 09 1.35 0.03  1.39 0.07
    d6 = unbroadcast(a6, y .* d7)       # 10 1.39 0.04  1.40 0.01
    d5 = d6 - exp(a6) .* sum(d6,1)	# 11 1.57 0.18  1.57 0.17
#    d5 = xentback(a5,y,1)	        # xx 1.63 0.31
    dw4 = unbroadcast(w[4], d5)         # 12 1.68 0.11  1.63 0.06
    d4 = unbroadcast(a4, d5)            # 13 1.69 0.01  1.63 0.00
    dw3 = A_mul_Bc(d4, a3)              # 14 1.85 0.16  1.79 0.16
    d3 = Ac_mul_B(w[3], d4)             # 15 1.96 0.11  1.89 0.10  1.90 0.11
    d2 = reluback(d3,a3)                # 16 1.99 0.03  1.99 0.10  1.93 0.03
    dw2 = unbroadcast(w[2], d2)         # 17 2.04 0.05  2.05 0.06  1.98 0.05
    d1 = unbroadcast(a1, d2)            # 18 2.04 0.00  2.05 0.00  1.99 0.01
    dw1 = A_mul_Bc(d1, x)               # 19 2.22 0.18  2.23 0.18  2.17 0.18
end                                     # 21 2.48 0.26  2.52 0.29  2.52 0.35 (grad)

forw = filter(p->isa(p,Expr), forw.args)
back = filter(p->isa(p,Expr), back.args)
fun = Any[]

for i=1:length(forw)
    f = Symbol("f$i")
    eval(Expr(:function, :($f(w,x,y)), Expr(:block,forw[1:i]...)))
    push!(fun, eval(f))
end
fforw = fun[end]
for i=1:length(back)
    j=i+length(forw)
    f = Symbol("f$j")
    eval(Expr(:function, :($f(w,x,y)), Expr(:block,forw...,back[1:i]...)))
    push!(fun, eval(f))
end

push!(fun,(w,x,y)->AutoGrad.forward_pass(fforw,(w,x,y),(),1))
push!(fun,grad(fforw))

function loop(f,w,d,t)
    for i in 1:t
        for (x,y) in d
            f(w,x,y)
        end
    end
end

@time for i=1:length(fun)
    loop(fun[i],model,data,1)   # to compile
end

function profmlp(fun,model,data,range=1:length(fun); samples=5, seconds=30, o...)
    for i in range
        f = fun[i]
        b = @benchmarkable (loop($f,$model,$data,10);gpusync())
        t = run(b; samples=samples, seconds=seconds, o...).times ./ 10^9
        println((:f,i,:n,length(t),:min,minimum(t),:med,median(t),:avg,mean(t),:max,maximum(t)))
    end
end

# (:f,1,:n,5,:min,0.674071033,:med,0.674303586,:avg,0.6742660326000001,:max,0.674379312)
# (:f,2,:n,5,:min,0.71386119,:med,0.714166844,:avg,0.714133768,:max,0.714420962)
# (:f,3,:n,5,:min,0.749161714,:med,0.749578344,:avg,0.7576167076000001,:max,0.789836767)
# (:f,4,:n,5,:min,0.812147525,:med,0.812692507,:avg,0.81255644,:max,0.812902328)
# (:f,5,:n,5,:min,0.850032182,:med,0.850729729,:avg,0.8505739558000001,:max,0.850794662)
# (:f,6,:n,5,:min,0.889973775,:med,0.890434318,:avg,0.8905268512,:max,0.890990265)
# (:f,7,:n,5,:min,0.916586753,:med,0.91684639,:avg,0.9174720356,:max,0.919871972)
# (:f,8,:n,5,:min,1.218039367,:med,1.223258585,:avg,1.222431379,:max,1.226436017)
# (:f,9,:n,5,:min,1.433074136,:med,1.437228564,:avg,1.4373717761999998,:max,1.441818492)
# (:f,10,:n,5,:min,2.967745663,:med,2.975174783,:avg,2.9851800988,:max,3.029562013)

ext: jl


text:
# Utilities for working with HSA signals

mutable struct HSASignal
    signal::Ref{Signal}
end

function HSASignal(init::Integer=1)
    signal = HSASignal(Ref{Signal}())
    HSA.signal_create(Int64(init), 0, C_NULL, signal.signal)
    hsaref!()
    finalizer(signal) do signal
        HSA.signal_destroy(signal.signal[]) |> check
        hsaunref!()
    end
    return signal
end
HSASignal(signal::HSA.Signal) = HSASignal(Ref(signal))

Adapt.adapt_structure(::Adaptor, sig::HSASignal) = sig.signal[]

"""
    Base.wait(signal::HSASignal; soft=true, minlat=0.000001)

Waits on an `HSASignal` to decrease below 1. If `soft=true` (default), uses
sleep polling; otherwise uses HSA's signal waiter. `minlat` sets the minimum
latency for the software waiter; lower values can decrease latency at the cost
of increased polling load. `timeout`, if not `nothing`, sets the timeout for
software waiting, after which the call will error.
"""
function Base.wait(signal::HSASignal; soft=true, minlat=0.000001, timeout=nothing)
    if soft
        start_time = time_ns()
        while true
            value = HSA.signal_load_scacquire(signal.signal[])
            if value < 1
                return
            end
            now_time = time_ns()
            diff_time = (now_time - start_time) / 10^9
            if timeout !== nothing && diff_time > timeout
                error("Timeout while waiting on signal")
            end
            if minlat < 0.001 && diff_time < 10^3
                # Use Libc.systemsleep for higher precision in the microsecond range
                Libc.systemsleep(minlat)
                yield()
            else
                sleep(minlat)
            end
        end
    else
        # Wait on the dispatch completion signal until the kernel is finished
        HSA.signal_wait_acquire(signal.signal[],
                                HSA.SIGNAL_CONDITION_LT, 1, typemax(UInt64),
                                HSA.WAIT_STATE_BLOCKED)
    end
end
Base.wait(signal::HSA.Signal; kwargs...) = wait(HSASignal(signal); kwargs...)

Base.notify(signal::HSASignal) = HSA.signal_store_screlease(signal.signal[], 0)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EncryptionScopeProperties <: SwaggerModel
    source::Any # spec type: Union{ Nothing, String } # spec name: source
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    creationTime::Any # spec type: Union{ Nothing, DateTime } # spec name: creationTime
    lastModifiedTime::Any # spec type: Union{ Nothing, DateTime } # spec name: lastModifiedTime
    keyVaultProperties::Any # spec type: Union{ Nothing, EncryptionScopeKeyVaultProperties } # spec name: keyVaultProperties

    function EncryptionScopeProperties(;source=nothing, state=nothing, creationTime=nothing, lastModifiedTime=nothing, keyVaultProperties=nothing)
        o = new()
        validate_property(EncryptionScopeProperties, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        validate_property(EncryptionScopeProperties, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(EncryptionScopeProperties, Symbol("creationTime"), creationTime)
        setfield!(o, Symbol("creationTime"), creationTime)
        validate_property(EncryptionScopeProperties, Symbol("lastModifiedTime"), lastModifiedTime)
        setfield!(o, Symbol("lastModifiedTime"), lastModifiedTime)
        validate_property(EncryptionScopeProperties, Symbol("keyVaultProperties"), keyVaultProperties)
        setfield!(o, Symbol("keyVaultProperties"), keyVaultProperties)
        o
    end
end # type EncryptionScopeProperties

const _property_map_EncryptionScopeProperties = Dict{Symbol,Symbol}(Symbol("source")=>Symbol("source"), Symbol("state")=>Symbol("state"), Symbol("creationTime")=>Symbol("creationTime"), Symbol("lastModifiedTime")=>Symbol("lastModifiedTime"), Symbol("keyVaultProperties")=>Symbol("keyVaultProperties"))
const _property_types_EncryptionScopeProperties = Dict{Symbol,String}(Symbol("source")=>"String", Symbol("state")=>"String", Symbol("creationTime")=>"DateTime", Symbol("lastModifiedTime")=>"DateTime", Symbol("keyVaultProperties")=>"EncryptionScopeKeyVaultProperties")
Base.propertynames(::Type{ EncryptionScopeProperties }) = collect(keys(_property_map_EncryptionScopeProperties))
Swagger.property_type(::Type{ EncryptionScopeProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EncryptionScopeProperties[name]))}
Swagger.field_name(::Type{ EncryptionScopeProperties }, property_name::Symbol) =  _property_map_EncryptionScopeProperties[property_name]

const _allowed_EncryptionScopeProperties_source = ["Microsoft.Storage", "Microsoft.KeyVault"]

const _allowed_EncryptionScopeProperties_state = ["Enabled", "Disabled"]

function check_required(o::EncryptionScopeProperties)
    true
end

function validate_property(::Type{ EncryptionScopeProperties }, name::Symbol, val)
    if name === Symbol("source")
        Swagger.validate_param(name, "EncryptionScopeProperties", :enum, val, _allowed_EncryptionScopeProperties_source)
    end
    if name === Symbol("state")
        Swagger.validate_param(name, "EncryptionScopeProperties", :enum, val, _allowed_EncryptionScopeProperties_state)
    end
end


text:
export  subpopulationPopReport

"""
```
subpopulationPopReport(
    mpSim::MPsim,
    timeGrid::Vector{Float64},
    subpopulations::Subpopulation... )
```
This function generates reports for the evolution of the population of the valid subpopulations in `subpopulations` of the manpower simulation `mpSim`. The reports are generated on the grid of time points `timeGrid`.

This function will issue a warning and not generate any report in the following two cases:
1. There are no positive time points in the time grid;
2. None of the entered subpopulations are valid subpopulations in the simulation, where a subpopulation is invalid if its souce nodes doesn't exist in the simulation.

This function returns a `DataFrame`, with the first column the time points and the other columns corresponding to the population counts at each time point for each subpopulation. In case the function issues a warning, its return value will be an empty `DataFrame`.
"""
function subpopulationPopReport( mpSim::MPsim, timeGrid::Vector{Float64},
    subpopulations::Subpopulation... )::DataFrame

    result = DataFrame()

    if now( mpSim ) == 0
        @warn "Simulation hasn't started yet, can't make report."
        return result
    end  # if now( mpSim ) == 0

    timeGrid = timeGrid[0.0 .<= timeGrid .<= now( mpSim )]
    timeGrid = unique( sort( timeGrid, rev = true ) )

    if isempty( timeGrid )
        @warn "No valid time points in time grid, cannot generate report."
        return result
    end  # if isempty( timeGrid )

    if timeGrid[end] > 0.0
        push!( timeGrid, 0.0 )
    end  # if timeGrid[end] > 0.0

    reverse!( timeGrid )

    subpopulations = filter( collect( subpopulations ) ) do subpopulation
        sourceNode = subpopulation.sourceNode
        return ( sourceNode == "active" ) ||
            haskey( mpSim.baseNodeList, sourceNode ) ||
            haskey( mpSim.compoundNodeList, sourceNode )
    end  # filter( ... ) do subpopulation

    if isempty( subpopulations )
        return result
    end  # if isempty( subpopulations )

    counts = zeros( Int, length( timeGrid ), length( subpopulations ) )

    for ii in eachindex( timeGrid )
        counts[ii, :] = length.( getSubpopulationAtTime( mpSim,
            timeGrid[ii], subpopulations ) )
    end  # for ii in eachindex( timeGrid )

    return DataFrame( hcat( timeGrid, counts ), vcat( :timePoint,
        map( subpopulation -> Symbol( subpopulation.name ), subpopulations ) ) )

end  # subpopulationPopReport( mpSim, timeGrid, subpopulations )

"""
```
subpopulationPopReport(
    mpSim::MPsim,
    timeRes::Real,
    subpopulations::Subpopulation... )
```
This function generates reports for the evolution of the population of the valid subpopulations in `subpopulations` of the manpower simulation `mpSim`.  The reports are generated on a grid of time points with resolution `timeRes`.

This function will issue a warning and not generate any report in the following two cases:
1. There are no positive time points in the time grid;
2. None of the entered subpopulations are valid subpopulations in the simulation, where a subpopulation is invalid if its souce nodes doesn't exist in the simulation.

This function returns a `DataFrame`, with the first column the time points and the other columns corresponding to the population counts at each time point for each subpopulation. In case the function issues a warning, its return value will be an empty `DataFrame`.
"""
subpopulationPopReport( mpSim::MPsim, timeRes::Real,
    subpopulations::Subpopulation... )::DataFrame =
    subpopulationPopReport( mpSim, generateTimeGrid( mpSim, timeRes ),
    subpopulations... )

include( joinpath( repPrivPath, "subpopulationpop.jl" ) )


text:
using ShaleDrillingModel
using ShaleDrillingData
using Test
using StatsFuns
using FileIO
using Interpolations
using Statistics
using SparseArrays
using LinearAlgebra
using BenchmarkTools
using MarkovTransitionMatrices
using Plots
using SparseArrays

# -----------------------------------------------------
# VF setup
# -----------------------------------------------------

# initial parameters
FF = :cheb3_cost_restr
θfull = [-3.58283, -13.9518*0.3, -12.5996*0.3, 11.6749*0.3, -6.40167*0.3, 2.0126*0.3, -1.0, -1.17545, 0.673459,]

# discount
β = (1.0234188 / 1.125) ^ (1.0/4.)

# problem sizes
wp = well_problem(8, 8, 30, 20, 8)

# -----------------------------------------------------
# Price matrix
# -----------------------------------------------------

minp = 1e-4

# make deviations
zrandwalk(x1::Real, mu::Real, σ::Real) = (x1 - mu) / σ
Elogc(x0::Real) = x0 + 0.017992141 - 0.004975763 * exp(x0)

N_YEARS = 14
nlogp = 17
nlogc = 17
nlogσ = 9

extrema_logp = [0.7923374  2.443669]
extrema_logc = [0.5520625  1.446025]
extrema_logσ = [-2.7975151 -1.579074]
sd_logp = 0.10037167
sd_logc = 0.05040059
sd_logσ = 0.18021598
cσ_to_pσ = 0.02777901 / 0.08488665

# make range for coefs
logσ_space = range(extrema_logσ[1] - 2*sd_logσ, stop=extrema_logσ[2] + 2*sd_logσ, length=nlogσ)
logp_space = range(extrema_logp[1] - log(2.0), stop=extrema_logp[2] + log(2.0), length=nlogp)
logc_space = range(extrema_logc[1] - log(2.0), stop=extrema_logc[2] + log(2.0), length=nlogc)
year_space = LinRange(-1, 1, N_YEARS)  # Using a LinRange is a bit hacky -- internally it goes to a BSpline(Constant()) instead of Bspline(Quadratic()), but then converts to StepRangeLen

# range for i,j block in matrix where block is n x n
blockrange(i,j,n) = ((i-1)*n+1):(i*n),  ((j-1)*n+1):(j*n)

# create log σ discretization
println("Discretizing log σ")
flush(stdout)
Pσ, JN, myΛ, L_p1, approxErr = discreteNormalApprox(logσ_space, logσ_space, (x::Real,st::Real) -> zrandwalk(x,st,sd_logσ), 7)
all(sum(Pσ,dims=2) .≈ 1.0) || throw(error("each row of π must sum to 1"))

# match moments for logp
p_from_sigma(logσ::Real) = Dict( [:P, :JN, :Λ, :L, :approxErr,] .=> discreteNormalApprox(logp_space, logp_space, (x1::Real,x0::Real) -> zrandwalk(x1, x0, sd_logp), 11) )
c_from_sigma(logσ::Real) = Dict( [:P, :JN, :Λ, :L, :approxErr,] .=> discreteNormalApprox(logc_space, logc_space, (x1::Real,x0::Real) -> zrandwalk(x1, Elogc(x0), sd_logc), 11) )
P_year = Matrix(Tridiagonal(zeros(N_YEARS-1), vcat(fill(0.75,N_YEARS-1), 1.0), fill(0.25, N_YEARS-1)))

# println("Discretizing log p and log c")
# ps_from_sigma = [p_from_sigma(logσ) for logσ in logσ_space]
# cs_from_sigma = [c_from_sigma(logσ) for logσ in logσ_space]
#
# # Fill giant transition matrix
# println("Filling transition matrix")
# P = Matrix{eltype(logp_space)}(undef, nlogp*nlogc*nlogσ, nlogp*nlogc*nlogσ)
# for (j, logσj) in enumerate(logσ_space)
#   for (i, logσi) in enumerate(logσ_space)
#     P[blockrange(i,j,nlogp*nlogc)...] .= Pσ[i, j] .* ps_from_sigma[i][:P] # kron( cs_from_sigma[i][:P],  ps_from_sigma[i][:P])
#   end
# end

PC = kron( c_from_sigma(log(sd_logc))[:P] ,  p_from_sigma(log(sd_logp))[:P] )
PCY = kron(P_year, c_from_sigma(log(sd_logc))[:P] ,  p_from_sigma(log(sd_logp))[:P] )
PY = kron(P_year, p_from_sigma(log(sd_logp))[:P])

# sparsify transigion
Πp = MarkovTransitionMatrices.sparsify!(PC, minp)
Πpc = MarkovTransitionMatrices.sparsify!(PCY, minp)
Πpyr = MarkovTransitionMatrices.sparsify!(PY, minp)
P = 0.0 # destroy P
PC = 0.0
PCY = 0.0

THIS_Π = Πpc
# zspace = (logp_space, logc_space,) # (logp_space, logc_space, logσ_space,)
# zspace = (logp_space, year_space,) # (logp_space, logc_space, logσ_space,)
zspace = (logp_space, logc_space, year_space,) # (logp_space, logc_space, logσ_space,)

# -----------------------------------------------------
# Types
# -----------------------------------------------------

# royalty rates
royalty_rates = [0.125, 1.0/6., 0.1875, 0.20, 0.225, 0.25]
royalty_types = 1:length(royalty_rates)
extrema(i[2] for i in data_full.ichars) == (1,6,) || throw(error("Don't have royalty rates 1:6!"))

# geology types
n_gt = 11
geology_types = range(1.3430409262656042, stop=5.194950729101042, length=n_gt)

# number of simulations
haltonbases = (2,3)
haltonskip = 5000

# grids for prices & unobserved heterogeneity
nψ = 11
ψspace = range(-4.0, stop=4.0, length=nψ)

# -----------------------------------------------------
# Types
# -----------------------------------------------------

# sparsify transigion
println("Πpyr has $(length(Πpyr.nzval)) values")
println("Πpc has $(length(Πpc.nzval)) values")

# k = 11
# Pxform = reshape(P, nlogp, nlogσ, nlogp, nlogσ)
# heatmap(P)
# heatmap(P[nlogp*k+1:nlogp*(k+1), nlogp*k+1:nlogp*(k+1)])
# heatmap(P[11:35:(35*nlogσ), 11:35:(35*nlogσ)])
# heatmap(Pxform[35,:,35,:])
# dropdims(sum(Pxform, dims=(1,3)), dims=(1,3))./35
# sum(Pxform[35,:,35,:],dims=2)

# plot moments matched for log P
# plot(logp_space, [ps_from_sigma[sig][:L] for sig in logσ_space ], yticks = 1:11, labels = ["$(round(exp(sig); digits=2))" for sig in logσ_space], xlabel="Log p", ylabel="Moments matched")
# savefig("D:/projects/royalty-rates-and-drilling/plots/moments-matched-logponly.pdf")

# -----------------------------------------------------
#    speed tests
# -----------------------------------------------------

function prim_tmp_evs(FFname, zzspace, bigp)
    prim = dcdp_primitives(FFname, β, wp, zzspace, bigp, ψspace)
    tmpv = dcdp_tmpvars(prim)
    evs = dcdp_Emax(prim)
    ShaleDrillingModel.check_size(prim, evs)
    zero!(evs)
    return evs, tmpv, prim
end

dograd = false
pids = [1,]
typidx = (11,3,)
itypes = (geology_types,royalty_rates,)
typs = getindex.(itypes, typidx)

# println("\n\ndoing NO vol")
# tup = prim_tmp_evs((pdat["pspace"],),       pdat["Πp1"])
# println("Πp has $(length(tup[3].Πz.nzval)) values")
# @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=30, maxit1=20, vftol=1e-9)
#
# println("\n\ndoing hi/lo vol")
# tup = prim_tmp_evs((pdat["pspace"],1:2,),   pdat["Πp"])
# println("Πp has $(length(tup[3].Πz.nzval)) values")
# @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=30, maxit1=20, vftol=1e-9)

println("\n\ndoing giant price/vol matrix")
tup = prim_tmp_evs(FF, zspace, Πpc)
println("Πz has $(length(tup[3].Πz.nzval)) values")
zero!(tup[1])

@show @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=20, maxit1=20, vftol=1e-8)
@show @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=30, maxit1=20, vftol=1e-8)
@show @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=40, maxit1=20, vftol=1e-8)
@show @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=50, maxit1=20, vftol=1e-8)

# tup = prim_tmp_evs(:dgt1_restr, zspace, Πp)
# println("Πz has $(length(tup[3].Πz.nzval)) values")
# zero!(tup[1])
# @btime solve_vf_all!(tup..., θfull, typs, dograd; maxit0=40, maxit1=20, vftol=1e-9)

# println("\nSolving for all types")
# sev = SharedEV(pids, tup[3], geology_types[1:3], royalty_rates[1:2])
# @btime serial_solve_vf_all!(sev, tup[2], tup[3], θfull, Val{dograd}; maxit0=30, maxit1=20, vftol=1e-9)


text:
using StanModels, CSV

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), delim=';');

# Define the Stan language model

m10_4sl = "
data {
  int N;
  int<lower=0, upper=1> L[N];
  vector[N] P;
  vector[N] C;

  int<lower=1, upper=N> N_chimps;
  int<lower=1, upper=N_chimps> chimp[N];
}
parameters {
  real a_chimp[N_chimps];
  real bp;
  real bpc;
}
model {
  vector[N] p;
  target += normal_lpdf(a_chimp | 0, 10);
  target += normal_lpdf(bp | 0, 10);
  target += normal_lpdf(bpc | 0, 10);
  for (i in 1:N) p[i] = a_chimp[chimp[i]] + (bp + bpc * C[i]) * P[i];
  target += binomial_logit_lpmf(L | 1, p);
}
generated quantities {
  vector[N] log_lik;
  {
    vector[N] p;
    for(n in 1:N) {
      p[n] = a_chimp[chimp[n]] + (bp + bpc * C[n]) * P[n];
      log_lik[n] = binomial_logit_lpmf(L[n] | 1, p[n]);
    }
  }
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m10.4sl", m10_4sl);

# Input data for cmdstan

m10_4_data = Dict("N" => size(df, 1), "N_chimps" => length(unique(df[!, :actor])),
"chimp" => df[!, :actor], "L" => df[!, :pulled_left],
"P" => df[!, :prosoc_left], "C" => df[!, :condition]);

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m10_4_data);

# Result rethinking

rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
bp    0.84 0.26  0.43  1.26  2271    1
bpC  -0.13 0.29 -0.59  0.34  2949    1

a[1] -0.74 0.27 -1.16 -0.31  3310    1
a[2] 10.88 5.20  4.57 20.73  1634    1
a[3] -1.05 0.28 -1.52 -0.59  4206    1
a[4] -1.05 0.28 -1.50 -0.60  4133    1
a[5] -0.75 0.27 -1.18 -0.32  4049    1
a[6]  0.22 0.27 -0.22  0.65  3877    1
a[7]  1.81 0.39  1.22  2.48  3807    1
";

# Update sections

if !(sample_file == nothing)
  chn = read_samples(sm)

  chn2 = set_section(chn, Dict(
    :parameters => ["bp", "bpc"],
    :pooled => ["a_chimp.$i" for i in 1:7],
    :generated => ["log_lik.$i" for i in 1:504],
    :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__",
      "divergent__", "energy__"]
    )
  )

  describe(chn2)
  describe(chn2, sections=[:pooled])
end


text:
using Documenter, ReverseDiff

makedocs(modules=[ReverseDiff], sitename="http://www.juliadiff.org/ReverseDiff.jl")


text:
@testset verbose=true "indexing" begin
    import SphericalFunctions: WignerHsize, WignerHrange, WignerHindex
    import SphericalFunctions: WignerDsize, WignerDrange, WignerDindex
    import SphericalFunctions: Ysize, Yrange, Yindex, deduce_limits
    import SphericalFunctions: theta_phi

    ell_max = 16
    ell_max_slow = ell_max ÷ 2

    @testset "WignerHrange" begin
        r(mp_max, ell_max) = hcat([
            [ell, mp, m] for ell in 0:ell_max
                for mp in -min(ell, mp_max):min(ell, mp_max)
                    for m in abs(mp):ell
                        ]...)'
        for ell_max in 0:ell_max
            a = WignerHrange(ell_max)  # Implicitly, mp_max=ell_max
            b = r(ell_max, ell_max)
            @test a == b
            for mp_max in 0:ell_max
                a = WignerHrange(mp_max, ell_max)
                b = r(mp_max, ell_max)
                @test a == b
            end
        end
    end

    @testset "WignerHsize" begin
        for ell_max in 0:ell_max
            @test WignerHsize(ell_max) == size(WignerHrange(ell_max, ell_max))[1]
            for mp_max in ell_max
                @test WignerHsize(mp_max, ell_max) == size(WignerHrange(mp_max, ell_max))[1]
            end
        end
    end

    @testset "WignerHindex" begin
        function fold_H_indices(ell, mp, m)
            if m < -mp
                if m < mp
                    return [ell, -mp, -m]
                else
                    return [ell, -m, -mp]
                end
            else
                if m < mp
                    return [ell, m, mp]
                else
                    return [ell, mp, m]
                end
            end
        end

        for ell_max_i in 0:ell_max_slow
            r = WignerHrange(ell_max_i)
            for ell in 0:ell_max_i
                for mp in -ell:ell
                    for m in -ell:ell
                        i = WignerHindex(ell, mp, m)
                        @test r[i, :] == fold_H_indices(ell, mp, m)
                    end
                end
            end
            for mp_max in 0:ell_max_i
                r = WignerHrange(mp_max, ell_max_i)
                for ell in 0:ell_max_i
                    for mp in -ell:ell
                        for m in -ell:ell
                            if abs(mp) > mp_max && abs(m) > mp_max
                                continue
                            end
                            i = WignerHindex(ell, mp, m, mp_max)
                            @test r[i, :] == fold_H_indices(ell, mp, m)
                        end
                    end
                end
            end
        end
    end

    @testset "WignerDrange" begin
        function r(ell_min, mp_max, ell_max)
            a = [
                [ℓ, mp, m]
                for ℓ in ell_min:ell_max
                for mp in -min(ℓ, mp_max):min(ℓ, mp_max)
                for m in -ℓ:ℓ
            ]
            collect(transpose(reshape(collect(Iterators.flatten(a)), (3, length(a)))))
        end

        for ell_max in 0:ell_max_slow÷2
            for ell_min in 0:ell_max
                a = WignerDrange(ell_min, ell_max)  # Implicitly, mp_max=ell_max
                b = r(ell_min, ell_max, ell_max)
                @test a == b
                for mp_max in 0:ell_max
                    a = WignerDrange(ell_min, mp_max, ell_max)
                    b = r(ell_min, mp_max, ell_max)
                    @test a == b
                end
            end
        end
    end

    @testset "WignerDsize" begin
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                for mp_max in 0:ell_max
                    a = WignerDsize(ell_min, mp_max, ell_max)
                    b = size(WignerDrange(ell_min, mp_max, ell_max))[1]
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            for ell_min in [0]
                for mp_max in 0:ell_max
                    a = WignerDsize(ell_min, mp_max, ell_max)
                    #a = WignerDsize_mpmax(ell_max, mp_max)
                    b = size(WignerDrange(ell_min, mp_max, ell_max))[1]
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                for mp_max in [ell_max]
                    a = WignerDsize(ell_min, mp_max, ell_max)
                    # a = WignerDsize_ellmin(ell_min, ell_max)
                    b = size(WignerDrange(ell_min, mp_max, ell_max))[1]
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            for ell_min in [0]
                for mp_max in [ell_max]
                    a = WignerDsize(ell_min, mp_max, ell_max)
                    # a = WignerDsize(ell_max)
                    b = size(WignerDrange(ell_min, mp_max, ell_max))[1]
                    @test a == b
                end
            end
        end
    end

    @testset "WignerDindex" begin
        function r(ell_min, mp_max, ell_max)
            a = [
                [ℓ, mp, m]
                for ℓ in ell_min:ell_max
                for mp in -min(ℓ, mp_max):min(ℓ, mp_max)
                for m in -ℓ:ℓ
            ]
            collect(transpose(reshape(collect(Iterators.flatten(a)), (3, length(a)))))
        end

        for ellmax in 0:ell_max_slow
            r = WignerDrange(0, ellmax)
            for ell in 0:ellmax
                for mp in -ell:ell
                    for m in -ell:ell
                        i = WignerDindex(ell, mp, m)
                        @test r[i, :] == [ell, mp, m]
                    end
                end
            end
            for ell_min in 0:ellmax
                r = WignerDrange(ell_min, ellmax)
                for ell in ell_min:ellmax
                    for mp in -ell:ell
                        for m in -ell:ell
                            i = WignerDindex(ell, mp, m, ell_min)
                            @test r[i, :] == [ell, mp, m]
                        end
                    end
                end
                for mp_max in 0:ellmax
                    r = WignerDrange(ell_min, mp_max, ellmax)
                    for ell in ell_min:ellmax
                        for mp in -min(ell, mp_max):min(ell, mp_max)
                            for m in -ell:ell
                                i = WignerDindex(ell, mp, m, ell_min, mp_max)
                                @test r[i, :] == [ell, mp, m]
                            end
                        end
                    end
                end
            end
        end
    end

    @testset "Yrange" begin
        function r(ell_min, ell_max)
            a = [
                [ℓ, m]
                for ℓ in ell_min:ell_max
                for m in -ℓ:ℓ
            ]
            collect(transpose(reshape(collect(Iterators.flatten(a)), (2, length(a)))))
        end
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                a = Yrange(ell_min, ell_max)
                b = r(ell_min, ell_max)
                @test a == b
            end
        end
    end

    @testset "Ysize" begin
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                a = Ysize(ell_min, ell_max)
                b = size(Yrange(ell_min, ell_max))[1]
                @test a == b
            end
        end
    end

    @testset "deduce_limits" begin
        for ℓmax in 0:4096
            for ℓmin in 0:min(2, ℓmax)
                deduced = deduce_limits(Ysize(ℓmin, ℓmax), nothing)
                @test deduced == (ℓmin, ℓmax)
            end
        end
        for ℓmax in 0:ell_max
            for ℓmin in [0]
                deduced = deduce_limits(Ysize(ℓmin, ℓmax))
                @test deduced == (ℓmin, ℓmax)
            end
        end
        for ℓmax in 0:ell_max
            for ℓmin in 0:ℓmax
                deduced = deduce_limits(Ysize(ℓmin, ℓmax), ℓmin)
                @test deduced == (ℓmin, ℓmax)
            end
        end
    end

    @testset "Yindex" begin
        for ell_max in 0:ell_max
            for ell_min in [0]
                r = Yrange(ell_min, ell_max)
                for ell in ell_min:ell_max
                    for m in -ell:ell
                        i = Yindex(ell, m)
                        @test r[i, :] == [ell, m]
                    end
                end
            end
            for ell_min in 0:ell_max
                r = Yrange(ell_min, ell_max)
                for ell in ell_min:ell_max
                    for m in -ell:ell
                        i = Yindex(ell, m, ell_min)
                        @test r[i, :] == [ell, m]
                    end
                end
            end
        end
    end

end


text:
using StructuralEquationModels, Test, FiniteDiff
# import StructuralEquationModels as SEM
include(
    joinpath(chop(dirname(pathof(StructuralEquationModels)), tail = 3),
    "test/examples/helper.jl")
    )

############################################################################################
### data
############################################################################################

dat = example_data("political_democracy")
dat_missing = example_data("political_democracy_missing")
solution_lav = example_data("political_democracy_solution")

############################################################################################
### specification - RAMMatrices
############################################################################################

# w.o. meanstructure -----------------------------------------------------------------------

x = Symbol.("x".*string.(1:31))

S =[:x1   0    0     0     0      0     0     0     0     0     0     0     0     0
    0     :x2  0     0     0      0     0     0     0     0     0     0     0     0
    0     0     :x3  0     0      0     0     0     0     0     0     0     0     0
    0     0     0     :x4  0      0     0     :x15  0     0     0     0     0     0
    0     0     0     0     :x5   0     :x16  0     :x17  0     0     0     0     0
    0     0     0     0     0     :x6  0      0     0     :x18  0     0     0     0
    0     0     0     0     :x16  0     :x7   0     0     0     :x19  0     0     0
    0     0     0     :x15 0      0     0     :x8   0     0     0     0     0     0
    0     0     0     0     :x17  0     0     0     :x9   0     :x20  0     0     0
    0     0     0     0     0     :x18 0      0     0     :x10  0     0     0     0
    0     0     0     0     0     0     :x19  0     :x20  0     :x11  0     0     0
    0     0     0     0     0     0     0     0     0     0     0     :x12  0     0
    0     0     0     0     0     0     0     0     0     0     0     0     :x13  0
    0     0     0     0     0     0     0     0     0     0     0     0     0     :x14]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1.0   0     0
    0  0  0  0  0  0  0  0  0  0  0     :x21  0     0
    0  0  0  0  0  0  0  0  0  0  0     :x22  0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1.0   0
    0  0  0  0  0  0  0  0  0  0  0     0     :x23  0
    0  0  0  0  0  0  0  0  0  0  0     0     :x24  0
    0  0  0  0  0  0  0  0  0  0  0     0     :x25  0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     :x26
    0  0  0  0  0  0  0  0  0  0  0     0     0     :x27
    0  0  0  0  0  0  0  0  0  0  0     0     0     :x28
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     :x29  0     0
    0  0  0  0  0  0  0  0  0  0  0     :x30  :x31  0]

spec = RAMMatrices(;
    A = A,
    S = S,
    F = F,
    parameters = x,
    colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65]
)

partable = ParameterTable(spec)

# w. meanstructure -------------------------------------------------------------------------

x = Symbol.("x".*string.(1:38))

M = [:x32; :x33; :x34; :x35; :x36; :x37; :x38; :x35; :x36; :x37; :x38; 0.0; 0.0; 0.0]

spec_mean = RAMMatrices(;
    A = A,
    S = S,
    F = F,
    M = M,
    parameters = x,
    colnames = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65])

partable_mean = ParameterTable(spec_mean)

start_test = [fill(1.0, 11); fill(0.05, 3); fill(0.05, 6); fill(0.5, 8); fill(0.05, 3)]
start_test_mean = [fill(1.0, 11); fill(0.05, 3); fill(0.05, 6); fill(0.5, 8); fill(0.05, 3); fill(0.1, 7)]

semoptimizer = SemOptimizerOptim
@testset "RAMMatrices | constructor | Optim" begin include("constructor.jl") end
semoptimizer = SemOptimizerNLopt
@testset "RAMMatrices | constructor | NLopt" begin include("constructor.jl") end

if ENV["JULIA_EXTENDED_TESTS"] == "true"
    semoptimizer = SemOptimizerOptim
    @testset "RAMMatrices | parts | Optim" begin include("by_parts.jl") end
    semoptimizer = SemOptimizerNLopt
    @testset "RAMMatrices | parts | NLopt" begin include("by_parts.jl") end
end

@testset "constraints | NLopt" begin include("constraints.jl") end

############################################################################################
### specification - RAMMatrices → ParameterTable
############################################################################################

spec = ParameterTable(spec)
spec_mean = ParameterTable(spec_mean)

partable = spec
partable_mean = spec_mean

semoptimizer = SemOptimizerOptim
@testset "RAMMatrices → ParameterTable | constructor | Optim" begin include("constructor.jl") end
semoptimizer = SemOptimizerNLopt
@testset "RAMMatrices → ParameterTable | constructor | NLopt" begin include("constructor.jl") end

if ENV["JULIA_EXTENDED_TESTS"] == "true"
    semoptimizer = SemOptimizerOptim
    @testset "RAMMatrices → ParameterTable | parts | Optim" begin include("by_parts.jl") end
    semoptimizer = SemOptimizerNLopt
    @testset "RAMMatrices → ParameterTable | parts | NLopt" begin include("by_parts.jl") end
end

############################################################################################
### specification - Graph
############################################################################################

observed_vars = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8]
latent_vars = [:ind60, :dem60, :dem65]

graph = @StenoGraph begin
    # loadings
    ind60 → fixed(1)*x1 + x2 + x3
    dem60 → fixed(1)*y1 + y2 + y3 + y4
    dem65 → fixed(1)*y5 + y6 + y7 + y8
    # latent regressions
    label(:a)*dem60 ← ind60
    dem65 ← dem60
    dem65 ← ind60
    # variances
    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars) ↔ _(latent_vars)
    # covariances
    y1 ↔ y5
    y2 ↔ y4 + y6
    y3 ↔ y7
    y8 ↔ y4 + y6
end

spec = ParameterTable(
    latent_vars = latent_vars,
    observed_vars = observed_vars,
    graph = graph)

sort!(spec)

partable = spec

# meanstructure
mean_labels = label.([:m1, :m2, :m3, :m4, :m5, :m6, :m7, :m4, :m5, :m6, :m7])

graph = @StenoGraph begin
    # loadings
    ind60 → fixed(1)*x1 + x2 + x3
    dem60 → fixed(1)*y1 + y2 + y3 + y4
    dem65 → fixed(1)*y5 + y6 + y7 + y8
    # latent regressions
    label(:a)*dem60 ← ind60
    dem65 ← dem60
    dem65 ← ind60
    # variances
    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars) ↔ _(latent_vars)
    # covariances
    y1 ↔ y5
    y2 ↔ y4 + y6
    y3 ↔ y7
    y8 ↔ y4 + y6
    # means
    Symbol("1") → _(mean_labels).*_(observed_vars)
end

spec_mean = ParameterTable(
    latent_vars = latent_vars,
    observed_vars = observed_vars,
    graph = graph)

sort!(spec_mean)

partable_mean = spec_mean

start_test = [fill(0.5, 8); fill(0.05, 3); fill(1.0, 11);  fill(0.05, 9)]
start_test_mean = [fill(0.5, 8); fill(0.05, 3); fill(1.0, 11); fill(0.05, 3); fill(0.05, 13)]

semoptimizer = SemOptimizerOptim
@testset "Graph → ParameterTable | constructor | Optim" begin include("constructor.jl") end
semoptimizer = SemOptimizerNLopt
@testset "Graph → ParameterTable | constructor | NLopt" begin include("constructor.jl") end

if ENV["JULIA_EXTENDED_TESTS"] == "true"
    semoptimizer = SemOptimizerOptim
    @testset "Graph → ParameterTable | parts | Optim" begin include("by_parts.jl") end
    semoptimizer = SemOptimizerNLopt
    @testset "Graph → ParameterTable | parts | NLopt" begin include("by_parts.jl") end
end


text:
using PyCall
using JLD

gym = pyimport("gym")

torch = pyimport("torch")
nn    = pyimport("torch.nn")
f     = pyimport("torch.nn.functional")

# %%

include("../../Modules/Gridcells/StateSpikeGaussGrid.jl")
include("animation.jl")

# %% define the constant cells

env = gym.make("CartPole-v1")

cells_tmp_pos = HexModule[]
cells_tmp_vel = HexModule[]

maxs = env.observation_space.high
mins = env.observation_space.low

borpos = [convert(Float64,mins[1]),convert(Float64,maxs[1])]
borvpos = [-5.0,5.0]
#bortheta = [convert(Float64,mins[3]),convert(Float64,maxs[3])]
bortheta = [-pi,pi]
borvtheta = [-15.0,15.0]

for neuron in 1:16
    res = 0.07
    dilat = rand()*(1.4-0.7)+0.7
    theta = rand()*(2pi/6)-pi/6
    mod_pos = HexModule(res,dilat = dilat, theta = theta, bor_x = borpos, bor_y = bortheta)
    mod_vel = HexModule(res,dilat = dilat, theta = theta, bor_x = borvpos, bor_y = borvtheta)
    push!(cells_tmp_pos,mod_pos)
    push!(cells_tmp_vel,mod_vel)
end

const cells_pos = cells_tmp_pos
const cells_vel = cells_tmp_vel

env.close()

# %% network

@pydef mutable struct DQLnet <: nn.Module
    function __init__(self)
        pybuiltin(:super)(DQLnet, self).__init__()
        self.fc1 = nn.Linear(4,16)
        self.fc2 = nn.Linear(16,16)
        self.fc3 = nn.Linear(16,2)
    end

    function forward(self, s)
        s = f.relu(self.fc1(s))
        s = f.relu(self.fc2(s))
        s = self.fc3(s)
        return s
    end
end

Qnet = DQLnet()
Qnet.load_state_dict(torch.load("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/saved_cartpole.net"))

# %% functions

function convert_rgb(img::Array{UInt8,3})
    x,y,_ = size(img)
    img_rgb = fill(RGB(1.0,1.0,1.0),x,y)
    for i in 1:x
        for j in 1:y
            img_rgb[i,j] = RGB(img[i,j,1]/256,img[i,j,2]/256,img[i,j,3]/256)
        end
    end
    return img_rgb
end

function scatter_spikes(spikes::Array{Any,1})
    x = Int64[]
    y = Int64[]
    for i in 1:length(spikes)
        for j in eachindex(spikes[i])
            push!(x,i)
            push!(y,spikes[i][j])
        end
    end
    return x,y
end

function play_env(n_steps::Int64, render::Bool=false, random::Bool=false)
    env = gym.make("CartPole-v1")
    env.seed(0)
    obs = env.reset()
    done = false
    display_spikes = []
    display_img = []
    for j in 1:n_steps
        if random == false
            action = torch.argmax(Qnet(torch.tensor(obs).float())).item()
        else
            action = rand(0:1)
        end
        obs_new, reward, done, _ = env.step(action)
        if render==true
            # img = env.render(mode = "rgb_array")
            # push!(display_img, convert_rgb(img))
            env.render()
        end
        @show j
        @show obs_new
        view_spikes =input_spikes(obs_new,cells_pos,cells_vel)
        @show view_spikes
        push!(display_spikes, view_spikes)
        obs = obs_new
    end

    x_spikes, y_spikes = scatter_spikes(display_spikes)
    env.close()
    env = nothing
    save("Julia/CartPole/Heatmaps/disp_spikesGaussGrid.jld", "display_spikes", display_spikes)
    # scatter(x_spikes,y_spikes,
    #         title ="HexGrid input viz",
    #         markersize = 2,
    #         xlabel = "Time step",
    #         ylabel = "Neuron",
    #         legend = false)
    # plot!([64], seriestype = "hline", color = :black)
    #savefig("input_gaussgrid_rand.png")
    # plot()
    # anim = @animate for i in 1:n_steps
    #     anim_viz(x_spikes,y_spikes,display_img[i],160,i)
    # end
    # gif(anim, "Julia/Cartpole/Input_Visual/anim_grid.gif", fps = 20)
end

# %% run visualization

play_env(200,true,true)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
mutable struct Prop2DAcoIsoDenQ_DEO2_FDTD
    p::Ptr{Cvoid}
end

function Prop2DAcoIsoDenQ_DEO2_FDTD(;
        nz=0,
        nx=0,
        nsponge=60,
        nbz=512,
        nbx=8,
        dz=5.0,
        dx=5.0,
        dt=0.001,
        nthreads=Sys.CPU_THREADS,
        freesurface=true,
        freqQ=5.0,
        qMin=0.1,
        qInterior=100.0)
    nz,nx,nsponge,nbz,nbx,nthreads = map(x->round(Int,x), (nz,nx,nsponge,nbz,nbx,nthreads))
    dz,dx,dt = map(x->Float32(x), (dz,dx,dt))

    @assert nx > 0
    @assert nz > 0
    @assert nsponge > 0
    @assert nthreads > 0
    @assert nbx > 0
    @assert nbz > 0

    fs = freesurface ? 1 : 0

    p = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_alloc, libprop2DAcoIsoDenQ_DEO2_FDTD),
        Ptr{Cvoid},
        (Clong, Clong,    Clong, Clong, Clong,   Cfloat, Cfloat, Cfloat, Clong, Clong),
         fs,    nthreads, nx,    nz,    nsponge, dx,     dz,     dt,     nbx,   nbz)

    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_SetupDtOmegaInvQ, libprop2DAcoIsoDenQ_DEO2_FDTD),
        Cvoid,
        (Ptr{Cvoid}, Cfloat, Cfloat, Cfloat),
         p,          freqQ,  qMin,   qInterior)

    Prop2DAcoIsoDenQ_DEO2_FDTD(p)
end

free(prop::Prop2DAcoIsoDenQ_DEO2_FDTD) = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_free, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

function size(prop::Prop2DAcoIsoDenQ_DEO2_FDTD)
    nz = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_getNz, libprop2DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nx = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_getNx, libprop2DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    (nz,nx)
end

for _f in (:V, :B, :PSpace, :PCur, :POld, :DtOmegaInvQ)
    symf = "Prop2DAcoIsoDenQ_DEO2_FDTD_get" * string(_f)
    @eval $(_f)(prop::Prop2DAcoIsoDenQ_DEO2_FDTD) = unsafe_wrap(Array, ccall(($symf, libprop2DAcoIsoDenQ_DEO2_FDTD), Ptr{Float32}, (Ptr{Cvoid},), prop.p), size(prop), own=false)
end

abstract type Prop2DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V <: Prop2DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B <: Prop2DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB <: Prop2DAcoIsoDenQ_DEO2_FDTD_Model end

propagateforward!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD) = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_TimeStep, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)
propagateadjoint!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD) = propagateforward!(prop) # self-adjoint

scale_spatial_derivatives!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD) =
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_ScaleSpatialDerivatives, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

# v in model-space
function forwardBornInjection!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_V, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"])
end

# v,b in model-space
function forwardBornInjection!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_VB, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},     Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# b in model-space
function forwardBornInjection!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_B, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},     Ptr{Cfloat}),
         prop.p,     dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# v in model-space
function adjointBornAccumulation!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_V, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"])
end

# v,b in model-space
function adjointBornAccumulation!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_VB, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},     Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# b in model-space
function adjointBornAccumulation!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_B, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},     Ptr{Cfloat}),
         prop.p,     dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# v in model-space with FWI wavefield separation
function adjointBornAccumulation!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionWaveFieldSeparationFWI, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},      Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], 1)
 end

# v in model-space with RTM wavefield separation
function adjointBornAccumulation!(prop::Prop2DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionWaveFieldSeparationRTM, dmodel, wavefields)
    ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep, libprop2DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},      Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], 0)
 end

function show(io::IO, prop::Prop2DAcoIsoDenQ_DEO2_FDTD)
    nx = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_getNx, libprop2DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nz = ccall((:Prop2DAcoIsoDenQ_DEO2_FDTD_getNz, libprop2DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    write(io, "Prop2DAcoIsoDenQ_DEO2_FDTD -- nx,nz; $nx,$nz")
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
mutable struct SGD{LR <: Real, GC <: Real}
    lr::LR
    gclip::GC
end

function update!(w, g, p::SGD)
end


text:
using jInv.Mesh;
using KrylovMethods
using Multigrid.DomainDecomposition

import jInv.LinearSolvers.copySolver;
import jInv.LinearSolvers.AbstractSolver;
import jInv.Utils.clear!
import Base.copy

export MGparam;
export getMGparam, MGsetup, clear!
export BlockFGMRES,hierarchyExists,copySolver,destroyCoarsestLU
export multilevelOperatorConstructor, getMultilevelOperatorConstructor

include("SpMatMul.jl");
include("FGMRES.jl");

"""
mutable struct Multigrid.multilevelOperatorConstructor

Fields:

	param - parameters for the PDE
	getOperator(mesh,param) - a function that creates the PDE from
	restrictParams(mesh_fine,mesh_coarse,param_fine,level) - returns param_coarse for level+1.

"""
mutable struct multilevelOperatorConstructor
	param
	getOperator		::Function
	restrictParams  ::Function
end

function getMultilevelOperatorConstructor(param,getOperator::Function,restrictParams)
if restrictParams == []
	## PDE with no params
	restrictParams2(mesh_fine,mesh_coarse,param_fine,level) = [];
	getOperator2(Mesh,param) = getOperator(Mesh);
	return multilevelOperatorConstructor([],getOperator2,restrictParams2);
else
	return multilevelOperatorConstructor(param,getOperator,restrictParams);
end
end
"""
mutable struct Multigrid.CYCLEmem

Fields:

	b::Array - memory for the right-hand-side
	r::Array - memory for the residual
	x::Array - memory for the iterated solution
"""
mutable struct CYCLEmem{VAL}
	b 					::Array{VAL}
	r					::Array{VAL}
	x					::Array{VAL}
end

"""
mutable struct Multigrid.MGparam

Fields:

	levels::Int64            - Maximum number of multigrid levels
	numCores::Int64          - Number of OMP cores to work with. Some operations (setup) are not parallelized.
	maxOuterIter::Int64      - Maximum outer iterations.
	relativeTol::Float64	 - Relative L2/Frobenius norm stopping criterion.
	relaxType:: String	 - Relax type. Can be "Jac", "Jac-GMRES" or "SPAI".
	relaxParam::Float64	     - Relax damping parameter.
	relaxPre::Function 	     - pre and post relaxation numbers
	relaxPost::Function	     - Can be 'V', 'F', 'W', 'K' (Krylov cycles are done with FGMRES).
	Ps::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	Rs::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	As::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	relaxPrecs -  an array of relaxation preconditioners for all levels.
	memCycle::Array{CYCLEmem} - Space for x,b and r for each level.
	memRelax::Union{Array{FGMRESmem},Array{BlockFGMRESmem}}  - This is used just in case of GMRES relaxation.
	memKcycle::Union{Array{FGMRESmem},Array{BlockFGMRESmem}} - Memory for the Krylov-cycle FGMRES. First field is ignored.
	coarseSolveType::String - Can be "MUMPS" or "NoMUMPS" for Julia backslash.
	LU - Factorization of coarsest level.
	doTranspose::Int64
	strongConnParam::Float64  - (for SA-AMG only) A threshold for determining a strong connection should >0.25, and <0.85.
	FilteringParam::Float64	  - (for SA-AMG only) A threshold for prolongation filtering >0.0, and <0.2.
	Meshes::Array{RegularMesh} 		 - Array of Regular Meshes for geometric multigrid.
	transferOperatorType:: String	- (for geometric MG only) May be "FullWeighting", "SystemsFacesLinear"
	singlePrecision		:: Bool - indicator for single precision computation. For internal use.
"""
mutable struct MGparam{VAL,IND}
	levels				:: Int64
	numCores			:: Int64
	maxOuterIter		:: Int64
	relativeTol			:: Float64
	relaxType		    :: String
	relaxParam			:: Union{Array{Float64},Float64,Any}
	relaxPre			:: Function
	relaxPost			:: Function
	cycleType			:: Char
	Ps					:: Array{SparseCSCTypes}
	Rs					:: Array{SparseCSCTypes}
	As					:: Array{SparseMatrixCSC{VAL,IND}}
	relaxPrecs
	memCycle			:: Array{CYCLEmem}
	memRelax			:: Array{FGMRESmem}
	memKcycle			:: Array{FGMRESmem}
	coarseSolveType		:: String
	LU
	doTranspose			:: Int64
	strongConnParam		:: Float64
	FilteringParam		:: Float64
	Meshes				:: Union{Array{RegularMesh},RegularMesh}
	transferOperatorType:: String
	singlePrecision		:: Bool
end

const spValType = ComplexF32

include("GeometricTransferOperators.jl")
include("MGsetup.jl");
include("SA-AMG.jl");
include("MGcycle.jl");
include("SolveFuncs.jl");
include("SAAMGWrapper.jl");
include("MGWrapper.jl");
include("Systems.jl");
include("parRelax.jl");
include("Vanka.jl");
include("ClassicalAMG.jl")
"""
function Multigrid.copySolver(MG::MGparam)

copies the solver parameters without the setup and allocated memory.
"""
function copySolver(MG::MGparam{VAL,IND}) where {VAL,IND}
return getMGparam(VAL,IND,MG.levels,MG.numCores,MG.maxOuterIter,MG.relativeTol,MG.relaxType,MG.relaxParam,
					MG.relaxPre,MG.relaxPost,MG.cycleType,MG.coarseSolveType,MG.strongConnParam,MG.FilteringParam,MG.transferOperatorType);
end

function getMGparam(VAL::Type,IND::Type,levels::Int64,numCores::Int64,maxIter::Int64,relativeTol:: Float64,relaxType::String,relaxParam,
					relaxPre::Function,relaxPost::Function,cycleType::Char='V',coarseSolveType::String="NoMUMPS",strongConnParam::Float64=0.4,FilteringParam::Float64 = 0.0,transferOperatorType = "FullWeighting")
singlePrecision = (VAL==Float32 || VAL==ComplexF32);
return MGparam{VAL,IND}(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,[],[],[],[],Array{CYCLEmem}(undef,0),
				Array{FGMRESmem{VAL}}(undef,0),Array{FGMRESmem{VAL}}(undef,0),coarseSolveType,[],0,strongConnParam,FilteringParam,Array{RegularMesh}(undef,0),transferOperatorType,singlePrecision);
end

function getMGparam(VAL::Type,IND::Type,levels::Int64=3,numCores::Int64=8,maxIter::Int64=20,relativeTol::Float64=1e-6,relaxType::String="SPAI",relaxParam=1.0,
					relaxPre::Int64=2,relaxPost::Int64=2,cycleType::Char='V',coarseSolveType::String="NoMUMPS",strongConnParam::Float64=0.4,FilteringParam::Float64 = 0.0,transferOperatorType = "FullWeighting")
relaxPreFun(x) = relaxPre;
relaxPostFun(x) = relaxPost;
return getMGparam(VAL,IND,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPreFun,relaxPostFun,cycleType,coarseSolveType,strongConnParam,FilteringParam,transferOperatorType);
end

function getCYCLEmem(n::Int64,m::Int64,T::Type,withB::Bool=true)
b = zeros(T,0);
if m==1
	if withB
		b = zeros(T,n);
	end
	return CYCLEmem(b,zeros(T,n),zeros(T,n));
else
	if withB
		b = zeros(T,n,m);
	end
	return CYCLEmem(b,zeros(T,n,m),zeros(T,n,m));
end
end

# import jInv.Utils.clear!
function clear!(param::MGparam{VAL,IND}) where {VAL,IND}
param.Ps = Array{SparseMatrixCSC{real(VAL),IND}}(undef,0);
param.Rs = Array{SparseMatrixCSC{real(VAL),IND}}(undef,0);
param.As = Array{SparseMatrixCSC{VAL,IND}}(undef,0);
param.relaxPrecs = [];
param.memCycle = Array{CYCLEmem{VAL}}(undef,0);
param.memRelax = Array{FGMRESmem{VAL}}(undef,0);
param.memKcycle = Array{FGMRESmem{VAL}}(undef,0);
param.Meshes = Array{RegularMesh}(undef,0);
destroyCoarsestLU(param);
end

function destroyCoarsestLU(param::MGparam{VAL,IND}) where {VAL,IND}
if param.LU==[]
	return;
end
if param.coarseSolveType=="MUMPS"
	if isa(param.LU,MUMPSfactorization)
		destroyMUMPS(param.LU);
		param.LU = [];
	end
elseif isa(param.LU,AbstractSolver)
	clear!(param.LU);
else
	param.LU = [];
end
return;
end

function hierarchyExists(param::MGparam{VAL,IND}) where {VAL,IND}
return length(param.As) > 0;
end


text:
using PyPlot

function plot_matrix(M, posneg=true,eps=1e-5)

    nrow = size(M,1)
    ncol = size(M,2)

    ntot = nrow * ncol

    A = reshape(M, ntot)

    irow = zeros(Int,ntot)
    icol = zeros(Int,ntot)

    er = maxabs(A)*eps

    cnt = 1
    for j = 1:ncol
        for i = 1:nrow
            irow[cnt] = i
            icol[cnt] = j
            cnt = cnt + 1
        end
    end

    plot([0, ncol+1, ncol+1, 0, 0], [0, 0, nrow+1, nrow+1, 0], color="black", linewidth=3)
    for i = 1:ncol
        axvline(i, color="lightgray")
    end
    for i = 1:nrow
        axhline(i, color="lightgray")
    end

    if posneg
        ipos = A .> er
        ineg = A .< -er

        plot(icol[ipos], nrow + 1 - irow[ipos], "ro")
        plot(icol[ineg], nrow + 1 - irow[ineg], "bo")
    else
        idx = abs(A) .> er
        plot(icol[idx], nrow + 1 - irow[idx], "o", color="black")
    end

    return
end


text:
function print_all_descriptions(root::String = ".")
    fs = joinpath.(root, readdir(root));
    for f in fs
        if isdir(f)
            try
                print_description(f);
            catch ignored
            end
        end
    end
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager