#ifndef LSHDERIVS_H
#define LSHDERIVS_H

#include "SETTINGS.h"

// generates a matrix [v x] where [v x] * u = v.cross(u)
MATRIX3 asCross(const VEC3& v);

// Gets a simple lsh derivative with respect to omega (D al/ D om)
MATRIX3 DdarbEqnDom(const VEC3& globalDarb, const VEC3& base, const VEC3& t0, const double& s);
// Gets a derivative with respect to local darb coordinates (uses DdarbEqnDom and chain rule)
// (D al / D q)
MATRIX3 Dlsh(const VEC3& lDarb, const VEC3& base, const MATRIX3& F0, const double& s);
// Gets derivative of lsh with respect to t0 (D al / D t0)
MATRIX3 DdarbEqnDt(const VEC3& globalDarb, const VEC3& t0, const double& s);

// Gets a simple (D al / D s) derivative with respect to omega (D2 al / D s D om)
MATRIX3 DFdarbEqnDom(const VEC3& globalDarb, const VEC3& t0, const double& s);
// Gets a simple (D al / D s) derivative with respect to local coords (D2 al / D s D q)
MATRIX3 DFlsh(const VEC3& localDarb, const MATRIX3& F0, const double& s, const int col = 0);
// Gets derivative of (D al / D s) with respect to t0 (D2 al / D s D t0)
MATRIX3 DFdarbEqnDt(const VEC3& globalDarb, const VEC3& t0, const double& s);

/**
 * Fills the given frames std::vector<> from lsh sampling
 * \param lshs the set of lshs 
 * \param frames the std::vector to fill the frames into
 * \param F0 the first frame
 */
void fillFrames(const std::vector<superHelix>& lshs, std::vector<MATRIX3>& frames, const MATRIX3& F0);

/**
 * Fills the given roots and frames std::vector<>s with roots and frames from lsh sampling
 * \param lshs the set of lshs 
 * \param roots the std::vector to fill the root points into
 * \param frames the std::vector to fill the frames into
 * \param p0 the first root point
 * \param F0 the first frame
 */
void fillRootsAndFrames(const std::vector<superHelix>& lshs, std::vector<VEC3>& roots, std::vector<MATRIX3>& frames, const VEC3& p0, const MATRIX3& F0);

/**
 * Recursively gets the 3x3 Jacobian (D2 alpha_i/ D s D q_j) at length s (relative to alpha_i's start point)
 * \param lshs the set of lshs
 * \param F0 the initial frame
 * \param s the length sample (relative to alpha_i's starting point)
 * \param i the index of which alpha_i to use
 * \param j which qj darboux are we applying a variation to?
 * \param col which column of F0 is getting transported?
 */
MATRIX3 D2aliDsDqj(const std::vector<superHelix>& lshs, const MATRIX3& F0, const double s, const int i, const int j, const int col);


/**
 * Recursively gets the 3x3 Jacobian (D alpha_i / D q_j) at length s (relative to alpha_i's start point)
 * \param lshs the set of lshs
 * \param p0 the initial point
 * \param F0 the initial frame
 * \param s the length sample (relative to alpha_i's starting point)
 * \param i the index of which alpha_i to use
 * \param j which qj darboux are we applying a variation to?
 */
MATRIX3 DaliDqj(const std::vector<superHelix>& lshs, const VEC3& p0, const MATRIX3& F0, const double s, const int i, const int j);

/**
 * Makes and fills helper tangent matrix up to the ith 3x3 block row 
 * \param lshs the current set of lshs
 * \param frames the set of all frames needed to perform the computation
 * \param i the index of the superhelix that the holders are filling to (inclusive)
 * \param tHolder, uHolder, vHolder the matrices that will contain the t, u, v information
 */
void MakeD2aliDsDqArrays(const std::vector<superHelix>& lshs, const std::vector<MATRIX3>& frames, const int i,
                          MATRIX& tHolder, MATRIX& uHolder, MATRIX& vHolder);

/**
 * Makes and fills helper tangent matrix up to the ith 3x3 block row 
 * \param lshs the current set of lshs
 * \param roots the set of all roots needed to perform the computation
 * \param frames the set of all frames needed to perform the computation
 * \param i the index of the superhelix that the holders are filling up to (inclusive)
 * \param tHolder, uHolder, vHolder the matrices that contain the t, u, v information (expected to already be filled)
 * \param difHolder the matrix that will contain all the endpoint derivative information
 */
void MakeDaliDqArrays(const std::vector<superHelix>& lshs, const std::vector<VEC3>& roots, const std::vector<MATRIX3>& frames, const int i, 
                      const MATRIX& tHolder, const MATRIX& uHolder, const MATRIX& vHolder, MATRIX& difHolder);

/**
 * Gets the big 3 x 3 * lshs.size() Jacobian of alpha_i at length s (optimizer must manage which length to land s in)
 * Uses memoization to not have constant recursion
 * \param lshs the current set of lshs 
 * \param p0 the initial point
 * \param F0 the initial frame
 * \param s the length sample (relative to alpha_i's starting point)
 * \param i the index of which alpha_i to use
 */
MATRIX fullConnectedStateDeriv(const std::vector<superHelix>& lshs, const VEC3& p0, const MATRIX3& F0, const double s, const int i);

/**
 * Gets multiple big 3 x 3 * lshs.size() Jacobian of strand at lengths s
 * Uses memoization to not have constant recursion
 * \param lshs the current set of lshs 
 * \param p0 the initial point
 * \param F0 the initial frame
 * \param sAs the length samples (relative to the full strand with domain rescaled to [0, 1])
 */
std::vector<MATRIX> fullConnectedStateDerivs(const std::vector<superHelix>& lshs, const VEC3& p0, const MATRIX3& F0, const std::vector<double>& sAs);

#endif