#ifndef PM_SUBTRACTION_KERNELS_H
#define PM_SUBTRACTION_KERNELS_H

#include <pslib.h>

/// Type of subtraction kernel
typedef enum {
    PM_SUBTRACTION_KERNEL_NONE,         ///< Nothing --- an error
    PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
    PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
} pmSubtractionKernelsType;

/// Modes --- specifies which image to convolve
typedef enum {
    PM_SUBTRACTION_MODE_ERR,            // Error in the mode
    PM_SUBTRACTION_MODE_1,              // Convolve image 1
    PM_SUBTRACTION_MODE_2,              // Convolve image 2
    PM_SUBTRACTION_MODE_UNSURE,         // Not sure yet which image to convolve so try to satisfy both
    PM_SUBTRACTION_MODE_DUAL,           // Dual convolution
} pmSubtractionMode;

/// Kernels specification
typedef struct {
    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    psString description;               ///< Description of the kernel parameters
    long num;                           ///< Number of kernel components (not including the spatial ones)
    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
    psVector *widths;                   ///< Gaussian FWHMs (ISIS)
    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
    float penalty;                      ///< Penalty for wideness
    psVector *penalties;                ///< Penalty for each kernel component
    int size;                           ///< The half-size of the kernel
    int inner;                          ///< The size of an inner region
    int spatialOrder;                   ///< The spatial order of the kernels
    int bgOrder;                        ///< The order for the background fitting
    pmSubtractionMode mode;             ///< Mode for subtraction
    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    psVector *solution1, *solution2;    ///< Solution for the PSF matching
} pmSubtractionKernels;

// Assertion to check pmSubtractionKernels
#define PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(KERNELS, RETURNVALUE) { \
    PS_ASSERT_PTR_NON_NULL(KERNELS, RETURNVALUE); \
    PS_ASSERT_STRING_NON_EMPTY((KERNELS)->description, RETURNVALUE); \
    PS_ASSERT_INT_POSITIVE((KERNELS)->num, RETURNVALUE); \
    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->u, RETURNVALUE); \
    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->v, RETURNVALUE); \
    PS_ASSERT_VECTOR_TYPE((KERNELS)->u, PS_TYPE_S32, RETURNVALUE); \
    PS_ASSERT_VECTOR_TYPE((KERNELS)->v, PS_TYPE_S32, RETURNVALUE); \
    PS_ASSERT_VECTOR_SIZE((KERNELS)->u, (KERNELS)->num, RETURNVALUE); \
    PS_ASSERT_VECTOR_SIZE((KERNELS)->v, (KERNELS)->num, RETURNVALUE); \
    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS) { \
        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
    } \
    if ((KERNELS)->uStop || (KERNELS)->vStop) { \
        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \
        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->vStop, RETURNVALUE); \
        PS_ASSERT_VECTOR_TYPE((KERNELS)->uStop, PS_TYPE_S32, RETURNVALUE); \
        PS_ASSERT_VECTOR_TYPE((KERNELS)->vStop, PS_TYPE_S32, RETURNVALUE); \
        PS_ASSERT_VECTOR_SIZE((KERNELS)->uStop, (KERNELS)->num, RETURNVALUE); \
        PS_ASSERT_VECTOR_SIZE((KERNELS)->vStop, (KERNELS)->num, RETURNVALUE); \
    } \
    if ((KERNELS)->preCalc) { \
        PS_ASSERT_ARRAY_NON_NULL((KERNELS)->preCalc, RETURNVALUE); \
        PS_ASSERT_ARRAY_SIZE((KERNELS)->preCalc, (KERNELS)->num, RETURNVALUE); \
    } \
    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->size, RETURNVALUE); \
    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->inner, RETURNVALUE); \
    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->spatialOrder, RETURNVALUE); \
    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->bgOrder, RETURNVALUE); \
}

// Assertion to check that the solution is attached
#define PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(KERNELS, RETURNVALUE) { \
    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->solution1, RETURNVALUE); \
    PS_ASSERT_VECTOR_TYPE((KERNELS)->solution1, PS_TYPE_F64, RETURNVALUE); \
    PS_ASSERT_VECTOR_SIZE((KERNELS)->solution1, \
                          (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder) + 1 + \
                              PM_SUBTRACTION_POLYTERMS((KERNELS)->bgOrder), \
                          RETURNVALUE); \
    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { \
        PS_ASSERT_VECTOR_NON_NULL(kernels->solution2, RETURNVALUE); \
        PS_ASSERT_VECTOR_TYPE((KERNELS)->solution2, PS_TYPE_F64, RETURNVALUE); \
        PS_ASSERT_VECTOR_SIZE((KERNELS)->solution2, \
                              (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder), \
                               RETURNVALUE); \
    } \
}

/// Generate a delta-function grid for subtraction kernels (like the POIS kernel)
bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to
                                   int start, ///< Index at which to start appending
                                   int size ///< Half-size of the grid
    );

/// General allocator for pmSubtractionKernels
///
/// Unlike the functions for the specific kernel type, this function does not set up the basis functions, but
/// merely allocates space for their storage.
pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, ///< Number of basis functions
                                                pmSubtractionKernelsType type, ///< Kernel type
                                                int size, ///< Half-size of kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate POIS kernels
pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
                                               int spatialOrder, ///< Order of spatial variations
                                               float penalty, ///< Penalty for wideness
                                               pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate ISIS kernels without the flux scaling built in
pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, ///< Half-size of the kernel
                                                    int spatialOrder, ///< Order of spatial variations
                                                    const psVector *fwhms, ///< Gaussian FWHMs
                                                    const psVector *orders, ///< Polynomial order of gaussians
                                                    float penalty, ///< Penalty for wideness
                                                    pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate ISIS kernels
pmSubtractionKernels *pmSubtractionKernelsISIS(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               const psVector *fwhms, ///< Gaussian FWHMs
                                               const psVector *orders, ///< Polynomial order of gaussians
                                               float penalty, ///< Penalty for wideness
                                               pmSubtractionMode mode ///< Mode for subtraction
                                               );

/// Generate SPAM kernels
pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               int inner, ///< Inner radius to preserve unbinned
                                               int binning, ///< Kernel binning factor
                                               float penalty, ///< Penalty for wideness
                                               pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate FRIES kernels
pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner, ///< Inner radius to preserve unbinned
                                                float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate GUNK kernels
pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               const psVector *fwhms, ///< Gaussian FWHMs
                                               const psVector *orders, ///< Polynomial order of gaussians
                                               int inner, ///< Inner radius containing grid of delta functions
                                               float penalty, ///< Penalty for wideness
                                               pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate RINGS kernels
pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, ///< Half-size of the kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner, ///< Inner radius to preserve unbinned
                                                int ringsOrder, ///< Polynomial order
                                                float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
    );


/// Generate a kernel of a specified type
pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type
                                                   int size, ///< Half-size of the kernel
                                                   int spatialOrder, ///< Order of spatial variations
                                                   const psVector *fwhms, ///< Gaussian FWHMs
                                                   const psVector *orders, ///< Polynomial order of gaussians
                                                   int inner, ///< Inner radius to preserve unbinned
                                                   int binning, ///< Kernel binning factor
                                                   int ringsOrder, ///< Polynomial order for RINGS
                                                   float penalty, ///< Penalty for wideness
                                                   pmSubtractionMode mode ///< Mode for subtraction
    );

/// Generate a kernel using the description
pmSubtractionKernels *pmSubtractionKernelsFromDescription(
    const char *description,            ///< Description of kernel
    int bgOrder,                        ///< Polynomial order for background fitting
    pmSubtractionMode mode              ///< Mode for subtraction
    );

/// Return the appropriate type from a string
pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *string // String name for kernel type
    );


#endif
