text:
################################################################################################################################
# video. Video Analysis
################################################################################################################################

# Motion Analysis and Object Tracking

#------------------------------------------------------------------------------------------------------------------------------#
# calcOpticalFlowPyrLK
# Calculates an optical flow for a sparse feature set using the iterative Lucas-Kanade method with pyramids
# Lucas-Kanade optical flow in pyramids (Bouguet 2000)

# void calcOpticalFlowPyrLK(InputArray prevImg, InputArray nextImg, InputArray prevPts, InputOutputArray nextPts,
# OutputArray status, OutputArray err, Size winSize=Size(21,21), int maxLevel=3,
# TermCriteria criteria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01), int flags=0, double minEigThreshold=1e-4 )

# parameters:
# prevImg   – first 8-bit input image (InputArray) or pyramid constructed by buildOpticalFlowPyramid()
# nextImg   – second input image (InputArray) or pyramid of the same size and the same type as prevImg
# prevPts   – vector of 2D points (InputArray single-precision fp X,Y) to find flow
# nextPts   – output vector of 2D points (InputOutputArray single-precision fp X,Y) with new positions of input features in the second image
#             if OPTFLOW_USE_INITIAL_FLOW => Set size of vector = size of input vector
# status    – output status vector (uchar)
#             status = 1 if flow for the corresponding features has been found, else status = 0
# err       – output vector of errors; each element of the vector is set to an error for the corresponding feature
#             type of the error measure can be set in flags parameter
#             if no flow, error is not defined (status = 0)
# winSize   – cv::Size (size) of the search window at each pyramid level, default = cv::Size(21,21)
# maxLevel  – int 0-based maximal pyramid level number, default = 3
#             if maxLevel = 0, single level, i.e., pyramids are not used
#             if maxLevel = 1, two levels; etc
#             pyramids < maxLevel
# criteria  – TermCriteria parameter of the iterative search algorithm
#             stop when max iterations criteria.maxCount || search window moves < criteria.epsilon
#             default = TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01)

# flags, int default =0
# OPTFLOW_USE_INITIAL_FLOW       =>  uses initial estimations, stored in nextPts
#                                    default: prevPts => nextPts
# OPTFLOW_LK_GET_MIN_EIGENVALS   =>  use minimum eigen values as an error measure (minEigThreshold);
#                                    default: L1 distance between patches around the original and a moved point,
#                                    divided by number of pixels in a window, is used as a error measure
# minEigThreshold (double)       =>  min eigen value of a 2x2 normal matrix of optical flow equations (spatial gradient matrix in [Bouguet00]/window pixels)
#                                    if eigen value < minEigThreshold, feature is filtered out and flow ignored
#                                    default = 1e-4
#-------------------------------------------------------------------------------------------------------------------#

# e.g., parameters
# prevImg::Mat(rows, columns, CV_8UC3)  || buildOpticalFlowPyramid()
# nextImg::Mat(rows, columns, CV_8UC3)
# prevPts = tostdvec([pts1::cvPoint2f, pts2::cvPoint2f, pts3::cvPoint2f, pts4::cvPoint2f, pts5::cvPoint2f . . .])
# nextPts = tostdvec([pts1::cvPoint2f, pts2::cvPoint2f, pts3::cvPoint2f, pts4::cvPoint2f, pts5::cvPoint2f . . .])

calcOpticalFlowPyrLK(prevImg, nextImg, prevPts, nextPts, status, err, winSize=Size(21,21), maxLevel=3,
    criteria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01), flags=0, minEigThreshold=1e-4) =
       @cxx cv::calcOpticalFlowPyrLK(prevImg, nextImg, prevPts, nextPts, status, err, winSize, maxLevel, criteria, flags, minEigThreshold)

# buildOpticalFlowPyramid
# Constructs the image pyramid which can be passed to calcOpticalFlowPyrLK().
# int buildOpticalFlowPyramid(InputArray img, OutputArrayOfArrays pyramid, Size winSize, int maxLevel, bool withDerivatives=true,
# int pyrBorder=BORDER_REFLECT_101, int derivBorder=BORDER_CONSTANT, bool tryReuseInputImage=true)
# Parameters:
# img                – 8-bit input image
# pyramid            – output pyramid
# winSize            – window size of optical flow algorithm, >= winSize in calcOpticalFlowPyrLK()
# maxLevel           – 0-based maximal pyramid level number
# withDerivatives    – set to precompute gradients for the every pyramid level
# pyrBorder          – the border mode for pyramid layers
# derivBorder        – the border mode for gradients
# tryReuseInputImage – put ROI of input image into the pyramid if possible
#                      if tryReuseInputImage = false, then forces data copying
# Output:
# number of levels in constructed pyramid, maybe < maxLevel

buildOpticalFlowPyramid(img, pyramid, winSize, maxLevel, withDerivatives=true, pyrBorder=BORDER_REFLECT_101,
    derivBorder=BORDER_CONSTANT, tryReuseInputImage=true) = @cxx cv::buildOpticalFlowPyramid(img, pyramid, winSize,
        maxLevel, withDerivatives, pyrBorder, derivBorder, tryReuseInputImage)

# calcOpticalFlowFarneback
# Computes a dense optical flow using the Gunnar Farneback’s algorithm
# void calcOpticalFlowFarneback(InputArray prev, InputArray next, InputOutputArray flow, double pyr_scale,
# int levels, int winsize, int iterations, int poly_n, double poly_sigma, int flags)

# Parameters:
# prev       – first 8-bit single-channel input image
# next       – second input image: size and type = prev
# flow       – CV_32FC2 computed flow image, size = prev
# pyr_scale  – image scale (<1) to build pyramids for each image
#              pyr_scale=0.5 => classical pyramid, next layer is 2X smaller than the previous one
# levels     – number of pyramid layers including the initial image
#              levels=1 => only the original images are used
# winsize    – averaging window size
#              larger => increase robustness, less noise and better fast motion detection), but yield more blurred motion field
# iterations – iterations/pyramid level
# poly_n     – typically poly_n =5 or 7
#              size of the pixel neighborhood used to find polynomial expansion in each pixel
#              larger => smoother surfaces, robust but more blurred motion field,
# poly_sigma – standard deviation of the Gaussian to smooth derivatives for the polynomial expansion
#              poly_n=5  => poly_sigma=1.1
#              poly_n=7  => poly_sigma=1.5
# flags      – same as calcOpticalFlowPyrLK

calcOpticalFlowFarneback(prev, next, flow, pyr_scale=0.5, levels=1, winsize=10, iterations=5, poly_n=5,
    poly_sigma=1.1, flags=0) = @cxx cv::calcOpticalFlowFarneback(prev, next, flow, pyr_scale, levels, winsize,
      iterations, poly_n, poly_sigma, flags)

# estimateRigidTransform
# Computes an optimal affine transformation between two 2D point sets
# Mat estimateRigidTransform(InputArray src, InputArray dst, bool fullAffine)

# Parameters:
# src – 2D point set stored in std::vector or Mat, or an image stored in Mat
# dst – 2D point set of the same size and the same type as A, or another image
# fullAffine =true  => affine transformation with no additional restrictions (6 dof) - degrees of freedom
# fullAffine =false => translation, rotation, and uniform scaling (5 dof)

estimateRigidTransform(src, dst, fullAffine=false) = @cxx cv::estimateRigidTransform(src, dst, fullAffine)

# findTransformECC
# Finds the geometric transform (warp) between two images in terms of the ECC criterion [EP08]
# double findTransformECC(InputArray templateImage, InputArray inputImage, InputOutputArray warpMatrix,
# int motionType=MOTION_AFFINE, TermCriteria criteria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 50, 0.001))

# Parameters:
# templateImage – single-channel template image, CV_8U or CV_32F array
# inputImage    – single-channel input image which should be warped with the final warpMatrix, type = temlateImage
# warpMatrix    – floating-point 2x3 or 3x3 mapping matrix (warp)
# motionType    – type of motion parameter
# MOTION_TRANSLATION  => translational motion model: warpMatrix is 2x3 with the first 2x2 part being the unity matrix and the rest two parameters being estimated
# MOTION_EUCLIDEAN    => Euclidean (rigid) transformation as motion model: 2 parameters are estimated, warpMatrix is 2x3
# MOTION_AFFINE       => affine motion model (DEFAULT):                    6 parameters are estimated, warpMatrix is 2x3
# MOTION_HOMOGRAPHY   => homography as a motion model:                     8 parameters are estimated, warpMatrix is 3x3
# criteria      – termination criteria of the ECC algorithm
#                 criteria.epsilon = threshold of the increment in the correlation coefficient between two iterations
#                 criteria.epsilon < 0 => criteria.maxcount is termination criterion

findTransformECC(templateImage, inputImage, warpMatrix, motionType=MOTION_AFFINE, criteria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 50, 0.001)) =
      @cxx cv::findTransformECC(templateImage, inputImage, warpMatrix, motionType, criteria)

#-------------------------------------------------------------------------------------------------------------------#
# CamShift
# Finds an object center, size, and orientation
# CAMSHIFT object tracking algorithm [Bradski98]
# RotatedRect CamShift(InputArray probImage, Rect& window, TermCriteria criteria)

# Parameters:
# probImage – Back projection of the object histogram. See calcBackProject() .
# window    – Initial search window.
# criteria  – Stop criteria for the underlying meanShift() .
# Returns:
# RotatedRect

# e.g., input
# probImage = calcBackProject(...) see OpenCV_imgproc.jl
# window = cvRect
# criteria = TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01)
CamShift(probImage, window, criteria) = @cxx cv::CamShift(probImage, window, criteria)

# meanShift
# Finds an object on a back projection image
# iterative object search algorithm
# int meanShift(InputArray probImage, Rect& window, TermCriteria criteria)

# Parameters:
# probImage  – Back projection of the object histogram. See calcBackProject().
# window     – Initial search window.
# criteria   – Stop criteria for the iterative search algorithm
# Returns:
# Number of iterations CAMSHIFT took to converge

meanShift(probImage, window, criteria) = @cxx cv::meanShift(probImage, window, criteria)

#-------------------------------------------------------------------------------------------------------------------#
# KalmanFilter
# class
# implements a standard Kalman filter http://en.wikipedia.org/wiki/Kalman_filter, [Welch95]
# however, can be extended by modifying transitionMatrix, controlMatrix, and measurementMatrix

# Parameters:
# dynamParams   – Dimensionality of the state
# measureParams – Dimensionality of the measurement
# controlParams – Dimensionality of the control vector
# type          – Type of the created matrices that should be CV_32F or CV_64F

# constructors, return KF
KalmanFilter() = @cxxnew cv::KalmanFilter::KalmanFilter()
KalmanFilter(dynamParams::Int,measureParams::Int, controlParams=0, cvtype=CV_32F) =
     @cxxnew cv::KalmanFilter::KalmanFilter(dynamParams,measureParams, controlParams, cvtype)

# Re-initializes Kalman filter
KalmanInit(KF, dynamParams::Int,measureParams::Int, controlParams=0, cvtype=CV_32F) =
      @cxx KF->init(dynamParams,measureParams, controlParams, cvtype)

# Compute a predicted state (optional= cv::Mat)
KalmanPredict(KF, control=Mat()) = @cxx KF->predict(control)

# Update the predicted state from the measurement (measurement = cv::Mat)
KalmanCorrect(KF,measurement) = @cxx KF->correct(measurement)
#-------------------------------------------------------------------------------------------------------------------#

# BackgroundSubtractor class
# BackgroundSubtractor::apply
# Computes a foreground mask
# void BackgroundSubtractor::apply(InputArray image, OutputArray fgmask, double learningRate=-1)
# Parameters:
# image        – Next video frame
# fgmask       – The output foreground mask as an 8-bit binary image
# learningRate – 0-1 that indicates how fast the background model is learnt
#                  < 0 automatically chosen learning rate
#                    0 not updated at all
#                    1 reinitialized from the last frame

# Create a pointer to the abstract class BackgroundSubtractorMOG
cxx"""
   cv::bgsegm::BackgroundSubtractorMOG* BackgroundSubtractorMOG()
   {
       cv::bgsegm::BackgroundSubtractorMOG* bkgsub;
       return bkgsub;
   }
"""

# Constructor
bkgsubMOG() = @cxx BackgroundSubtractorMOG()

# Apply
bkgsubApply(bkgsubMOG, Iimage, fgmask, learningRate=-1.0) = @cxx bkgsubMOG->apply(image,fgmask, learningRate)

# Computes a background image
# input: backgroundImage = Mat(r,c,cvtype)
getBackgroundImage(bkgsubMOG, backgroundImage) = @cxx bkgsubMOG->getBackgroundImage(backgroundImage)

# BackgroundSubtractorMOG class: algorithm described in [KB2001]
# Create a mixture-of-gaussian background subtractor
# Ptr<BackgroundSubtractorMOG> createBackgroundSubtractorMOG(int history=200, int nmixtures=5, double backgroundRatio=0.7, double noiseSigma=0)
# Parameters:
# history           – length of the history
# nmixtures         – Number of Gaussian mixtures
# backgroundRatio   – Background ratio
# noiseSigma        – Noise strength (SD for brightness, each color channel). 0 = automatic value

createBackgroundSubtractorMOG(history=200, nmixtures=5, backgroundRatio=0.7, noiseSigma=0) =
     @cxx cv::bgsegm::createBackgroundSubtractorMOG(history, nmixtures, backgroundRatio, noiseSigma)

# BackgroundSubtractorMOG2
# Creates MOG2 Background Subtractor
# Ptr<BackgroundSubtractorMOG2> createBackgroundSubtractorMOG2(int history=500, double varThreshold=16, bool detectShadows=true )
# Parameters:
# history       – length of the history
# varThreshold  – Threshold on the squared Mahalanobis distance between the pixel
# detectShadows – if =true, detect shadows and mark them. slows down!, so set to false.

# http://www.developark.com/5098_19238765/

# Create a pointer to the abstract class BackgroundSubtractorMOG2
# cxx"""
# cv::BackgroundSubtractorMOG2* BackgroundSubtractorMOG2()
# {
#      cv::BackgroundSubtractorMOG2* bkgsub;
#      return bkgsub;
# }
# """

createBackgroundSubtractorMOG2(history=500, varThreshold=16, detectShadows=true) =
    @cxx cv::createBackgroundSubtractorMOG(history, varThreshold, detectShadows)

# Returns the number of last frames that affect the background model
getHistory(bkgsubMOG2) = @cxx bkgsubMOG2->getHistory()

# Sets the number of last frames that affect the background model
setHistory(bkgsubMOG2, history::Int) = @cxx bkgsubMOG2->setHistory(history)

# Returns number of gaussian components in the background model
getNMixtures(bkgsubMOG2) = @cxx bkgsubMOG2->getNMixtures()

# Sets number of gaussian components in the background model, needs to be reinitialized
setNMixtures(bkgsubMOG2, nmixtures::Int) = @cxx bkgsubMOG2->setNMixtures(nmixtures)

# Returns the “background ratio” parameter of the algorithm
getBackgroundRatio(bkgsubMOG2) = @cxx bkgsubMOG2->getBackgroundRatio()

# Sets the “background ratio” parameter of the algorithm
setBackgroundRatio(bkgsubMOG2, ratio::Float64) = @cxx bkgsubMOG2->setBackgroundRatio(ratio)

# Returns the variance threshold for the pixel-model match
getVarThreshold(bkgsubMOG2) = @cxx bkgsubMOG2->getVarThreshold()

# Sets the variance threshold for the pixel-model match
setVarThreshold(bkgsubMOG2, varThreshold::Float64) = @cxx bkgsubMOG2-> setVarThreshold(varThreshold)

# Returns the variance threshold for the pixel-model match used for new mixture component generation
getVarThresholdGen(bkgsubMOG2) = @cxx bkgsubMOG2-> getVarThresholdGen()

# Set the variance threshold for the pixel-model match used for new mixture component generation
setVarThresholdGen(bkgsubMOG2, VarThreshold::Float64) = @cxx bkgsubMOG2-> setVarThresholdGen(VarThreshold)

# Returns the initial variance of each gaussian component
getVarInit(bkgsubMOG2) = @cxx bkgsubMOG2->getVarInit()

# Sets the initial variance of each gaussian component
setVarInit(bkgsubMOG2, varInit::Float64) = @cxx bkgsubMOG2->setVarInit(varInit)

# Returns the complexity reduction threshold
getComplexityReductionThreshold(bkgsubMOG2) = @cxx bkgsubMOG2->getComplexityReductionThreshold()

# Sets the complexity reduction threshold
setComplexityReductionThreshold(bkgsubMOG2, ct::Float64) = @cxx bkgsubMOG2->setComplexityReductionThreshold(ct)

# Returns the shadow detection flag
getDetectShadows(bkgsubMOG2) = @cxx bkgsubMOG2->getDetectShadows()

# Enables or disables shadow detection
setDetectShadows(bkgsubMOG2, detectShadows::Bool) = @cxx bkgsubMOG2->setDetectShadows(detectShadows)

# Returns the shadow value, Default value is 127
getShadowValue(bkgsubMOG2) = @cxx bkgsubMOG2->getShadowValue()

# Sets the shadow value
setShadowValue(bkgsubMOG2, value::Int) = @cxx bkgsubMOG2->setShadowValue(bkgsubMOG2, value)

# Returns the shadow threshold
getShadowThreshold(bkgsubMOG2) = @cxx bkgsubMOG2->getShadowThreshold()

# Sets the shadow threshold
setShadowThreshold(bkgsubMOG2, threshold::Float64) = @cxx bkgsubMOG2->setShadowThreshold(threshold)

# BackgroundSubtractorKNN
# K-nearest neigbours - based Background/Foreground Segmentation Algorithm
# The class implements the K-nearest neigbours background subtraction described in [Zivkovic2006]
# Very efficient if number of foreground pixels is low.

# C++: Ptr<BackgroundSubtractorKNN> createBackgroundSubtractorKNN(int history=500,
#     double dist2Threshold=400.0, bool detectShadows=true )

# Parameters:
# history        – Length of the history
# dist2Threshold – Threshold on the squared distance between the pixel
#                  This parameter does not affect the background update
# detectShadows  – If true, the algorithm will detect shadows and mark them
#                  set to false to speed up

createBackgroundSubtractorKNN(history=500, dist2Threshold=400.0,detectShadows=true) =
    @cxx cv::createBackgroundSubtractorKNN(history, dist2Threshold, detectShadows)

# Returns the number of last frames that affect the background model
getHistory() = @cxx bkgsubKNN->getHistory()

# Sets the number of last frames that affect the background model
setHistory(history::Int) = @cxx bkgsubKNN->setHistory(history)

# Returns the number of data samples in the background model
getNSamples() = @cxx bkgsubKNN->getNSamples()

# Sets the number of data samples in the background model. The model needs to be reinitalized to reserve memory
setNSamples(_nN::Int) = @cxx bkgsubKNN->setNSamples(_nN)

# Returns the threshold on the squared distance between the pixel and the sample
getDist2Threshold() = @cxx bkgsubKNN->getDist2Threshold()

# Sets the threshold on the squared distance
setDist2Threshold(_dist2Threshold::Float64) = @cxx bkgsubKNN->setDist2Threshold(_dist2Threshold)

# Returns the number of neighbours, the k in the kNN. K is the number of samples that need to be within
# dist2Threshold in order to decide that that pixel is matching the kNN background mode
getkNNSamples() = @cxx bkgsubKNN->getkNNSamples()

# Sets the k in the kNN. How many nearest neigbours need to match
setkNNSamples(_nkNN::Int) = @cxx bkgsubKNN->setkNNSamples(_nkNN)

# Returns the shadow detection flag
getDetectShadows() = @cxx bkgsubKNN->getDetectShadows()

# Enables or disables shadow detection
setDetectShadows(detectShadows::Bool) = @cxx bkgsubKNN->setDetectShadows(detectShadows)

# Returns the shadow value
getShadowValue() = @cxx bkgsubKNN->getShadowValue()

# Sets the shadow value
setShadowValue(value::Int) = @cxx bkgsubKNN->setShadowValue(value)

# Returns the shadow threshold
getShadowThreshold() = @cxx bkgsubKNN->getShadowThreshold()

# Sets the shadow threshold
setShadowThreshold(threshold::Float64) = @cxx bkgsubKNN->setShadowThreshold(threshold)

# BackgroundSubtractorGMG
# Background Subtractor module based on the algorithm given in [Gold2012]
# Creates a GMG Background Subtractor
# Ptr<BackgroundSubtractorGMG> createBackgroundSubtractorGMG(int initializationFrames=120,
# double decisionThreshold=0.8)

# Parameters:
# initializationFrames – number of frames used to initialize the background models
# decisionThreshold    – Threshold value, above which it is marked foreground, else background

createBackgroundSubtractorGMG(initializationFrames=120, decisionThreshold=0.8) =
   @cxx cv::bgsegm::createBackgroundSubtractorGMG(initializationFrames, decisionThreshold)

# Returns the number of frames used to initialize background model
getNumFrames(bkgsubGMG) = @cxx bkgsubGMG->getNumFrames()

# Sets the number of frames used to initialize background model
setNumFrames(nframes::Int) = @cxx bkgsubGMG->setNumFrames(nframes)

# Returns the learning rate of the algorithm
# It lies between 0.0 and 1.0. It determines how quickly features are “forgotten” from histograms
getDefaultLearningRate() = @cxx bkgsubGMG->getDefaultLearningRate()

# Sets the learning rate of the algorithm
setDefaultLearningRate(lr::Float64) = @cxx bkgsubGMG->setDefaultLearningRate(lr)

# Returns the value of decision threshold
# Decision value is the value above which pixel is determined to be FG
getDecisionThreshold() = @cxx bkgsubGMG->getDecisionThreshold()

# Sets the value of decision threshold
setDecisionThreshold(thresh::Float64) = @cxx bkgsubGMG->setDecisionThreshold(thresh)

# Returns total number of distinct colors to maintain in histogram
getMaxFeatures() = @cxx bkgsubGMG->getMaxFeatures()

# Sets total number of distinct colors to maintain in histogram
setMaxFeatures(maxFeatures::Int) = @cxx bkgsubGMG->setMaxFeatures(maxFeatures)

# Returns the parameter used for quantization of color-space.
# It is the number of discrete levels in each channel to be used in histograms.
getQuantizationLevels() = @cxx bkgsubGMG->getQuantizationLevels()

# Sets the parameter used for quantization of color-space
setQuantizationLevels(nlevels::Int) = @cxx bkgsubGMG->setQuantizationLevels(nlevels)

# Returns the kernel radius used for morphological operations
getSmoothingRadius() = @cxx bkgsubGMG->getSmoothingRadius()

# Sets the kernel radius used for morphological operations
setSmoothingRadius(radius::Int) = @cxx bkgsubGMG->setSmoothingRadius(radius)

# Returns the status of background model update
getUpdateBackgroundModel() = @cxx bkgsubGMG->getUpdateBackgroundModel()

# Sets the status of background model update
setUpdateBackgroundModel(update::Bool) = @cxx bkgsubGMG->setUpdateBackgroundModel(update)

# Returns the minimum value taken on by pixels in image sequence. Usually 0
getMinVal() = @cxx bkgsubGMG->getMinVal()

# Sets the minimum value taken on by pixels in image sequence
setMinVal(val::Float64) = @cxx bkgsubGMG->setMinVal(val)

# Returns the maximum value taken on by pixels in image sequence. e.g. 1.0 or 255
getMaxVal() = @cxx bkgsubGMG->getMaxVal()

# Sets the maximum value taken on by pixels in image sequence
setMaxVal(val::Float64) = @cxx bkgsubGMG->setMaxVal(val)

# Returns the prior probability that each individual pixel is a background pixel
getBackgroundPrior() = @cxx bkgsubGMG->getBackgroundPrior()

# Sets the prior probability that each individual pixel is a background pixel
setBackgroundPrior(bgprior::Float64) = @cxx bkgsubGMG->setBackgroundPrior(bgprior)

#-------------------------------------------------------------------------------------------------------------------#
# createOptFlow_DualTVL1
# “Dual TV L1” optical flow algorithm described in [Zach2007] and [Javier2012]
#  See http://stackoverflow.com/questions/19309567/speeding-up-optical-flow-createoptflow-dualtvl1
# Ptr<DenseOpticalFlow> createOptFlow_DualTVL1()
createOptFlow_DualTVL1() = @cxxnew cv::createOptFlow_DualTVL1()

# DenseOpticalFlow::calc
# void DenseOpticalFlow::calc(InputArray I0, InputArray I1, InputOutputArray flow)
# Parameters:
# I0     – first 8-bit single-channel input image
# I1     – second input image of the same size and the same type as I0
# flow   – computed flow image that has the same size as I0 and type CV_32FC2

calcDenseOpticalFLow(tvl1, I0, I1, flow) = @cxx tvl1->calc(I0, I1, flow)

# DenseOpticalFlow::collectGarbage
# void DenseOpticalFlow::collectGarbage()
collectGarbage(tvl1) = tvl1->collectGarbage()

# how to use:
# tvl1 = createOptFlow_DualTVL1()
# tvl1->calc(Previous_Gray_Frame, Current_Gray_Frame, Optical_Flow)

# createOptFlow_DualTVL1 parameters (class members)
# double tau
    # Time step of the numerical scheme
    # tvl1->getDouble("tau");
    # tvl1->set("tau",0.125);

# double lambda
    # Weight parameter for the data term, attachment parameter.
    # This is the most relevant parameter, which determines the smoothness of the output.
    # The smaller this parameter is, the smoother the solutions we obtain.
    # It depends on the range of motions of the images, so its value should be adapted to each image sequence.
    # tvl1->getDouble("lambda");
    # tvl1->set("lambda",0.125);

# double theta
    # Weight parameter for (u - v)^2, tightness parameter.
    # It serves as a link between the attachment and the regularization terms.
    # In theory, it should have a small value in order to maintain both parts in correspondence.
    # The method is stable for a large range of values of this parameter.
    # tvl1->getDouble("theta");
    # tvl1->set("theta",0.125);

# int nscales
    # Number of scales used to create the pyramid of images.
    # tvl1->getInt("nscales");
    # tvl1->set("nscales",2);

# int warps
    # Number of warpings per scale. Represents the number of times that I1(x+u0) and grad( I1(x+u0) ) are computed per scale.
    # This is a parameter that assures the stability of the method.
    # It also affects the running time, so it is a compromise between speed and accuracy.
    # tvl1->getInt("warps");
    # tvl1->set("warps",2);

# double epsilon
    # Stopping criterion threshold used in the numerical scheme, which is a trade-off between precision and running time.
    # A small value will yield more accurate solutions at the expense of a slower convergence.
    # tvl1->getDouble("epsilon");
    # tvl1->set("epsilon",0.5);

# int iterations
    # Stopping criterion iterations number used in the numerical scheme.
    # tvl1->getInt("iterations");
    # tvl1->set("iterations",20);


text:
#Abstract supertypes
"BootMethod <- Abstract supertype for all dependent bootstrap methods"
abstract type BootMethod end
"BlockLengthMethod <- Abstract supertype for all block length selection methods"
abstract type BlockLengthMethod end
"BandwidthMethod <- Abstract supertype for all bandwidth selection methods"
abstract type BandwidthMethod end
"KernelFunctionMethod <- Abstract supertype for all kernel function methods"
abstract type KernelFunctionMethod end

#Kernel function types
"KernelDummy <- Dummy type for kernel functions"
struct KernelDummy <: KernelFunctionMethod ; end
struct KernelTrap <: KernelFunctionMethod ; end
struct KernelSmooth <: KernelFunctionMethod ; end
Base.show(io::IO, x::KernelDummy) = print(io, "Dummy kernel function method")
Base.show(io::IO, x::KernelTrap) = print(io, "Trapezoidal kernel function method from Paparoditis, Politis (2002)")
Base.show(io::IO, x::KernelSmooth) = print(io, "Smooth kernel function method from Paparoditis, Politis (2002)")
kernel_function_dict_input = Pair{Symbol,KernelFunctionMethod}[
	:dummy => KernelDummy(),
	:trap => KernelTrap(),
	:smooth => KernelSmooth()
]::Vector{Pair{Symbol,KernelFunctionMethod}}
"KERNEL_FUNCTION_DICT <- Dictionary for converting string or symbol inputs into kernel function types"
const KERNEL_FUNCTION_DICT = Dict{Union{Symbol,String},KernelFunctionMethod}()::Dict{Union{Symbol,String},KernelFunctionMethod}
sizehint!(KERNEL_FUNCTION_DICT, 2*length(kernel_function_dict_input))
for kf in kernel_function_dict_input ; KERNEL_FUNCTION_DICT[kf[1]] = kf[2] ; end
for kf in kernel_function_dict_input ; KERNEL_FUNCTION_DICT[string(kf[1])] = kf[2] ; end

#Bootstrap method types
"BootDummy <- Dummy type used within the module. Should never be seen by the end user"
struct BootDummy <: BootMethod ; end
"BootIID <- Type for using multiple dispatch to get the IID boostrap"
struct BootIID <: BootMethod ; end
"BootStationary <- Type for using multiple dispatch to get the stationary bootstrap"
struct BootStationary <: BootMethod ; end
"BootMoving <- Type for using multiple dispatch to get the moving blocks bootstrap"
struct BootMoving <: BootMethod ; end
"BootNoOverlap <- Type for using multiple dispatch to get the non-overlapping blocks bootstrap"
struct BootNoOverlap <: BootMethod ; end
"BootCircular <- Type for using multiple dispatch to get the circular blocks bootstrap"
struct BootCircular <: BootMethod ; end
# "BootTapered <- Type for using multiple dispatch to get the tapered block bootstrap"
# struct BootTapered{Tkf<:KernelFunctionMethod} <: BootMethod
# 	kernelfunction::Tkf
# 	function BootTapered(kf::Tkf) where {Tkf}
# 		new{Tkf}(kf)
# 	end
# end
# BootTapered()::BootTapered = BootTapered(KernelTrap())
Base.show(io::IO, bm::BootDummy) = print(io, "Dummy bootstrap method")
Base.show(io::IO, bm::BootIID) = print(io, "IID bootstrap")
Base.show(io::IO, bm::BootStationary) = print(io, "Stationary bootstrap")
Base.show(io::IO, bm::BootMoving) = print(io, "Moving block bootstrap")
Base.show(io::IO, bm::BootNoOverlap) = print(io, "Non-overlapping block bootstrap")
Base.show(io::IO, bm::BootCircular) = print(io, "Circular block bootstrap")
#Base.show(io::IO, bm::BootTapered) = print(io, "Tapered block bootstrap")
boot_method_dict_input = Pair{Symbol, BootMethod}[
	:iid => BootIID(),
	:efron => BootIID(),
	:stationary => BootStationary(),
	:movingblock => BootMoving(),
	:moving => BootMoving(),
	:nonoverlappingblock => BootNoOverlap(),
	:nooverlap => BootNoOverlap(),
	:circularblock => BootCircular(),
	:circular => BootCircular()]::Vector{Pair{Symbol,BootMethod}}
"BOOT_METHOD_DICT <- Dictionary for converting string or symbol inputs into bootstrap methods"
const BOOT_METHOD_DICT = Dict{Union{Symbol,String},BootMethod}()::Dict{Union{Symbol,String},BootMethod}
sizehint!(BOOT_METHOD_DICT, 2*length(boot_method_dict_input))
for bm in boot_method_dict_input ; BOOT_METHOD_DICT[bm[1]] = bm[2] ; end
for bm in boot_method_dict_input ; BOOT_METHOD_DICT[string(bm[1])] = bm[2] ; end

#bandwidth method types
"P2003 <- Type for using multiple dispatch to get the bandwidth selection procedure of Politis (2003)"
struct P2003 <: BandwidthMethod ; end
Base.show(io::IO, bw::P2003) = print(io, "Bandwidth selection of Politis (2003)")
bandwidth_method_dict_input = Pair{Symbol,BandwidthMethod}[
	:politis2003 => P2003(),
	:p2003 => P2003()
]::Vector{Pair{Symbol,BandwidthMethod}}
"BANDWIDTH_METHOD_DICT <- Dictionary for converting string or symbol inputs into bandwidth methods"
const BANDWIDTH_METHOD_DICT = Dict{Union{Symbol,String},BandwidthMethod}()::Dict{Union{Symbol,String},BandwidthMethod}
sizehint!(BANDWIDTH_METHOD_DICT, 2*length(bandwidth_method_dict_input))
for bw in bandwidth_method_dict_input ; BANDWIDTH_METHOD_DICT[bw[1]] = bw[2] ; end
for bw in bandwidth_method_dict_input ; BANDWIDTH_METHOD_DICT[string(bw[1])] = bw[2] ; end

#blocklength method types
"BLDummy <- Dummy type for block length method"
struct BLDummy <: BlockLengthMethod ; end
# "BLPP2002 <- Type for using multiple dispatch to get the block length selection procedure of Paparoditis and Politis (2002)"
# struct BLPP2002{Tbw<:BandwidthMethod,Tkf<:KernelFunctionMethod} <: BlockLengthMethod
# 	bandwidthmethod::Tbw
# 	kernelfunction::Tkf
# end
"BLPPW2009 <- Type for using multiple dispatch to get the block length selection procedure of Patton, Politis, and White (2009)"
struct BLPPW2009{Tbw<:BandwidthMethod} <: BlockLengthMethod
	bandwidthmethod::Tbw
end
BLPPW2009()::BLPPW2009{P2003} = BLPPW2009(P2003())
# Base.show(io::IO, bl::BLPP2002) = print(io, "Block length selection of Paparoditis and Politis (2002)")
Base.show(io::IO, bl::BLPPW2009) = print(io, "Block length selection of Patton, Politis, and White (2009)")
blocklength_method_dict_input = Pair{Symbol,BlockLengthMethod}[
	:ppw2009 => BLPPW2009(P2003())]::Vector{Pair{Symbol,BlockLengthMethod}}
"BLOCKLENGTH_METHOD_DICT <- Dictionary for converting string or symbol inputs into bandwidth methods"
const BLOCKLENGTH_METHOD_DICT = Dict{Union{Symbol,String},BlockLengthMethod}()::Dict{Union{Symbol,String},BlockLengthMethod}
sizehint!(BLOCKLENGTH_METHOD_DICT, 2*length(blocklength_method_dict_input))
for bl in blocklength_method_dict_input ; BLOCKLENGTH_METHOD_DICT[bl[1]] = bl[2] ; end
for bl in blocklength_method_dict_input ; BLOCKLENGTH_METHOD_DICT[string(bl[1])] = bl[2] ; end

"bootmethod_to_blocklengthmethod <- Convert a bootstrap method to the most appropriate block length method"
(bootmethod_to_blocklengthmethod(bm::T)::BLPPW2009{P2003}) where {T<:Union{BootStationary,BootMoving,BootCircular,BootNoOverlap,BootIID}} = BLPPW2009(P2003())
#(bootmethod_to_blocklengthmethod(bm::BootTapered)::BLPP2002{P2003,KernelDummy}) = BLPP2002(P2003(), KernelDummy())

"""
	BootInput

Core type that defines all parameters needed to perform a bootstrap procedure. The
vast majority of users should use the keyword argument constructor that has the
method signature:

    BootInput(data ; kwargs...)

where data is the dataset to be bootstrapped, and kwargs denotes a set of keyword arguments
(defined below) that are used for every exported function in the DependentBootstrap package.
The following keyword arguments and default values follow: \n
    - blocklength         <- Block length for bootstrapping procedure. Default value is 0. Set to <= 0 to auto-estimate the optimal block length from the dataset. Float64 inputs allowed.
	- numresample         <- Number of times to resample the input dataset. Default value is the module constant NUM_RESAMPLE, currently set to 1000.
    - bootmethod          <- Bootstrapping methodology to use. Default value is the Symbol :stationary (for the stationary bootstrap).
	- blocklengthmethod   <- Block length selection procedure to use if user wishes to auto-estimate the block length. Default value is the Symbol :ppw2009 (use the method described in Patton, Politis, and White (2009)).
	- flevel1             <- A function that converts the input dataset to the estimator that the user wishes to bootstrap. Default value is the sample mean.
	- flevel2             <- A function that converts a vector of estimators constructed by flevel1 into a distributional parameter. Default value is sample variance.
	- numobsperresample   <- Number of observations to be drawn (with replacement) per resample. The default value is the number of observations in the dataset (the vast majority of users will want this default value).
	- fblocklengthcombine <- A function for converting a Vector{Float64} of estimated blocklengths to a single Float64 blocklength estimate. Default value is median.

The constructor will attempt to convert all provided keyword arguments to appropriate types,
and will notify the user via an error if a supplied keyword argument is not valid.

Note that the bootmethod and blocklengthmethod keyword arguments will accept both
Symbol and String inputs, and will convert them to BootMethod and BlockLengthMethod
types internally. To see a list of acceptable Symbol or String values for the bootmethod and
blocklengthmethod keyword arguments, use: \n
	- collect(keys(DependentBootstrap.BOOT_METHOD_DICT))
	- collect(keys(DependentBootstrap.BLOCKLENGTH_METHOD_DICT))

respectively. A small proportion of users may need the fine-grained control that comes
from constructing BootMethod and BlockLengthMethod types explicitly and then providing
them to the keyword constructor. These users should use ?BootMethod and ?BlockLengthMethod
at the REPL for more info.

BootInput is not mutable, but the type is near instantaneous to construct, so if a user wishes
to amend a BootInput it is recommended to just construct another one. A special constructor
is provided to facilitate this process that has the method definition: \n
	- BootInput(data, bootinput::BootInput ; kwargs...)

where the new BootInput draws its fields from the keyword arguments that are provided, and then
the input BootInput for any keyword arguments that are not provided.

Note that all exported functions in the DependentBootstrap package exhibit the
method signature: \n
    - exported_func(data ; kwargs...)

which in practice just wraps the keyword argument constructor for a BootInput, and
then calls the method signature: \n
     - exported_func(data, bootinput::BootInput)
"""
struct BootInput{Tbm<:BootMethod,Tbl<:BlockLengthMethod,Tf1<:Function,Tf2<:Function,Tfc<:Function}
	numobs::Int
    blocklength::Float64
    numresample::Int
	bootmethod::Tbm
	blocklengthmethod::Tbl
	flevel1::Tf1
	flevel2::Tf2
    numobsperresample::Int
	fblocklengthcombine::Tfc
	function BootInput(numobs::Int, blocklength::Float64, numresample::Int, bootmethod::Tbm, blmethod::Tbl,
                       flevel1::Tf1, flevel2::Tf2, numobsperresample::Int, fblocklengthcombine::Tfc) where {Tbm<:BootMethod,Tbl<:BlockLengthMethod,Tf1<:Function,Tf2<:Function,Tfc<:Function}
		typeof(bootmethod) <: BootDummy && error("bootmethod is set to BootDummy. You should not have been able to accidentally reach this point. Please file an issue.")
		numobs < 2 && error("Number of observations input to BootInput must be 2 or greater: $(numobs)")
		(isnan(blocklength) || isinf(blocklength)) && error("Invalid blocklength: $(blocklength)")
		numresample < 1 && error("Number of resamples must be strictly positive: $(numresample)")
		numobsperresample < 1 && error("Number of observations per resample must be strictly positive: $(numobsperresample)")
		blocklength <= 0.0 && error("Block length must be strictly positive. BootInput outer constructors that include the dataset as the first argument will automatically estimate the block length if you do not specify it")
		Tbm == BootDummy && error("Do not use BootDummy as an input type. It is for internal module use only")
		Tbl == BLDummy && error("Do not use BLDummy as an input type. It is for internal module use only")
		new{Tbm,Tbl,Tf1,Tf2,Tfc}(numobs, blocklength, numresample, bootmethod, blmethod, flevel1, flevel2, numobsperresample, fblocklengthcombine)
	end
end
#BootInput empty constructor
BootInput() = BootInput(2, 1.0, 1, BootIID(), BLPPW2009(), mean, var, 1, identity)
#BootInput constructor that ensures blocklength is auto-detected if need be
function BootInput(data, numobs::Int, blocklength::Number, numresample::Int, bootmethod::Tbm, blocklengthmethod::Tbl, flevel1::Tf1,
                   flevel2::Tf2, numobsperresample::Int, fblocklengthcombine::Tfc) where {Tbm<:BootMethod,Tbl<:BlockLengthMethod,Tf1<:Function,Tf2<:Function,Tfc<:Function}
	blocklength = Float64(blocklength)
    Tbm <: BootIID && (blocklength = 1.0)
    blocklength <= 0.0 && (blocklength = optblocklength(data, blocklengthmethod, bootmethod, fblocklengthcombine))
	return BootInput(numobs, blocklength, numresample, bootmethod, blocklengthmethod, flevel1, flevel2, numobsperresample, fblocklengthcombine)
end
#BootInput constructors that use keyword arguments
function BootInput(data ; blocklength=0, numresample=NUM_RESAMPLE, bootmethod=:stationary, blocklengthmethod=:dummy,
				   flevel1=mean, flevel2=var, numobsperresample=num_obs(data), fblocklengthcombine=median)
	numobs = num_obs(data)
	blocklength = bootinput_get_blocklength(blocklength)
	numresample = bootinput_get_numresample(numresample)
	bootmethod = bootinput_get_bootmethod(bootmethod)
	blocklengthmethod = bootinput_get_blocklengthmethod(blocklengthmethod, bootmethod)
	flevel1 = bootinput_get_flevel1(flevel1)
	flevel2 = bootinput_get_flevel2(flevel2)
	numobsperresample = bootinput_get_numobsperresample(numobsperresample)
	fblocklengthcombine = bootinput_get_fblocklengthcombine(fblocklengthcombine)
	typeof(bootmethod) <: BootIID && (blocklength = 1.0)
	blocklength <= 0.0 && (blocklength = optblocklength(data, blocklengthmethod, bootmethod, fblocklengthcombine))
	return BootInput(numobs, blocklength, numresample, bootmethod, blocklengthmethod, flevel1, flevel2, numobsperresample, fblocklengthcombine)
end
#Constructor for building a new BootInput using keyword arguments where possible, and, failing that, the fields
#of an existing BootInput
_db_bi_dummy_f() = error("This function is designed to never be called. It is used as a default value for keyword arguments so as to check whether the user has specified them")
function BootInput(data, bi::BootInput ; blocklength=-9, numresample=-9, bootmethod=:z, blocklengthmethod=:z,
				   flevel1=_db_bi_dummy_f, flevel2=_db_bi_dummy_f, numobsperresample=-9,
				   fblocklengthcombine=_db_bi_dummy_f)
    blocklength == -9 && (blocklength = bi.blocklength)
	numresample == -9 && (numresample = bi.numresample)
	bootmethod == :z && (bootmethod = bi.bootmethod)
	blocklengthmethod == :z && (blocklengthmethod = bi.blocklengthmethod)
	flevel1 == _db_bi_dummy_f && (flevel1 = bi.flevel1)
	flevel2 == _db_bi_dummy_f && (flevel2 = bi.flevel2)
	numobsperresample == -9 && (numobsperresample = bi.numobsperresample)
	fblocklengthcombine == _db_bi_dummy_f && (fblocklengthcombine = bi.fblocklengthcombine)
	return BootInput(data, blocklength=blocklength, numresample=numresample, bootmethod=bootmethod, blocklengthmethod=blocklengthmethod,
					 flevel1=flevel1, flevel2=flevel2, numobsperresample=numobsperresample, fblocklengthcombine=fblocklengthcombine)
end
bootinput_get_blocklength(x::Number)::Float64 = Float64(x)
bootinput_get_blocklength(x) = error("Invalid type for blocklength input. Use a subtype of Number, eg Int or Float64.")
bootinput_get_numresample(x::Number)::Int = Int(x)
bootinput_get_numresample(x) = error("Invalid type for numresample input. Use a subtype of Number, eg Int or Float64.")
function bootinput_get_bootmethod(bootmethod::T) where {T<:Union{Symbol,String}}
	bm = get(BOOT_METHOD_DICT, bootmethod, BootDummy())
	typeof(bm) <: BootDummy && error("No matching entry found in dictionary for bootmethod input $(bootmethod). Please use collect(keys(DependentBootstrap.BOOT_METHOD_DICT)) at the REPL to see a list of valid keyword arguments.")
	return bm
end
bootinput_get_bootmethod(x) = typeof(x) <: BootMethod ? x : error("Invalid type for bootmethod input. Use Symbol, String, or a subtype of BootMethod")
function bootinput_get_blocklengthmethod(blocklengthmethod::T, bootmethod::Tbm) where {T<:Union{Symbol,String},Tbm<:BootMethod}
	string(blocklengthmethod) == "dummy" && return bootmethod_to_blocklengthmethod(bootmethod)
	blm = get(BLOCKLENGTH_METHOD_DICT, blocklengthmethod, BLDummy())
	typeof(blm) <: BLDummy && error("No matching entry found in dictionary for blocklengthmethod input $(blocklengthmethod). Please use collect(keys(DependentBootstrap.BLOCKLENGTH_METHOD_DICT)) at the REPL to see a list of valid keyword arguments.")
	return blm
end
bootinput_get_blocklengthmethod(x, bootmethod) = typeof(x) <: BlockLengthMethod ? x : error("Invalid type for bootmethod input. Use Symbol, String, or a subtype of BootMethod")
(bootinput_get_flevel1(f::Tf)::Tf) where {Tf<:Function} = f
bootinput_get_flevel1(f) = error("Invalid type for flevel1 input. Use a subtype of Function")
(bootinput_get_flevel2(f::Tf)::Tf) where {Tf<:Function} = f
bootinput_get_flevel2(f) = error("Invalid type for flevel2 input. Use a subtype of Function")
bootinput_get_numobsperresample(x::Number)::Int = Int(x)
bootinput_get_numobsperresample(x) = error("Invalid type for numobsperresample input. Use a subtype of Number, eg Int or Float64.")
(bootinput_get_fblocklengthcombine(f::Tf)::Tf) where {Tf<:Function} = f
bootinput_get_fblocklengthcombine(f) = error("Invalid type for fblocklengthcombine input. Use a subtype of Function")
function Base.show(io::IO, x::BootInput)
	println(io, "Dependent bootstrap input:")
	println(io, "    Number of observations in dataset = $(x.numobs)")
    println(io, "    Current block length = $(x.blocklength)")
	println(io, "    Number of resamples = $(x.numresample)")
	println(io, "    Bootstrap method = $(x.bootmethod)")
	println(io, "    Block length method = $(x.blocklengthmethod)")
	println(io, "    Level 1 statistic = $(x.flevel1)")
	println(io, "    Level 2 statistic = $(x.flevel2)")
	println(io, "    Number of observations per resample = $(x.numobsperresample)")
	println(io, "    Block length combine function = $(x.fblocklengthcombine)")
end

ext: jl


text:
"""
wage Income

any fiIncome provided by
from a(ny) Busniess to Personal account
for the
(can  includeable to  cuts or bonuses)
"""

mutable struct wage :> Income

function  wage(amount)

end

"""
cut can durablee over a few periods of tim

"""
function wagecut(amount)

end

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# The policy type
"""
    ValueIterationPolicy <: Policy

The policy type. Contains the Q-Matrix, the Utility function and an array of indices corresponding to optimal actions.
There are three ways to initialize the policy type:

    `policy = ValueIterationPolicy(mdp)`
    `policy = ValueIterationPolicy(mdp, utility_array)`
    `policy = ValueIterationPolicy(mdp, qmatrix)`

The Q-matrix is nxm, where n is the number of states and m is the number of actions.

# Fields
- `qmat`  Q matrix storing Q(s,a) values
- `util` The value function V(s)
- `policy` Policy array, maps state index to action index
- `action_map` Maps the action index to the concrete action type
- `include_Q` Flag for including the Q-matrix
- `mdp`  uses the model for indexing in the action function
"""
struct ValueIterationPolicy{Q<:AbstractMatrix, U<:AbstractVector, P<:AbstractVector, A, M<:MDP} <: Policy
    qmat::Q
    util::U
    policy::P
    action_map::Vector{A}
    include_Q::Bool
    mdp::M
end

# constructor with an optinal initial value function argument
function ValueIterationPolicy(mdp::Union{MDP,POMDP};
                              utility::AbstractVector{Float64}=zeros(length(states(mdp))),
                              policy::AbstractVector{Int64}=zeros(Int64, length(states(mdp))),
                              include_Q::Bool=true)
    ns = length(states(mdp))
    na = length(actions(mdp))
    @assert length(utility) == ns "Input utility dimension mismatch"
    @assert length(policy) == ns "Input policy dimension mismatch"
    action_map = ordered_actions(mdp)
    include_Q ? qmat = zeros(ns,na) : qmat = zeros(0,0)
    return ValueIterationPolicy(qmat, utility, policy, action_map, include_Q, mdp)
end

# constructor for solved q, util and policy
function ValueIterationPolicy(mdp::Union{MDP,POMDP}, q::AbstractMatrix{F}, util::AbstractVector{F}, policy::Vector{Int64}) where {F}
    action_map = ordered_actions(mdp)
    include_Q = true
    return ValueIterationPolicy(q, util, policy, action_map, include_Q, mdp)
end

# constructor for default Q-matrix
function ValueIterationPolicy(mdp::Union{MDP,POMDP}, q::AbstractMatrix{F}) where {F}
    (ns, na) = size(q)
    p = zeros(Int64, ns)
    u = zeros(ns)
    for i = 1:ns
        p[i] = argmax(q[i,:])
        u[i] = maximum(q[i,:])
    end
    action_map = ordered_actions(mdp)
    include_Q = true
    return ValueIterationPolicy(q, u, p, action_map, include_Q, mdp)
end

# returns the fields of the policy type
function locals(p::ValueIterationPolicy)
    return (p.qmat,p.util,p.policy,p.action_map)
end

function action(policy::ValueIterationPolicy, s::S) where S
    sidx = stateindex(policy.mdp, s)
    aidx = policy.policy[sidx]
    return policy.action_map[aidx]
end

function value(policy::ValueIterationPolicy, s::S) where S
    sidx = stateindex(policy.mdp, s)
    policy.util[sidx]
end

value(policy::ValueIterationPolicy, s, a) = actionvalues(policy, s)[actionindex(policy.mdp, a)]

function POMDPPolicies.actionvalues(policy::ValueIterationPolicy, s::S) where S
    if !policy.include_Q
        error("ValueIterationPolicy does not contain the Q matrix. Use the include_Q=true keyword argument in the solver.")
    else
        sidx = stateindex(policy.mdp, s)
        return policy.qmat[sidx,:]
    end
end

function Base.show(io::IO, mime::MIME"text/plain", p::ValueIterationPolicy)
    println(io, "ValueIterationPolicy:")
    ds = get(io, :displaysize, displaysize(io))
    ioc = IOContext(io, :displaysize=>(first(ds)-1, last(ds)))
    showpolicy(ioc, mime, p.mdp, p)
end


text:
@testset "get_variables_str" begin
    client = Client("url","ws"; introspect=false)

    # Simple test of field
    client.query_to_args_map["Query"] = Dict("fieldName" => "Int")
    args = Dict("fieldName" => 1)
    arg_names = Dict("fieldName" => "fieldName1")
    @test GraphQLClient.get_variables_str(client, args, arg_names, "Query") == "\$fieldName1: Int,"

    # Simple test of input object
    client.query_to_args_map["Query"] = Dict("fieldName" => "InputObject")
    client.input_object_fields_to_type_map["InputObject"] = Dict("subfieldName" => "Int")
    args = Dict("fieldName" => Dict("subfieldName" => 1))
    arg_names = Dict("fieldName" => Dict("subfieldName" => "subfieldName1"))
    @test GraphQLClient.get_variables_str(client, args, arg_names, "Query") == "\$subfieldName1: Int,"

    # Simple test of non null input object
    client.query_to_args_map["Query"] = Dict("fieldName" => "InputObject!")
    client.input_object_fields_to_type_map["InputObject"] = Dict("subfieldName" => "Int")
    args = Dict("fieldName" => Dict("subfieldName" => 1))
    arg_names = Dict("fieldName" => Dict("subfieldName" => "subfieldName1"))
    @test GraphQLClient.get_variables_str(client, args, arg_names, "Query") == "\$subfieldName1: Int,"

    # Simple test of vector of input objects
    client.query_to_args_map["Query"] = Dict("fieldName" => "[InputObject]")
    client.input_object_fields_to_type_map["InputObject"] = Dict("subfieldName" => "Int")
    args = Dict("fieldName" => [Dict("subfieldName" => 1), Dict("subfieldName" => 1)])
    arg_names = Dict("fieldName" => [
        Dict("subfieldName" => "subfieldName1"),
        Dict("subfieldName" => "subfieldName2")
    ])
    @test GraphQLClient.get_variables_str(client, args, arg_names, "Query") == "\$subfieldName1: Int,\$subfieldName2: Int,"
end

@testset "get_query_variables" begin
    # Simple test of field
    args = Dict("fieldName" => 1)
    arg_names = Dict("fieldName" => "fieldName1")
    variables = GraphQLClient.get_query_variables(args, arg_names)
    @test haskey(variables, "fieldName1")
    @test variables["fieldName1"] == 1

    # Simple test of input object
    args = Dict("fieldName" => Dict("subfieldName" => true))
    arg_names = Dict("fieldName" => Dict("subfieldName" => "subfieldName1"))
    variables = GraphQLClient.get_query_variables(args, arg_names)
    @test haskey(variables, "subfieldName1")
    @test variables["subfieldName1"] == true

    # Simple test of vector of input objects
    args = Dict("fieldName" => [Dict("subfieldName" => 1), Dict("subfieldName" => 2)])
    arg_names = Dict("fieldName" => [
        Dict("subfieldName" => "subfieldName1"),
        Dict("subfieldName" => "subfieldName2")
    ])
    variables = GraphQLClient.get_query_variables(args, arg_names)
    @test haskey(variables, "subfieldName1")
    @test haskey(variables, "subfieldName2")
    @test variables["subfieldName1"] == 1
    @test variables["subfieldName2"] == 2

    # Test lots of repeated fields
    args = Dict(
        "arg1" => 1,
        "arg2" => Dict(
            "arg1" => 2,
            "arg2" => [
                Dict(
                    "arg1" => 6.0,
                    "arg2" => 5.0,
                ),
                Dict(
                    "arg1" => 7.0,
                    "arg2" => 2.0,
                ),
            ]
        )
    )
    _, arg_names = GraphQLClient.get_query_args_str(args)
    variables = GraphQLClient.get_query_variables(args, arg_names)
    @test variables[arg_names["arg1"]] == 1
    @test variables[arg_names["arg2"]["arg1"]] == 2
    @test variables[arg_names["arg2"]["arg2"][1]["arg1"]] == 6.0
    @test variables[arg_names["arg2"]["arg2"][1]["arg2"]] == 5.0
    @test variables[arg_names["arg2"]["arg2"][2]["arg1"]] == 7.0
    @test variables[arg_names["arg2"]["arg2"][2]["arg2"]] == 2.0
end


text:
using Cairo
using LSystems
using DataStructures

c = CairoRGBSurface(1024, 1024);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 1024.0, 1024.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 260., 760., 0, 2.0
set_line_width(cr, 0.6);
move_to(cr, x, y);

koch_start = "F+F+F+F"
koch_trans = Dict([('F', "F+F-F-FF+F+F-F")])

function determine_new_position(x, y, a, d)
    x -= sind(a) * d
    y -= cosd(a) * d
    return x, y
end

s = Stack(Vector{Float64}) # [x, y, a]

level = 4
drawing_instructions = @task lindenmayer(koch_start, koch_trans, level)

for v in drawing_instructions
    if v == 'F'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= 90
    elseif v == '-'
        a += 90
    elseif v == '['
        push!(s, [x, y, a])
    elseif v == ']'
        x, y, a = pop!(s)
        move_to(cr, x, y)
    end
end

stroke(cr);
write_to_png(c, "koch_$level.png");


text:
@doc raw"""
    AbstractMetric

Abstract type for the pseudo-Riemannian metric tensor ``g``, a family of smoothly
varying inner products on the tangent space. See [`inner`](@ref).

# Functor

    (metric::Metric)(M::Manifold)

Generate the `MetricManifold` that wraps the manifold `M` with given `metric`.
This works for both a variable containing the metric as well as a subtype `T<:AbstractMetric`,
where a zero parameter constructor `T()` is availabe.
"""
abstract type AbstractMetric end

# piping syntax for decoration
(metric::AbstractMetric)(M::AbstractManifold) = MetricManifold(M, metric)
(::Type{T})(M::AbstractManifold) where {T<:AbstractMetric} = MetricManifold(M, T())

"""
    MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: AbstractDecoratorManifold{𝔽}

Equip a [`AbstractManifold`](@ref) explicitly with a [`AbstractMetric`](@ref) `G`.

For a Metric AbstractManifold, by default, assumes, that you implement the linear form
from [`local_metric`](@ref) in order to evaluate the exponential map.

If the corresponding [`AbstractMetric`](@ref) `G` yields closed form formulae for e.g.
the exponential map and this is implemented directly (without solving the ode),
you can of course still implement that directly.

# Constructor

    MetricManifold(M, G)

Generate the [`AbstractManifold`](@ref) `M` as a manifold with the [`AbstractMetric`](@ref) `G`.
"""
struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <:
       AbstractDecoratorManifold{𝔽}
    manifold::M
    metric::G
end

@doc raw"""
    RiemannianMetric <: AbstractMetric

Abstract type for Riemannian metrics, a family of positive definite inner
products. The positive definite property means that for ``X  ∈ T_p \mathcal M``, the
inner product ``g(X, X) > 0`` whenever ``X`` is not the zero vector.
"""
abstract type RiemannianMetric <: AbstractMetric end

@doc raw"""
    christoffel_symbols_first(
        M::MetricManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Compute the Christoffel symbols of the first kind in local coordinates of basis `B`.
The Christoffel symbols are (in Einstein summation convention)

```math
Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr],
```

where ``g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}`` is the coordinate
derivative of the local representation of the metric tensor. The dimensions of
the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
christoffel_symbols_first(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_first(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    ∂g = local_metric_jacobian(M, p, B; backend=backend)
    n = size(∂g, 1)
    Γ = allocate(∂g, Size(n, n, n))
    @einsum Γ[i, j, k] = 1 / 2 * (∂g[k, j, i] + ∂g[i, k, j] - ∂g[i, j, k])
    return Γ
end
@decorator_transparent_signature christoffel_symbols_first(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(christoffel_symbols_first),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

@doc raw"""
    christoffel_symbols_second(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Compute the Christoffel symbols of the second kind in local coordinates of basis `B`.
The Christoffel symbols are (in Einstein summation convention)

````math
Γ^{l}_{ij} = g^{kl} Γ_{ijk},
````

where ``Γ_{ijk}`` are the Christoffel symbols of the first kind, and
``g^{kl}`` is the inverse of the local representation of the metric tensor.
The dimensions of the resulting multi-dimensional array are ordered ``(l,i,j)``.
"""
christoffel_symbols_second(::AbstractManifold, ::AbstractBasis, ::Any)
function christoffel_symbols_second(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Γ₁ = christoffel_symbols_first(M, p, B; backend=backend)
    Γ₂ = allocate(Γ₁)
    @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k]
    return Γ₂
end
@decorator_transparent_signature christoffel_symbols_second(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(christoffel_symbols_second),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

@doc raw"""
    christoffel_symbols_second_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Get partial derivatives of the Christoffel symbols of the second kind
for manifold `M` at `p` with respect to the coordinates of `B`, i.e.

```math
\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}.
```

The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k,l)``.
"""
christoffel_symbols_second_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_second_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    ∂Γ = reshape(
        _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend),
        n,
        n,
        n,
        n,
    )
    return ∂Γ
end
@decorator_transparent_signature christoffel_symbols_second_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(christoffel_symbols_second_jacobian),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

Base.copyto!(M::MetricManifold, q, p) = copyto!(M.manifold, q, p)
Base.copyto!(M::MetricManifold, Y, p, X) = copyto!(M.manifold, Y, p, X)

decorator_transparent_dispatch(::typeof(exp), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(exp!), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(exp!), M::MetricManifold, q, p, X, t) = Val(:parent)
decorator_transparent_dispatch(::typeof(flat), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(flat!), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(
    ::typeof(get_coordinates),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(get_coordinates!),
    M::MetricManifold,
    args...,
)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(::typeof(get_vector), M::MetricManifold, args...)
    return Val(:parent)
end
function decorator_transparent_dispatch(::typeof(get_vector!), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(::typeof(get_basis), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(::typeof(inner), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(
    ::typeof(inverse_retract),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(inverse_retract!),
    M::MetricManifold,
    args...,
)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(
    ::typeof(inverse_retract!),
    M::MetricManifold,
    X,
    p,
    q,
    m::LogarithmicInverseRetraction,
)
    return Val(:parent)
end
decorator_transparent_dispatch(::typeof(log), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(log!), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(mean), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(mean!), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(median), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(median!), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(::typeof(norm), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(project), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(project!), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(sharp), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(sharp!), M::MetricManifold, args...)
    return Val(:intransparent)
end
decorator_transparent_dispatch(::typeof(retract), M::MetricManifold, args...) = Val(:parent)
function decorator_transparent_dispatch(::typeof(retract!), M::MetricManifold, args...)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(
    ::typeof(retract!),
    M::MetricManifold,
    q,
    p,
    X,
    m::ExponentialRetraction,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_along),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_along!),
    M::MetricManifold,
    args...,
)
    return Val(:intransparent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_direction),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_direction!),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_to),
    M::MetricManifold,
    args...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(vector_transport_to!),
    M::MetricManifold,
    args...,
)
    return Val(:intransparent)
end

@doc raw"""
    det_local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the determinant of local matrix representation of the metric tensor ``g``, i.e. of the
matrix ``G(p)`` representing the metric in the tangent space at ``p`` with as a matrix.

See also [`local_metric`](@ref)
"""
det_local_metric(::AbstractManifold, p, ::AbstractBasis)
function det_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return det(local_metric(M, p, B))
end
@decorator_transparent_signature det_local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)
function decorator_transparent_dispatch(
    ::typeof(det_local_metric),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end
"""
    einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Einstein tensor of the manifold `M` at the point `p`, see [https://en.wikipedia.org/wiki/Einstein_tensor](https://en.wikipedia.org/wiki/Einstein_tensor)
"""
einstein_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function einstein_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ric = ricci_tensor(M, p, B; backend=backend)
    g = local_metric(M, p, B)
    Ginv = inverse_local_metric(M, p, B)
    S = sum(Ginv .* Ric)
    G = Ric - g .* S / 2
    return G
end
@decorator_transparent_signature einstein_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(einstein_tensor),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

@doc raw"""
    exp(N::MetricManifold{M,G}, p, X)

Copute the exponential map on the [`AbstractManifold`](@ref) `M` equipped with the [`AbstractMetric`](@ref) `G`.

If the metric was declared the default metric using [`is_default_metric`](@ref), this method
falls back to `exp(M, p, X)`.

Otherwise it numerically integrates the underlying ODE, see [`solve_exp_ode`](@ref).
Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.
"""
exp(::MetricManifold, ::Any...)

@decorator_transparent_fallback function exp!(M::MetricManifold, q, p, X)
    tspan = (0.0, 1.0)
    A = get_default_atlas(M)
    i = get_chart_index(M, A, p)
    B = induced_basis(M, A, i, TangentSpace)
    sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0])
    n = length(p)
    return copyto!(q, sol.u[1][(n + 1):end])
end

@doc raw"""
    flat(N::MetricManifold{M,G}, p, X::FVector{TangentSpaceType})

Compute the musical isomorphism to transform the tangent vector `X` from the
[`AbstractManifold`](@ref) `M` equipped with [`AbstractMetric`](@ref) `G` to a cotangent by
computing

````math
X^♭= G_p X,
````
where ``G_p`` is the local matrix representation of `G`, see [`local_metric`](@ref)
"""
flat(::MetricManifold, ::Any...)

@decorator_transparent_fallback function flat!(
    M::MetricManifold,
    ξ::CoTFVector,
    p,
    X::TFVector,
)
    g = local_metric(M, p, ξ.basis)
    copyto!(ξ.data, g * X.data)
    return ξ
end

"""
    gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`.
"""
gaussian_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    return ricci_curvature(M, p, B; kwargs...) / 2
end
@decorator_transparent_signature gaussian_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(gaussian_curvature),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

function injectivity_radius(M::MetricManifold, p)
    return injectivity_radius(base_manifold(M), p)
end
function injectivity_radius(M::MetricManifold, m::AbstractRetractionMethod)
    return injectivity_radius(base_manifold(M), m)
end
function injectivity_radius(M::MetricManifold, m::ExponentialRetraction)
    return injectivity_radius(base_manifold(M), m)
end
function injectivity_radius(M::MetricManifold, p, m::AbstractRetractionMethod)
    return injectivity_radius(base_manifold(M), p, m)
end
function injectivity_radius(M::MetricManifold, p, m::ExponentialRetraction)
    return injectivity_radius(base_manifold(M), p, m)
end

@doc raw"""
    inverse_local_metric(M::AbstractcManifold, p, B::AbstractBasis)

Return the local matrix representation of the inverse metric (cometric) tensor, usually
written ``g^{ij}``.

See also [`local_metric`](@ref)
"""
inverse_local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
function inverse_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return inv(local_metric(M, p, B))
end
@decorator_transparent_signature inverse_local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)
function decorator_transparent_dispatch(
    ::typeof(inverse_local_metric),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

default_decorator_dispatch(M::MetricManifold) = default_metric_dispatch(M)

"""
    is_default_metric(M, G)

Indicate whether the [`AbstractMetric`](@ref) `G` is the default metric for
the [`AbstractManifold`](@ref) `M`. This means that any occurence of
[`MetricManifold`](@ref)(M,G) where `typeof(is_default_metric(M,G)) = true`
falls back to just be called with `M` such that the [`AbstractManifold`](@ref) `M`
implicitly has this metric, for example if this was the first one implemented
or is the one most commonly assumed to be used.
"""
function is_default_metric(M::AbstractManifold, G::AbstractMetric)
    return _extract_val(default_metric_dispatch(M, G))
end

default_metric_dispatch(::AbstractManifold, ::AbstractMetric) = Val(false)
function default_metric_dispatch(M::MetricManifold)
    return default_metric_dispatch(base_manifold(M), metric(M))
end

"""
    is_default_metric(MM::MetricManifold)

Indicate whether the [`AbstractMetric`](@ref) `MM.G` is the default metric for
the [`AbstractManifold`](@ref) `MM.manifold,` within the [`MetricManifold`](@ref) `MM`.
This means that any occurence of
[`MetricManifold`](@ref)`(MM.manifold, MM.G)` where `is_default_metric(MM.manifold, MM.G)) = true`
falls back to just be called with `MM.manifold,` such that the [`AbstractManifold`](@ref) `MM.manifold`
implicitly has the metric `MM.G`, for example if this was the first one
implemented or is the one most commonly assumed to be used.
"""
function is_default_metric(M::MetricManifold)
    return _extract_val(default_metric_dispatch(M))
end

function Base.convert(::Type{MetricManifold{𝔽,MT,GT}}, M::MT) where {𝔽,MT,GT}
    return _convert_with_default(M, GT, default_metric_dispatch(M, GT()))
end

function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{true},
) where {MT<:AbstractManifold}
    return MetricManifold(M, T())
end
function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{false},
) where {MT<:AbstractManifold}
    return error(
        "Can not convert $(M) to a MetricManifold{$(MT),$(T)}, since $(T) is not the default metric.",
    )
end

@doc raw"""
    inner(N::MetricManifold{M,G}, p, X, Y)

Compute the inner product of `X` and `Y` from the tangent space at `p` on the
[`AbstractManifold`](@ref) `M` using the [`AbstractMetric`](@ref) `G`. If `G` is the default
metric (see [`is_default_metric`](@ref)) this is done using `inner(M, p, X, Y)`,
otherwise the [`local_metric`](@ref)`(M, p)` is employed as

````math
g_p(X, Y) = ⟨X, G_p Y⟩,
````
where ``G_p`` is the loal matrix representation of the [`AbstractMetric`](@ref) `G`.
"""
inner(::MetricManifold, ::Any, ::Any, ::Any)

@decorator_transparent_fallback :intransparent function inner(
    M::MetricManifold,
    p,
    X::TFVector,
    Y::TFVector,
)
    X.basis === Y.basis ||
        error("calculating inner product of vectors from different bases is not supported")
    return dot(X.data, local_metric(M, p, X.basis) * Y.data)
end

@doc raw"""
    local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the local matrix representation at the point `p` of the metric tensor ``g`` with
respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`, usually written ``g_{ij}``.
The matrix has the property that ``g(X, Y)=X^\mathrm{T} [g_{ij}] Y = g_{ij} X^i Y^j``,
where the latter expression uses Einstein summation convention.
The metric tensor is such that the formula works for the given [`AbstractBasis`](@ref) `B`.
"""
local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
@decorator_transparent_signature local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(::typeof(local_metric), ::MetricManifold, args...)
    return Val(:parent)
end

@doc raw"""
    local_metric_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Get partial derivatives of the local metric of `M` at `p` in basis `B` with respect to the
coordinates of `p`, ``\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}``. The
dimensions of the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
function local_metric_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    ∂g = reshape(_jacobian(q -> local_metric(M, q, B), p, backend), n, n, n)
    return ∂g
end
@decorator_transparent_signature local_metric_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(local_metric_jacobian),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

@doc raw"""
    log(N::MetricManifold{M,G}, p, q)

Copute the logarithmic map on the [`AbstractManifold`](@ref) `M` equipped with the [`AbstractMetric`](@ref) `G`.

If the metric was declared the default metric using [`is_default_metric`](@ref), this method
falls back to `log(M,p,q)`. Otherwise, you have to provide an implementation for the non-default
[`AbstractMetric`](@ref) `G` metric within its [`MetricManifold`](@ref)`{M,G}`.
"""
log(::MetricManifold, ::Any...)

@doc raw"""
    log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)

Return the natural logarithm of the metric density ``ρ`` of `M` at `p`, which
is given by ``ρ = \log \sqrt{|\det [g_{ij}]|}`` for the metric tensor expressed in basis `B`.
"""
log_local_metric_density(::AbstractManifold, ::Any, ::AbstractBasis)
function log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)
    return log(abs(det_local_metric(M, p, B))) / 2
end
@decorator_transparent_signature log_local_metric_density(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)
function decorator_transparent_dispatch(
    ::typeof(log_local_metric_density),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end

@doc raw"""
    metric(M::MetricManifold)

Get the metric ``g`` of the manifold `M`.
"""
metric(::MetricManifold)

function metric(M::MetricManifold)
    return M.metric
end

"""
    ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`.
"""
ricci_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_curvature(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Ric = ricci_tensor(M, p, B; backend=backend)
    S = sum(Ginv .* Ric)
    return S
end
@decorator_transparent_signature ricci_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(
    ::typeof(ricci_curvature),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end
"""
    ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Ricci tensor, also known as the Ricci curvature tensor,
of the manifold `M` at the point `p` using basis `B`,
see [https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition](https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition).
"""
ricci_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    R = riemann_tensor(M, p, B; kwargs...)
    n = size(R, 1)
    Ric = allocate(R, Size(n, n))
    @einsum Ric[i, j] = R[l, i, l, j]
    return Ric
end
@decorator_transparent_signature ricci_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(::typeof(ricci_tensor), ::MetricManifold, args...)
    return Val(:parent)
end
@doc raw"""
    riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature
tensor, at the point `p`. The dimensions of the resulting multi-dimensional
array are ordered ``(l,i,j,k)``.
"""
riemann_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function riemann_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    Γ = christoffel_symbols_second(M, p, B; backend=backend)
    ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n
    R = allocate(∂Γ, Size(n, n, n, n))
    @einsum R[l, i, j, k] =
        ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k]
    return R
end
@decorator_transparent_signature riemann_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)
function decorator_transparent_dispatch(::typeof(riemann_tensor), ::MetricManifold, args...)
    return Val(:parent)
end

@doc raw"""
    sharp(N::MetricManifold{M,G}, p, ξ::FVector{CotangentSpaceType})

Compute the musical isomorphism to transform the cotangent vector `ξ` from the
[`AbstractManifold`](@ref) `M` equipped with [`AbstractMetric`](@ref) `G` to a tangent by
computing

````math
ξ^♯ = G_p^{-1} ξ,
````
where ``G_p`` is the local matrix representation of `G`, i.e. one employs
[`inverse_local_metric`](@ref) here to obtain ``G_p^{-1}``.
"""
sharp(::MetricManifold, ::Any, ::CoTFVector)

function sharp!(M::N, X::TFVector, p, ξ::CoTFVector) where {N<:MetricManifold}
    Ginv = inverse_local_metric(M, p, X.basis)
    copyto!(X.data, Ginv * ξ.data)
    return X
end

function Base.show(io::IO, M::MetricManifold)
    return print(io, "MetricManifold($(M.manifold), $(M.metric))")
end

@doc raw"""
    solve_exp_ode(
        M::MetricManifold,
        p,
        X,
        tspan,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
        solver = AutoVern9(Rodas5()),
        kwargs...,
    )

Approximate the exponential map on the manifold over the provided timespan
assuming the Levi-Civita connection by solving the ordinary differential
equation

```math
\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,
```

where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and
the Einstein summation convention is assumed. The arguments `tspan` and
`solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword
arguments that will be passed to `OrdinaryDiffEq.solve`.

Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.

!!! note
    This function only works when
    [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) is loaded with
    ```julia
    using OrdinaryDiffEq
    ```
"""
function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...)
    return error(
        "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.",
    )
end


text:
module RvSpectraKitLearn

using Compat
import Compat.view

export fit_gen_pca_rv, est_rvs_from_pc
include("alt_pca.jl")

export fit_pca_default, fit_pca_eford_proj_only, fit_pca_eford, test_pca_eford
include("generalized_pca.jl")

export calc_doppler_component_simple, calc_doppler_quadratic_term_simple
include("deriv_spectra_simple.jl")
export calc_doppler_component_gp
include("deriv_spectra_gp.jl")

export fit_arp_slr_2d, refit_arp_slr_2d
include("alg_regression_path/sparse_linear_regression.jl")
export print_alr_path, compute_alr_path_stats, compute_alr_path_chisq!
include("alg_regression_path/util.jl")

#export fit_arp_rrr
#include("alg_regression_path/reduced_rank_multitask_learning.jl")  # Implemented, but untested

export get_filenames, read_filelist, make_noisy_spectrum, make_noisy_spectra
include("read_inputs.jl")
#include("rand_junk.jl")

end # module

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn