//+------------------------------------------------------------------+
//|                                                        slspq.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "num_diff.mqh"
#include<Object.mqh>
//---
#define SLSQP_EPS 2.220446049250313e-16
//---
//+------------------------------------------------------------------+
//| copy from vector to array                                        |
//+------------------------------------------------------------------+
template<typename T>
bool Copy(T& dest_array[],vector<T>& src_vec,int count,int dest_start = 0, int src_start = 0)
  {
// Bounds verification checking if requested copy counts exceed dimensions of destination or source
   if(int(dest_array.Size())<fabs(count+dest_start) || int(src_vec.Size())<fabs(count+src_start))
     {
      Print(__FUNCTION__, " invalid inputs ");
      return false;
     }
//---
// Element-by-element buffer assignment loop
   for(int i = 0; i<(count); ++i)
      dest_array[dest_start+i] = src_vec[src_start+i];
//---
   return true;
  }
//+------------------------------------------------------------------+
//| copy from array to vector                                        |
//+------------------------------------------------------------------+
template<typename T>
bool Copy(vector<T>& dest_vec, T& src_array[],int count,int dest_start = 0, int src_start = 0)
  {
// Bounds verification matching dynamic vector tracking against plain reference buffer indices
   if(dest_vec.Size()<ulong(count+dest_start) || int(src_array.Size())<fabs(count+src_start))
     {
      Print(__FUNCTION__, " invalid inputs ");
      return false;
     }
//---
// Element-by-element vector insertion loop
   for(int i = 0; i<(count); ++i)
      dest_vec[dest_start+i] = src_array[src_start+i];
//---
   return true;
  }

//+------------------------------------------------------------------+
//| constraint                                                       |
//+------------------------------------------------------------------+
// Structure defining non-linear constraints for the SLSQP optimizer,
// handling both single-valued (scalar) and vector-valued constraints.
struct slsqp_constraint
  {
   int               m;      // Number of constraint functions (dimension of the constraint vector)
   IObjective*       f_data; // Pointer to the underlying interface tracking the user's objective/constraint logic

   // Default Constructor
                     slsqp_constraint(void)
     {
      m=0;
     }

   // Copy Constructor
                     slsqp_constraint(slsqp_constraint& other)
     {
      m = other.m;
      f_data = other.f_data;
     }

   // Destructor
                    ~slsqp_constraint(void)
     {
     }

   // Assignment Operator Overload
   void              operator=(slsqp_constraint& other)
     {
      m = other.m;
      f_data = other.f_data;
     }

   // Evaluates a single scalar constraint value and its optional gradients
   double            f(int n, double& x[],double& gradient[],IObjective* func_data, int &offset)
     {
      vector vx = vector::Zeros(n); // Initialize MQL5 native vector container matching decision parameters
      vector out;
      Copy(vx,x,n);                 // Convert and pass the primitive array values into the native vector

      // If the optimizer expects gradient evaluation
      if(gradient.Size())
        {
         // Simultaneously calculate the constraint value and its Jacobian/Gradient profile
         ObjReturn or = func_data.fun_and_grad(vx);

         // Flat-map the matrix gradients into the 1D gradient buffer expected by the solver core
         for(int j = 0; j<n; ++j)
            for(int i = 0; i<1; ++i)
               gradient[offset + (i+j*1)] = or.mg[i,j];
         out = or.mf;
        }
      else
         // Value-only evaluation mode without executing heavy gradient tracks
         out = func_data.objective_function(vx);
      ++offset;
      return out[0]; // Return the singular scalar constraint function result
     }

   // Evaluates vector-valued multi-constraints (m > 1) with their multi-dimensional Jacobians
   void              mf(int mm, double& result[], int n, double& x[],double& gradient[],IObjective* func_data, int &offset)
     {
      vector vx = vector::Zeros(n);
      Copy(vx,x,n); // Use the global helper function to move array inputs into vector format

      // Joint vector evaluation block mapping multiple outputs and multi-gradient dimensions
      if(gradient.Size())
        {
         ObjReturn or = func_data.fun_and_grad(vx);
         Copy(result,or.mf,mm,offset);
         for(int j = 0; j<n; ++j)
           {
            // Map the constraint Jacobian matrix rows into the flattened vector expected by the optimization loop
            for(int i = 0; i<mm; ++i)
               gradient[offset + (i+j*mm)] = or.mg[i,j];
           }
        }
      else
        {
         // Multi-constraint value evaluation without Jacobian configurations
         vector res = func_data.objective_function(vx);
         Copy(result,res,mm,offset); // Export local vector allocations back to the flat output array at the designated offset
        }
      offset+=mm;
     }
  };

//+------------------------------------------------------------------+
//| count the number of constraints                                  |
//+------------------------------------------------------------------+
// Iterates across multiple constraint elements to compute the aggregate constraints dimension sizing
int                slsqp_count_constraints(int p, slsqp_constraint& c[])
  {
   int i, count = 0;
   for(i = 0; i < p; ++i)
      count += c[i].m;
   return count;
  }

//+------------------------------------------------------------------+
//| constraints parsing utility                                      |
//+------------------------------------------------------------------+
// Scans multi-constraint setups to locate the largest single vector dimension block
int                slsqp_max_constraint_dim(int p, slsqp_constraint& c[])
  {
   int i, max_dim = 0;
   for(i = 0; i < p; ++i)
      if(c[i].m > max_dim)
         max_dim = c[i].m;
   return max_dim;
  }

//+------------------------------------------------------------------+
//| evaluates a or all constraints                                   |
//+------------------------------------------------------------------+
// Directs single scalar vs multi-dimensional constraint vector metrics reporting to targeted workspace targets
void              slsqp_eval_constraint(double& result[], double& grad[], slsqp_constraint& c, int n, double& x[], int offset)
  {
   if(c.m == 1)
      result[offset] = c.f(n, x, grad, c.f_data,offset); // Direct scalar assignment mapping
   else
      c.mf(c.m, result, n, x, grad, c.f_data,offset); // Structured array segment expansion pass
  }
//+------------------------------------------------------------------+
//| Small utilities                                                  |
//+------------------------------------------------------------------+
double SLSQP_NaN()
  {
   double zero = 0.0;
   return(MathSqrt(-1.0 + zero)); // NaN sentinel, used exactly like SciPy's convention
// for "this bound is absent" (xl[i]/xu[i] = NaN).
  }

//+------------------------------------------------------------------+
//| Is not a number wrapper                                          |
//+------------------------------------------------------------------+
bool SLSQP_IsNaN(const double v)
  {
   return(!MathIsValidNumber(v));
  }

//+------------------------------------------------------------------+
//| Apply sign change to input                                       |
//+------------------------------------------------------------------+
double SLSQP_CopySign(const double mag, const double sgn)
  {
   double a = MathAbs(mag);
   return(sgn < 0.0 ? -a : a);
  }

//+------------------------------------------------------------------+
//| Hypotenuse computation                                           |
//+------------------------------------------------------------------+
double SLSQP_Hypot(const double a, const double b)
  {
   double aa = MathAbs(a), ab = MathAbs(b);
   if(aa == 0.0 && ab == 0.0)
      return(0.0);
   if(aa > ab)
      return(aa*MathSqrt(1.0 + (ab/aa)*(ab/aa)));
   return(ab*MathSqrt(1.0 + (aa/ab)*(aa/ab)));
  }

//+------------------------------------------------------------------+
//| Level-1 / Level-2 "BLAS" primitives                              |
//| Every array-taking function below uses an (array[], offset)      |
//| pair in place of a C pointer: arr[off+i] instead of ptr[i].      |
//| These are close, deliberately unoptimized transliterations of    |
//| the reference BLAS algorithms -- validated in C against a full   |
//| SciPy SLSQP fuzz suite before being ported here.                 |
//+------------------------------------------------------------------+
void SLSQP_daxpy(const int n, const double alpha,
                 const double &x[], const int xo, const int incx,
                 double &y[], const int yo, const int incy)
  {
   for(int i = 0; i < n; i++)
      y[yo + i*incy] += alpha*x[xo + i*incx];
  }

//+------------------------------------------------------------------+
//|scales vector X(n) by constant da                                 |
//+------------------------------------------------------------------+
void SLSQP_dscal(const int n, const double alpha, double &x[], const int xo, const int incx)
  {
   for(int i = 0; i < n; i++)
      x[xo + i*incx] *= alpha;
  }

//+------------------------------------------------------------------+
//|compute the L2 norm of array DX of length N, stride INCX          |
//+------------------------------------------------------------------+
double SLSQP_dnrm2(const int n, const double &x[], const int xo, const int incx)
  {
   double s = 0.0;
   for(int i = 0; i < n; i++)
     {
      double v = x[xo + i*incx];
      s += v*v;
     }
   return(MathSqrt(s));
  }

//+------------------------------------------------------------------+
//| dot product dx dot dy                                            |
//+------------------------------------------------------------------+
double SLSQP_ddot(const int n, const double &x[], const int xo, const int incx,
                  const double &y[], const int yo, const int incy)
  {
   double s = 0.0;
   for(int i = 0; i < n; i++)
      s += x[xo + i*incx]*y[yo + i*incy];
   return(s);
  }
//+------------------------------------------------------------------------+
//| y = alpha*op(A)*x + beta*y ; A is m x n, column-major, leading dim lda.|
//| trans=false -> op(A)=A (y has length m, x has length n)                |
//| trans=true  -> op(A)=A^T (y has length n, x has length m)              |
//+------------------------------------------------------------------------+
void SLSQP_dgemv(const bool trans, const int m, const int n, const double alpha,
                 const double &a[], const int ao, const int lda,
                 const double &x[], const int xo, const int incx, const double beta,
                 double &y[], const int yo, const int incy)
  {
   if(!trans)
     {
      for(int i = 0; i < m; i++)
        {
         double s = 0.0;
         for(int j = 0; j < n; j++)
            s += a[ao + i + j*lda]*x[xo + j*incx];
         y[yo + i*incy] = beta*y[yo + i*incy] + alpha*s;
        }
     }
   else
     {
      for(int j = 0; j < n; j++)
        {
         double s = 0.0;
         for(int i = 0; i < m; i++)
            s += a[ao + i + j*lda]*x[xo + i*incx];
         y[yo + j*incy] = beta*y[yo + j*incy] + alpha*s;
        }
     }
  }

//+-----------------------------------------------------------------------------------------+
//|0-based packed lower-triangular column index (i>=j), matching LAPACK's packed 'L' layout.|
//+-----------------------------------------------------------------------------------------+
int SLSQP_LtpIdx(const int i, const int j, const int n)
  {
   return(j*n - (j*(j-1))/2 + (i-j));
  }


//+-----------------------------------------------------------------------------------------+
//|x := op(L)*x, L unit lower-triangular, packed. trans=false -> L*x, trans=true -> L^T*x.  |
//| (Only uplo='L', diag='U' are ever needed by SLSQP.)                                     |
//+-----------------------------------------------------------------------------------------+
void SLSQP_dtpmv(const bool trans, const int n, const double &ap[], const int apo,
                 double &x[], const int xo, const int incx)
  {
   if(!trans)
     {
      for(int j = n-1; j >= 0; j--)
        {
         double temp = x[xo + j*incx];
         for(int i = n-1; i > j; i--)
            x[xo + i*incx] += temp*ap[apo + SLSQP_LtpIdx(i,j,n)];
        }
     }
   else
     {
      for(int j = 0; j < n; j++)
        {
         double temp = x[xo + j*incx];
         for(int i = j+1; i < n; i++)
            temp += ap[apo + SLSQP_LtpIdx(i,j,n)]*x[xo + i*incx];
         x[xo + j*incx] = temp;
        }
     }
  }


//+------------------------------------------------------------------------------+
//|Solve L*x = b in place (unit lower-triangular, packed, forward substitution). |
//|(Only uplo='L', trans='N', diag='U' is ever needed by SLSQP.)                 |
//+------------------------------------------------------------------------------+
void SLSQP_dtpsv(const int n, const double &ap[], const int apo, double &x[], const int xo, const int incx)
  {
   for(int j = 0; j < n; j++)
     {
      double xj = x[xo + j*incx];
      if(xj != 0.0)
        {
         for(int i = j+1; i < n; i++)
            x[xo + i*incx] -= xj*ap[apo + SLSQP_LtpIdx(i,j,n)];
        }
     }
  }


//+---------------------------------------------------------------------------------------+
//|Solve op(A)*x = b in place, A upper-triangular n x n, non-unit diag.                   |
//|trans=false -> A*x=b (back substitution), trans=true -> A^T*x=b (forward substitution).|
//+---------------------------------------------------------------------------------------+
void SLSQP_dtrsv(const bool trans, const int n, const double &a[], const int ao, const int lda,
                 double &x[], const int xo, const int incx)
  {
   if(!trans)
     {
      for(int j = n-1; j >= 0; j--)
        {
         if(x[xo + j*incx] != 0.0)
           {
            x[xo + j*incx] /= a[ao + j + j*lda];
            double xj = x[xo + j*incx];
            for(int i = 0; i < j; i++)
               x[xo + i*incx] -= xj*a[ao + i + j*lda];
           }
        }
     }
   else
     {
      for(int j = 0; j < n; j++)
        {
         double s = x[xo + j*incx];
         for(int i = 0; i < j; i++)
            s -= a[ao + i + j*lda]*x[xo + i*incx];
         x[xo + j*incx] = s / a[ao + j + j*lda];
        }
     }
  }

//+---------------------------------------------------------------------------------------------------+
//|Solve X*A = alpha*B for X, A n x n upper-triangular (lda), B is m x n (ldb=m), overwritten with X. |
//| This is the one specific pattern SLSQP's lsi() needs (BLAS dtrsm("R","U","N","N",...)).           |
//+---------------------------------------------------------------------------------------------------+
void SLSQP_dtrsm_RUNN(const int m, const int n, const double alpha,
                      const double &a[], const int ao, const int lda,
                      double &b[], const int bo, const int ldb)
  {
   for(int i = 0; i < m; i++)
      for(int j = 0; j < n; j++)
         b[bo + i + j*ldb] *= alpha;
   for(int j = 0; j < n; j++)
     {
      for(int k = 0; k < j; k++)
        {
         double akj = a[ao + k + j*lda];
         if(akj != 0.0)
           {
            for(int i = 0; i < m; i++)
               b[bo + i + j*ldb] -= b[bo + i + k*ldb]*akj;
           }
        }
      double ajj = a[ao + j + j*lda];
      for(int i = 0; i < m; i++)
         b[bo + i + j*ldb] /= ajj;
     }
  }

//+------------------------------------------------------------------+
//| Householder / Givens helpers                                     |
//+------------------------------------------------------------------+

// Standard Householder reflector generation: reflects (alpha,x) -> (beta,0,...,0),
// beta = -sign(alpha)*hypot(alpha,||x||). Overwrites alpha_arr[alpha_idx] with beta
// and x[] with the reflector tail; returns tau via tau_arr[tau_idx].
void SLSQP_dlarfg_std(const int n, double &alpha_arr[], const int alpha_idx,
                      double &x[], const int xo, const int incx,
                      double &tau_arr[], const int tau_idx)
  {
   if(n <= 1)
     {
      tau_arr[tau_idx] = 0.0;
      return;
     }
   int nm1 = n-1;
   double xnorm = SLSQP_dnrm2(nm1, x, xo, incx);
   if(xnorm == 0.0)
     {
      tau_arr[tau_idx] = 0.0;
      return;
     }
   double a0 = alpha_arr[alpha_idx];
   double beta = -SLSQP_CopySign(SLSQP_Hypot(a0, xnorm), a0);
   tau_arr[tau_idx] = (beta - a0)/beta;
   double scale = 1.0/(a0 - beta);
   for(int i = 0; i < nm1; i++)
      x[xo + i*incx] *= scale;
   alpha_arr[alpha_idx] = beta;
  }

//+------------------------------------------------------------------------+
//|"Positive" variant: like dlarfg_std but forces beta >= 0. Used by NNLS. |
//+------------------------------------------------------------------------+
void SLSQP_dlarfgp(const int n, double &alpha_arr[], const int alpha_idx,
                   double &x[], const int xo, const int incx,
                   double &tau_arr[], const int tau_idx)
  {
   if(n <= 1)
     {
      tau_arr[tau_idx] = 0.0;
      return;
     }
   int nm1 = n-1;
   double xnorm = SLSQP_dnrm2(nm1, x, xo, incx);
   double a0 = alpha_arr[alpha_idx];
   if(xnorm == 0.0 && a0 >= 0.0)
     {
      tau_arr[tau_idx] = 0.0;
      return;
     }
   double beta = SLSQP_Hypot(a0, xnorm); // always >= 0
   tau_arr[tau_idx] = (beta - a0)/beta;
   double abdiff = (a0 - beta);
   double scale = abdiff?1.0/abdiff:double("inf");
   for(int i = 0; i < nm1; i++)
      x[xo + i*incx] *= scale;
   alpha_arr[alpha_idx] = beta;
  }

//+-------------------------------------------------------------------------------+
//|"Positive" Givens rotation generator: [cs sn;-sn cs]*[f;g] = [r;0], r >= 0.    |
//| f, g are plain inputs (already-read scalar values); the caller is responsible |
//| for writing r back into wherever f came from and zeroing wherever g came from |
//| (mirroring exactly how nnls.c uses BLAS's dlartgp).                           |
//+-------------------------------------------------------------------------------+
void SLSQP_dlartgp(const double f, const double g, double &cs, double &sn, double &r)
  {
   if(g == 0.0)
     {
      cs = 1.0;
      sn = 0.0;
      r = f;
      return;
     }
   if(f == 0.0)
     {
      cs = 0.0;
      sn = 1.0;
      r = g;
      return;
     }
   double rr = SLSQP_Hypot(f,g);
   r = rr;
   cs = f/rr;
   sn = g/rr;
  }

//+----------------------------------------------------------------------------------------+
//|Apply a Householder reflector H = I - tau*v*v^T to C (m x n, leading dim ldc).          |
//| side_left=true: C := H*C (v has length m). side_left=false: C := C*H (v has length n). |
//| v[vo] (with stride incv) must already hold its own explicit leading "1" element        |
//| (the caller temporarily pokes a 1.0 into the matrix position that represents it --     |
//| this mirrors exactly how the reference LAPACK routines use DLARF internally).          |
//+----------------------------------------------------------------------------------------+
void SLSQP_dlarf(const bool side_left, const int m, const int n,
                 const double &v[], const int vo, const int incv, const double tau,
                 double &c[], const int co, const int ldc,
                 double &work[], const int wo)
  {
   if(tau == 0.0)
      return;
   if(side_left)
     {
      for(int j = 0; j < n; j++)
        {
         double s = 0.0;
         for(int i = 0; i < m; i++)
            s += c[co + i + j*ldc]*v[vo + i*incv];
         work[wo + j] = s;
        }
      for(int j = 0; j < n; j++)
        {
         double wj = tau*work[wo + j];
         for(int i = 0; i < m; i++)
            c[co + i + j*ldc] -= v[vo + i*incv]*wj;
        }
     }
   else
     {
      for(int i = 0; i < m; i++)
        {
         double s = 0.0;
         for(int j = 0; j < n; j++)
            s += c[co + i + j*ldc]*v[vo + j*incv];
         work[wo + i] = s;
        }
      for(int i = 0; i < m; i++)
        {
         double wi = tau*work[wo + i];
         for(int j = 0; j < n; j++)
            c[co + i + j*ldc] -= wi*v[vo + j*incv];
        }
     }
  }

//+------------------------------------------------------------------+
//| Unblocked Householder QR (dgeqr2) and applying Q/Q^T (dorm2r)    |
//+------------------------------------------------------------------+

// QR-factorize A (m x n, leading dim lda) in place: R in the upper triangle
// (including diagonal), reflectors stored below the diagonal, tau[0..min(m,n)-1].
void SLSQP_dgeqr2(const int m, const int n, double &a[], const int ao, const int lda,
                  double &tau[], const int tauo, double &work[], const int wo)
  {
   int k = MathMin(m,n);
   for(int i = 0; i < k; i++)
     {
      int len = m - i;
      SLSQP_dlarfg_std(len, a, ao + i + i*lda, a, ao + (i+1) + i*lda, 1, tau, tauo + i);
      if(i < n-1)
        {
         double aii = a[ao + i + i*lda];
         a[ao + i + i*lda] = 1.0;
         int mm = m - i, nn = n - i - 1;
         SLSQP_dlarf(true, mm, nn, a, ao + i + i*lda, 1, tau[tauo+i], a, ao + i + (i+1)*lda, lda, work, wo);
         a[ao + i + i*lda] = aii;
        }
     }
  }

//+--------------------------------------------------------------------------------------------+
//|Apply Q or Q^T (from dgeqr2, side='L' only -- the only case SLSQP needs) to C (m x n, ldc). |
//|trans=true -> C := Q^T*C ; trans=false -> C := Q*C. k = number of reflectors used.          |
//+--------------------------------------------------------------------------------------------+
void SLSQP_dorm2r(const bool trans, const int m, const int n, const int k,
                  double &a[], const int ao, const int lda, double &tau[], const int tauo,
                  double &c[], const int co, const int ldc, double &work[], const int wo)
  {
   int start, stop, step;
   if(trans)
     {
      start = 0;   // Q^T*C = H(k-1)...H(0)*C -> apply ascending
      stop = k;
      step = 1;
     }
   else
     {
      start = k-1;   // Q*C   = H(0)...H(k-1)*C -> apply descending
      stop = -1;
      step = -1;
     }
   for(int i = start; i != stop; i += step)
     {
      double aii = a[ao + i + i*lda];
      a[ao + i + i*lda] = 1.0;
      int mm = m - i, nn = n;
      SLSQP_dlarf(true, mm, nn, a, ao + i + i*lda, 1, tau[tauo+i], c, co + i, ldc, work, wo);
      a[ao + i + i*lda] = aii;
     }
  }

//+------------------------------------------------------------------+
//| Unblocked Householder RQ (dgerq2) and applying Q/Q^T (dormr2)    |
//| RQ decomposition places R as the LAST m columns of the m x n     |
//| input matrix (upper triangular there); Q is n x n, represented   |
//| implicitly as a product of m elementary reflectors.              |
//+------------------------------------------------------------------+
void SLSQP_dgerq2(const int m, const int n, double &a[], const int ao, const int lda,
                  double &tau[], const int tauo, double &work[], const int wo)
  {
   int k = MathMin(m,n);
   for(int idx = 0; idx < k; idx++)
     {
      int i = k - 1 - idx;
      int row = m - k + i;
      int pcol = n - k + i;              // pivot column (rightmost of the swept range)
      int len = pcol + 1;
      // x = row 'row', columns 0..pcol-1, stepping across columns => stride = lda
      SLSQP_dlarfg_std(len, a, ao + row + pcol*lda, a, ao + row + 0*lda, lda, tau, tauo + idx);
      if(row > 0)
        {
         double aii = a[ao + row + pcol*lda];
         a[ao + row + pcol*lda] = 1.0;
         int mm = row, nn = pcol + 1;
         SLSQP_dlarf(false, mm, nn, a, ao + row + 0*lda, lda, tau[tauo+idx], a, ao, lda, work, wo);
         a[ao + row + pcol*lda] = aii;
        }
     }
  }

//+-----------------------------------------------------------------------------------------------+
//|Apply Q (from dgerq2) to C, side in {L,R}, trans in {N,T}, using k of the reflectors.          |
//| LAPACK convention: the reflector-storage matrix 'a' is k-by-M when side=L, k-by-N when side=R |
//| (M,N = dimensions of C) -- i.e. the "true" column count of 'a' to use when locating each      |
//| reflector's pivot column depends on which side we're applying from.                           |
//+-----------------------------------------------------------------------------------------------+
void SLSQP_dormr2(const bool side_left, const bool trans, const int m, const int n, const int k,
                  double &a[], const int ao, const int lda, double &tau[], const int tauo,
                  double &c[], const int co, const int ldc, double &work[], const int wo)
  {
   int mrq = k; // dgerq2 is always called with its own "m" parameter == k in this codebase
   bool leftMultQT  = side_left && trans;
   bool rightMultQ  = (!side_left) && (!trans);
   bool useDescending = (leftMultQT || rightMultQ);
   int ncols_of_a = side_left ? m : n;   // LAPACK: A is k-by-M (side=L) or k-by-N (side=R)

   for(int pass = 0; pass < k; pass++)
     {
      int idx = useDescending ? (k-1-pass) : pass;
      int i = k - 1 - idx;
      int row = mrq - k + i;
      int pcol = ncols_of_a - k + i;
      int len = pcol + 1;

      double aii = a[ao + row + pcol*lda];
      a[ao + row + pcol*lda] = 1.0;
      if(side_left)
        {
         int mm = len, nn = n;
         SLSQP_dlarf(true, mm, nn, a, ao + row + 0*lda, lda, tau[tauo+idx], c, co, ldc, work, wo);
        }
      else
        {
         int mm = m, nn = len;
         SLSQP_dlarf(false, mm, nn, a, ao + row + 0*lda, lda, tau[tauo+idx], c, co, ldc, work, wo);
        }
      a[ao + row + pcol*lda] = aii;
     }
  }

//+------------------------------------------------------------------+
//| dgelsy: rank-revealing least squares via column-pivoted           |
//| Householder QR (Businger-Golub pivoting). Solves min|Ax-b| for x  |
//| (n-vector), A is m x n (m>=n expected). Basic (not minimum-norm)  |
//| solution when rank-deficient -- SLSQP only reaches this routine   |
//| when there are no inequality/bound constraints at all, and only   |
//| the rank itself (not the min-norm refinement) affects the mode    |
//| SLSQP reports in that corner case.                                |
//+------------------------------------------------------------------+
void SLSQP_dgelsy(const int m, const int n, double &a[], const int ao, const int lda,
                  double &b[], const int bo, const double rcond, int &rank)
  {
   int pcols[];
   ArrayResize(pcols, n);
   double colnorm[];
   ArrayResize(colnorm, n);
   int k = MathMin(m,n);
   double tau[];
   ArrayResize(tau, MathMax(k,1));

   for(int j = 0; j < n; j++)
     {
      pcols[j] = j;
      double s = 0.0;
      for(int i = 0; i < m; i++)
        {
         double v = a[ao+i+j*lda];
         s += v*v;
        }
      colnorm[j] = s;
     }

   double firstdiag = 0.0;
   for(int kk = 0; kk < k; kk++)
     {
      int piv = kk;
      double best = colnorm[kk];
      for(int j = kk+1; j < n; j++)
         if(colnorm[j] > best)
           {
            best = colnorm[j];
            piv = j;
           }
      if(piv != kk)
        {
         for(int i = 0; i < m; i++)
           {
            double t = a[ao+i+kk*lda];
            a[ao+i+kk*lda] = a[ao+i+piv*lda];
            a[ao+i+piv*lda] = t;
           }
         double t2 = colnorm[kk];
         colnorm[kk] = colnorm[piv];
         colnorm[piv] = t2;
         int ti = pcols[kk];
         pcols[kk] = pcols[piv];
         pcols[piv] = ti;
        }
      int len = m - kk;
      double work_scratch[];
      ArrayResize(work_scratch, MathMax(n,1));
      SLSQP_dlarfg_std(len, a, ao+kk+kk*lda, a, ao+(kk+1)+kk*lda, 1, tau, kk);
      if(kk < n-1)
        {
         double aii = a[ao+kk+kk*lda];
         a[ao+kk+kk*lda] = 1.0;
         int mm = m-kk, nn = n-kk-1;
         SLSQP_dlarf(true, mm, nn, a, ao+kk+kk*lda, 1, tau[kk], a, ao+kk+(kk+1)*lda, lda, work_scratch, 0);
         a[ao+kk+kk*lda] = aii;
        }
      for(int j = kk+1; j < n; j++)
        {
         double v = a[ao+kk+j*lda];
         colnorm[j] -= v*v;
         if(colnorm[j] < 0)
            colnorm[j] = 0;
        }
      if(kk == 0)
         firstdiag = MathAbs(a[ao]);
     }

   double thresh = rcond * (firstdiag > 0 ? firstdiag : 1.0);
   int krank = 0;
   for(int kk = 0; kk < k; kk++)
     {
      if(MathAbs(a[ao+kk+kk*lda]) > thresh)
         krank++;
      else
         break;
     }
   rank = krank;

// Apply Q^T to b
   double w1[1];
   for(int kk = 0; kk < k; kk++)
     {
      double aii = a[ao+kk+kk*lda];
      a[ao+kk+kk*lda] = 1.0;
      int mm = m-kk;
      SLSQP_dlarf(true, mm, 1, a, ao+kk+kk*lda, 1, tau[kk], b, bo+kk, m, w1, 0);
      a[ao+kk+kk*lda] = aii;
     }

// Solve R11 (krank x krank) z = b[0:krank]; remaining unknowns = 0 (basic solution)
   double z[];
   ArrayResize(z, n);
   for(int j = 0; j < n; j++)
      z[j] = 0.0;
   for(int j = krank-1; j >= 0; j--)
     {
      double s = b[bo+j];
      for(int i = j+1; i < krank; i++)
         s -= a[ao+j+i*lda]*z[i];
      z[j] = s / a[ao+j+j*lda];
     }
   for(int j = 0; j < n; j++)
      b[bo + pcols[j]] = z[j];
  }

//+------------------------------------------------------------------+
//| NNLS: Lawson-Hanson non-negative least squares (min |Ax-b|,      |
//| x>=0), via Householder QR + active-set updates. a is m x n,      |
//| overwritten; b is length m, overwritten; x is length n (output); |
//| w, zz are length-n / length-m scratch. All double arrays use the |
//| (array[], offset) convention since callers slice them out of a   |
//| shared scratch buffer; indices is always a dedicated int array   |
//| (never sliced from a shared buffer anywhere in this codebase).   |
//| info: 1=success, 2=bad dims, 3=hit maxiter.                      |
//+------------------------------------------------------------------+
void SLSQP_nnls(const int m, const int n,
                double &a[], const int ao, double &b[], const int bo,
                double &x[], const int xo, double &w[], const int wo,
                double &zz[], const int zzo,
                int &indices[], const int maxiter, double &rnorm, int &info)
  {
   int i=0, ii=0, ip=0, indz=0, iteration=0, iz=0, izmax=0;
   int j=0, jj=0, k=0;
   double tau=0.0, unorm=0.0, ztest=0.0, alpha=0.0, cc=0.0, ss=0.0, wmax=0.0, T=0.0;
   double pivot=1.0, pivot2=0.0, tmp=0.0, spacing=0.0;

   info = 1;
   if(m <= 0 || n <= 0)
     {
      info = 2;
      return;
     }

   for(i = 0; i < n; i++)
      indices[i] = i;
   for(i = 0; i < n; i++)
      x[xo+i] = 0.0;

   bool terminate = false;
   while(indz < MathMin(m,n) && !terminate)
     {
      for(i = indz; i < n; i++)
        {
         j = indices[i];
         int tmpint = m - indz;
         w[wo+j] = SLSQP_ddot(tmpint, a, ao + indz + j*m, 1, b, bo + indz, 1);
        }

      bool found_pivot = false;
      while(!found_pivot && !terminate)
        {
         wmax = 0.0;
         for(k = indz; k < n; k++)
           {
            j = indices[k];
            if(w[wo+j] > wmax)
              {
               wmax = w[wo+j];
               izmax = k;
              }
           }
         if(wmax <= 0.0)
           {
            terminate = true;
            break;
           }
         iz = izmax;
         j = indices[iz];

         pivot = a[ao + indz + j*m];
         int tmpint = m - indz;
         double tauArr[1];
         tauArr[0] = tau;
         double pivArr[1];
         pivArr[0] = pivot;
         SLSQP_dlarfgp(tmpint, pivArr, 0, a, ao + indz + 1 + j*m, 1, tauArr, 0);
         pivot = pivArr[0];
         tau = tauArr[0];

         unorm = (indz > 0 ? SLSQP_dnrm2(indz, a, ao + j*m, 1) : 0.0);
         spacing = (unorm > 0.0 ? unorm*SLSQP_EPS : 0.0); // approximates nextafter(unorm,2*unorm)-unorm

         if(MathAbs(pivot) > 100.0*spacing)
           {
            for(i = 0; i < m; i++)
               zz[zzo+i] = b[bo+i];
            tmpint = m - indz;
            pivot2 = a[ao + indz + j*m];
            a[ao + indz + j*m] = 1.0;
            double workv[1];
            SLSQP_dlarf(true, tmpint, 1, a, ao + indz + j*m, 1, tau, zz, zzo + indz, tmpint, workv, 0);
            ztest = zz[zzo+indz] / pivot;
            if(ztest > 0.0)
              {
               found_pivot = true;
               break;
              }
            else
              {
               a[ao + indz + j*m] = pivot2;
              }
           }
         w[wo+j] = 0.0;
        }
      if(terminate)
         break;

      for(i = 0; i < m; i++)
         b[bo+i] = zz[zzo+i];
      indices[iz] = indices[indz];
      indices[indz] = j;
      indz++;
      if(indz < n)
        {
         int tmpint = m - indz + 1;
         for(k = indz; k < n; k++)
           {
            jj = indices[k];
            double workv[1];
            SLSQP_dlarf(true, tmpint, 1, a, ao + indz - 1 + j*m, 1, tau, a, ao + indz - 1 + jj*m, tmpint, workv, 0);
           }
        }
      a[ao + indz - 1 + j*m] = pivot;
      if(indz < m)
        {
         for(i = indz; i < m; i++)
            a[ao + j*m + i] = 0.0;
        }
      w[wo+j] = 0.0;

      for(k = 0; k < indz; k++)
        {
         ip = indz - 1 - k;
         if(k != 0)
           {
            for(i = 0; i <= ip; i++)
               zz[zzo+i] = zz[zzo+i] - a[ao + i + jj*m]*zz[zzo+ip+1];
           }
         jj = indices[ip];
         zz[zzo+ip] = zz[zzo+ip] / a[ao + ip + jj*m];
        }

      while(true)
        {
         iteration++;
         if(iteration >= maxiter)
           {
            info = 3;
            terminate = true;
            break;
           }

         alpha = 2.0;
         for(ip = 0; ip < indz; ip++)
           {
            k = indices[ip];
            if(zz[zzo+ip] <= 0.0)
              {
               T = -x[xo+k] / (zz[zzo+ip] - x[xo+k]);
               if(alpha > T)
                 {
                  alpha = T;
                  jj = ip;
                 }
              }
           }
         if(alpha == 2.0)
            break;

         for(ip = 0; ip < indz; ip++)
           {
            k = indices[ip];
            x[xo+k] = x[xo+k] + alpha*(zz[zzo+ip]-x[xo+k]);
           }

         i = indices[jj];
         while(true)
           {
            x[xo+i] = 0.0;

            if(jj != indz-1)
              {
               jj++;
               for(j = jj; j < indz; j++)
                 {
                  ii = indices[j];
                  indices[j-1] = ii;
                  double csv, ssv, rv;
                  SLSQP_dlartgp(a[ao + j-1+ii*m], a[ao + j+ii*m], csv, ssv, rv);
                  cc = csv;
                  ss = ssv;
                  a[ao + j-1+ii*m] = rv;
                  a[ao + j+ii*m] = 0.0;
                  for(k = 0; k < n; k++)
                    {
                     if(k != ii)
                       {
                        tmp = a[ao + j-1+k*m];
                        a[ao + j-1+k*m] = cc*tmp + ss*a[ao + j+k*m];
                        a[ao + j+k*m]   = -ss*tmp + cc*a[ao + j+k*m];
                       }
                    }
                  tmp = b[bo+j-1];
                  b[bo+j-1] = cc*tmp + ss*b[bo+j];
                  b[bo+j]   = -ss*tmp + cc*b[bo+j];
                 }
              }
            indz--;
            indices[indz] = i;

            bool nobreak = false;
            for(jj = 0; jj < indz; jj++)
              {
               i = indices[jj];
               if(x[xo+i] <= 0.0)
                 {
                  break;
                 }
               if(jj == indz-1)
                  nobreak = true;
              }
            if(nobreak)
               break;
           }

         for(i = 0; i < m; i++)
            zz[zzo+i] = b[bo+i];
         for(k = 0; k < indz; k++)
           {
            ip = indz - 1 - k;
            if(k != 0)
              {
               for(i = 0; i <= ip; i++)
                  zz[zzo+i] = zz[zzo+i] - a[ao + i+jj*m]*zz[zzo+ip+1];
              }
            jj = indices[ip];
            zz[zzo+ip] = zz[zzo+ip] / a[ao + ip+jj*m];
           }
        }
      if(terminate)
         break;

      for(k = 0; k < indz; k++)
        {
         i = indices[k];
         x[xo+i] = zz[zzo+k];
        }
     }

   if(indz < m)
     {
      int tmpint = m - indz;
      rnorm = SLSQP_dnrm2(tmpint, b, bo + indz, 1);
     }
   else
     {
      for(i = 0; i < n; i++)
         w[wo+i] = 0.0;
      rnorm = 0.0;
     }
  }

//+------------------------------------------------------------------+
//| ldl_update: rank-1 update of a packed lower-triangular LDL'      |
//| factorization (BFGS Hessian approximation storage): a := LDL'    |
//| update of a +/- sigma*z*z^T. w is length-n scratch (only used    |
//| for sigma<0). a holds L off-diagonal, D on the diagonal.         |
//+------------------------------------------------------------------+
void SLSQP_ldl_update(const int n, double &a[], const int ao, double &z[], const int zo,
                      const double sigma, double &w[], const int wo)
  {
   int j, ij = 0;
   if(sigma == 0.0)
      return;
   double alpha, beta, delta, gamma, u, v, tp, t = 1.0/sigma;

   if(sigma <= 0.0)
     {
      for(int i = 0; i < n; i++)
         w[wo+i] = z[zo+i];
      for(int i = 0; i < n; i++)
        {
         v = w[wo+i];
         t = t + v*v/a[ao+ij];
         for(j = i+1; j < n; j++)
           {
            ij++;
            w[wo+j] = w[wo+j] - v*a[ao+ij];
           }
         ij++;
        }
      if(t >= 0.0)
         t = SLSQP_EPS/sigma;

      for(int i = 0; i < n; i++)
        {
         j = n - i - 1;
         ij -= i + 1;
         u = w[wo+j];
         w[wo+j] = t;
         t = t - u*u/a[ao+ij];
        }
     }

   for(int i = 0; i < n; i++)
     {
      v = z[zo+i];
      delta = v / a[ao+ij];
      tp = (sigma < 0.0 ? w[wo+i] : t + delta*v);
      alpha = tp / t;
      a[ao+ij] = alpha*a[ao+ij];
      if(i == n-1)
         return;
      beta = delta / tp;
      if(alpha <= 4.0)
        {
         for(j = i+1; j < n; j++)
           {
            ij++;
            z[zo+j] = z[zo+j] - v*a[ao+ij];
            a[ao+ij] = a[ao+ij] + beta*z[zo+j];
           }
        }
      else
        {
         gamma = t / tp;
         for(j = i+1; j < n; j++)
           {
            ij++;
            u = a[ao+ij];
            a[ao+ij] = gamma*u + beta*z[zo+j];
            z[zo+j] = z[zo+j] - v*u;
           }
        }
      ij++;
      t = tp;
     }
  }

//+------------------------------------------------------------------+
//| ldp: least-distance programming. min |x| s.t. G*x >= h            |
//| (via the NNLS dual, Lawson & Hanson).                             |
//| g is m x n, h is length m, x is length n (output).                |
//| buffer must have room for (m+2)*(n+1)+m+n doubles; indices needs  |
//| length >= n+1. mode: 1=solved, 2=bad dims, 3=nnls iter exceeded,  |
//| 4=inequalities incompatible/infeasible.                           |
//+------------------------------------------------------------------+
void SLSQP_ldp(const int m, const int n, double &g[], const int go,
               double &h[], const int ho, double &x[], const int xo,
               double &buffer[], const int bo, int &indices[],
               double &xnorm, int &mode)
  {
   if(n <= 0)
     {
      mode = 2;
      return;
     }
   for(int i = 0; i < n; i++)
      x[xo+i] = 0.0;
   if(m == 0)
     {
      mode = 1;
      return;
     }

   int a_o  = bo;
   int b_o  = bo + m*(n+1);
   int zz_o = bo + (m+1)*(n+1);
   int y_o  = bo + (m+2)*(n+1);
   int w_o  = bo + (m+2)*(n+1) + m;

   for(int j = 0; j < m; j++)
     {
      for(int i = 0; i < n; i++)
         buffer[a_o + i + j*(n+1)] = g[go + j + i*m];
      buffer[a_o + n + j*(n+1)] = h[ho+j];
     }
   for(int i = 0; i < n; i++)
      buffer[b_o+i] = 0.0;
   buffer[b_o+n] = 1.0;

   double rnorm = 0.0;
   SLSQP_nnls(n+1, m, buffer, a_o, buffer, b_o, buffer, y_o, buffer, w_o, buffer, zz_o,
              indices, 3*m, rnorm, mode);
   if(mode != 1)
      return;
   mode = 4;
   if(rnorm <= 0.0)
      return;

   double fac = 1.0 - SLSQP_ddot(m, h, ho, 1, buffer, y_o, 1);
   if(!((1.0+fac) - 1.0 > 0.0))
      return;
   mode = 1;
   fac = 1.0/fac;
   SLSQP_dgemv(true, m, n, fac, g, go, m, buffer, y_o, 1, 0.0, x, xo, 1);
   xnorm = SLSQP_dnrm2(n, x, xo, 1);

   for(int i = 0; i < m; i++)
      buffer[bo+i] = fac*buffer[y_o+i];
  }

//+-------------------------------------------------------------------+
//| lsi: least-squares with inequality constraints.                   |
//|   min |A*x - b|  s.t.  G*x >= h                                   |
//| A is ma x n (overwritten with its QR factors), b is length ma     |
//| (overwritten), g is mg x n (overwritten), h is length mg          |
//| (overwritten), x is length n (output).                            |
//| buffer needs room for the QR scratch (>= n) plus ldp's requirement|
//| ((mg+2)*(n+1)+2*mg); jw (indices) needs length >= n+1.            |
//| mode: 1 ok, 2 bad dims, 3 nnls iter exceeded, 4 incompatible,     |
//| 5 rank-deficient A.                                               |
//+-------------------------------------------------------------------+
void SLSQP_lsi(const int ma, const int mg, const int n,
               double &a[], const int ao, double &b[], const int bo,
               double &g[], const int go, double &h[], const int ho,
               double &x[], const int xo, double &buffer[], const int bufo,
               int &jw[], double &xnorm, int &mode)
  {
   int tmp_int = MathMin(ma,n);
   double work_dummy[1];

   SLSQP_dgeqr2(ma, n, a, ao, ma, buffer, bufo, buffer, bufo + tmp_int);
   SLSQP_dorm2r(true, ma, 1, tmp_int, a, ao, ma, buffer, bufo, b, bo, ma, buffer, bufo + tmp_int);

   mode = 5;
   xnorm = 0.0;
   for(int i = 0; i < tmp_int; i++)
      if(!(MathAbs(a[ao + i + i*ma]) >= SLSQP_EPS))
         return;

   SLSQP_dtrsm_RUNN(mg, n, 1.0, a, ao, ma, g, go, mg);
   SLSQP_dgemv(false, mg, n, -1.0, g, go, mg, b, bo, 1, 1.0, h, ho, 1);

   SLSQP_ldp(mg, n, g, go, h, ho, x, xo, buffer, bufo, jw, xnorm, mode);
   if(mode != 1)
      return;

   SLSQP_daxpy(n, 1.0, b, bo, 1, x, xo, 1);
   SLSQP_dtrsv(false, n, a, ao, ma, x, xo, 1);

   tmp_int = ma - n;
   int btail = (n+1 > ma ? ma : n+1) - 1;
   double tmp_dbl = SLSQP_dnrm2(tmp_int, b, bo + btail, 1);
   xnorm = SLSQP_Hypot(xnorm, tmp_dbl);
  }

//+------------------------------------------------------------------+
//| lsei: least-squares with equality AND inequality constraints.    |
//|   min |A*x - b|  s.t.  E*x = f,  G*x >= h                        |
//| a is ma x n, b is length ma, e is me x n, f is length me,        |
//| g is mg x n, h is length mg (all overwritten); x is length n.    |
//| buffer needs (mg+2)*(n-me+1) + 3*mg + 2*me + ma + (ma+mg)*(n-me);|
//| jw (indices) needs length >= n - me + 1.                         |
//| mode: 1 ok, 2 over-constrained, 4 incompatible, 5 rank-deficient |
//| A (in the mg==0 sub-case), 6 rank-deficient E, 7 rank-deficient A|
//| (mg==0 case, unsolvable).                                        |
//+------------------------------------------------------------------+
void SLSQP_lsei(const int ma, const int me, const int mg, const int n,
                double &a[], const int ao, double &b[], const int bo,
                double &e[], const int eo, double &f[], const int fo,
                double &g[], const int go, double &h[], const int ho,
                double &x[], const int xo, double &buffer[], const int bufo,
                int &jw[], double &xnorm, int &mode)
  {
   for(int i = 0; i < n; i++)
      x[xo+i] = 0.0;
   if(me > n)
     {
      mode = 2;
      return;
     }

   int nvars = n - me;
   int gmults_o = bufo;
   int emults_o = bufo + mg;
   int wb_o     = bufo + me + mg;
   int tau_o    = bufo + me + mg + ma;
   int a2_o     = bufo + mg + 2*me + ma;
   int g2_o     = bufo + mg + 2*me + ma + ma*nvars;
   int lsis_o   = bufo + mg + 2*me + ma + (ma+mg)*nvars;

   int lde = (me > 0 ? me : 1);
   int ldg = (mg > 0 ? mg : 1);
   SLSQP_dgerq2(me, n, e, eo, lde, buffer, tau_o, buffer, lsis_o);

   SLSQP_dormr2(false, true, ma, n, me, e, eo, lde, buffer, tau_o, a, ao, ma, buffer, lsis_o);
   SLSQP_dormr2(false, true, mg, n, me, e, eo, lde, buffer, tau_o, g, go, ldg, buffer, lsis_o);

   for(int i = 0; i < me; i++)
      if(!(MathAbs(e[eo + i + (nvars+i)*me]) >= SLSQP_EPS))
        {
         mode = 6;
         return;
        }

   for(int i = 0; i < me; i++)
      x[xo + nvars+i] = f[fo+i];
   SLSQP_dtrsv(false, me, e, eo + nvars*me, lde, x, xo + nvars, 1);

   mode = 1;
   for(int i = 0; i < mg; i++)
      buffer[gmults_o+i] = 0.0;

   bool skip_rest = (me == n);

   if(!skip_rest)
     {
      for(int i = 0; i < ma; i++)
         buffer[wb_o+i] = b[bo+i];
      SLSQP_dgemv(false, ma, me, -1.0, a, ao + ma*nvars, ma, x, xo + nvars, 1, 1.0, buffer, wb_o, 1);

      for(int j = 0; j < nvars; j++)
        {
         for(int i = 0; i < ma; i++)
            buffer[a2_o + i + j*ma] = a[ao + i + j*ma];
         for(int i = 0; i < mg; i++)
            buffer[g2_o + i + j*mg] = g[go + i + j*mg];
        }

      if(mg == 0)
        {
         int lwork = ma*nvars + 3*nvars + 1;
         int wborig_o = lsis_o + lwork;
         for(int i = 0; i < ma; i++)
            buffer[wborig_o+i] = buffer[wb_o+i];

         int krank = 0;
         double t = MathSqrt(SLSQP_EPS);
         SLSQP_dgelsy(ma, nvars, buffer, a2_o, ma, buffer, wb_o, t, krank);

         for(int i = 0; i < nvars; i++)
            x[xo+i] = buffer[wb_o+i];

         SLSQP_dgemv(false, ma, nvars, 1.0, a, ao, ma, x, xo, 1, -1.0, buffer, wborig_o, 1);
         xnorm = SLSQP_dnrm2(ma, buffer, wborig_o, 1);

         mode = 7;
         if(krank < nvars)
            return;
         mode = 1;
         skip_rest = true;
        }

      if(!skip_rest)
        {
         SLSQP_dgemv(false, mg, me, -1.0, g, go + mg*nvars, ldg, x, xo + nvars, 1, 1.0, h, ho, 1);

         SLSQP_lsi(ma, mg, nvars, buffer, a2_o, buffer, wb_o, buffer, g2_o, h, ho, x, xo,
                   buffer, lsis_o, jw, xnorm, mode);

         for(int i = 0; i < mg; i++)
            buffer[gmults_o+i] = buffer[lsis_o+i];

         if(me == 0)
            return;

         double t2 = SLSQP_dnrm2(me, x, xo + nvars, 1);
         xnorm = SLSQP_Hypot(xnorm, t2);
         if(mode != 1)
            return;
        }
     }

// ORIGINAL_BASIS: convert the solution and multipliers back to the original basis.
   SLSQP_dgemv(false, ma, n, 1.0, a, ao, ma, x, xo, 1, -1.0, b, bo, 1);
   SLSQP_dgemv(true, ma, me, 1.0, a, ao + nvars*ma, ma, b, bo, 1, 0.0, f, fo, 1);
   SLSQP_dgemv(true, mg, me, -1.0, g, go + nvars*mg, ldg, buffer, gmults_o, 1, 1.0, f, fo, 1);

   SLSQP_dormr2(true, true, n, 1, me, e, eo, lde, buffer, tau_o, x, xo, n, buffer, lsis_o);

   for(int i = 0; i < me; i++)
      buffer[emults_o+i] = f[fo+i];
   SLSQP_dtrsv(true, me, e, eo + (n-me)*me, lde, buffer, emults_o, 1);
  }

//+------------------------------------------------------------------+
//| lsq: forms and solves the SQP sub-problem's quadratic-program     |
//| direction-finding step as an equivalent LSEI (least squares with  |
//| equality/inequality constraints) problem, using the packed LDL'   |
//| BFGS factor Lf as the QP's Hessian square root.                   |
//|                                                                    |
//|   Lf     : packed LDL' factor of the BFGS Hessian approx (n*(n+1)/2)|
//|   gradx  : gradient of the objective at the current point (n)     |
//|   C      : constraint Jacobian, m x n (meq equality rows first)   |
//|   d      : constraint values (m)                                  |
//|   xl,xu  : bounds (n; extended to n+1 in-place if augment)         |
//|   x      : output search direction (n, or n+1 if augment)         |
//|   y      : output Lagrange multipliers (m + 2*n)                  |
//|   jw     : integer scratch (length >= n - meq + 2)                 |
//|   augment: if true, solves the "inconsistent linearization"        |
//|            relaxation with one extra slack variable weighted by   |
//|            aug_weight, per Kraft's SQP paper section 2.2.3.        |
//| buffer must satisfy the (generous) size formula documented in the |
//| SLSQPBody buffer-sizing helper further down this file.             |
//+------------------------------------------------------------------+
void SLSQP_lsq(int m, const int meq, int n, const bool augment, const double aug_weight,
               double &Lf[], const int Lfo, double &gradx[], const int gradxo,
               double &C[], const int Co, double &d[], const int d_o,
               double &xl[], const int xlo, double &xu[], const int xuo,
               double &x[], const int xo, double &y[], const int yo,
               double &buffer[], const int bufo, int &jw[], int &mode)
  {
   int orign = n;
   int mineq = m - meq;
   double xnorm = 0.0;
   int cursor = 0;
   int ld = n;

   if(augment)
     {
      ld = n + 1;
      x[xo+n] = 1.0;
      xl[xlo+n] = 0.0;
      xu[xuo+n] = 1.0;
     }

   for(int i = 0; i < (ld+2)*ld; i++)
      buffer[bufo+i] = 0.0;
   int wA_o = bufo;
   int wb_o = bufo + ld*(ld+1);

   for(int j = 0; j < n; j++)
     {
      double diag = MathSqrt(Lf[Lfo+cursor]);
      cursor++;
      buffer[wA_o + j + j*ld] = diag;
      for(int i = j+1; i < n; i++)
        {
         buffer[wA_o + j + i*ld] = Lf[Lfo+cursor]*diag;
         cursor++;
        }
     }

   for(int i = 0; i < n; i++)
      buffer[wb_o+i] = gradx[gradxo+i];
   SLSQP_dtpsv(n, Lf, Lfo, buffer, wb_o, 1);
   cursor = 0;
   for(int i = 0; i < n; i++)
     {
      buffer[wb_o+i] /= -MathSqrt(Lf[Lfo+cursor]);
      cursor += n - i;
     }

   if(augment)
      buffer[wA_o + ld*ld - 1] = aug_weight;
   if(augment)
      n++;

   int wE_o = bufo + n*(n+1) + n;
   int wf_o = bufo + n*(n+1) + n + n*meq;
   if(meq > 0)
     {
      for(int j = 0; j < n-1; j++)
         for(int i = 0; i < meq; i++)
            buffer[wE_o + i + j*meq] = C[Co + i + j*m];

      if(augment)
        {
         for(int i = 0; i < meq; i++)
            buffer[wE_o + i + (n-1)*meq] = -d[d_o+i];
        }
      else
        {
         for(int i = 0; i < meq; i++)
            buffer[wE_o + i + (n-1)*meq] = C[Co + i + (n-1)*m];
        }
      for(int i = 0; i < meq; i++)
         buffer[wf_o+i] = -d[d_o+i];
     }

   int wG_o = bufo + n*(n+1) + n + n*meq + meq;
   int wh_o = bufo + n*(n+1) + n + n*meq + meq + (mineq + 2*n)*ld;
   for(int i = 0; i < (mineq + 2*n)*(ld+1); i++)
      buffer[wG_o+i] = 0.0;

   int nancount = 0;
   int nrow = mineq;
   if(m > meq)
      for(int i = 0; i < mineq; i++)
         buffer[wh_o+i] = -d[d_o + meq + i];

   for(int i = 0; i < n; i++)
     {
      if(SLSQP_IsNaN(xl[xlo+i]))
         nancount++;
      else
        {
         buffer[wh_o+nrow] = xl[xlo+i];
         nrow++;
        }
     }
   for(int i = 0; i < n; i++)
     {
      if(SLSQP_IsNaN(xu[xuo+i]))
         nancount++;
      else
        {
         buffer[wh_o+nrow] = -xu[xuo+i];
         nrow++;
        }
     }

   int n_wG_rows = mineq + 2*n - nancount;

   if(m > meq)
     {
      for(int j = 0; j < orign; j++)
         for(int i = 0; i < mineq; i++)
            buffer[wG_o + i + j*n_wG_rows] = C[Co + meq + i + j*m];
     }

   if(augment)
     {
      for(int i = 0; i < mineq; i++)
         buffer[wG_o + i + orign*n_wG_rows] = MathMax(-d[d_o + meq + i], 0.0);
     }

   nrow = mineq;
   for(int i = 0; i < n; i++)
     {
      if(!SLSQP_IsNaN(xl[xlo+i]))
        {
         buffer[wG_o + nrow + i*n_wG_rows] = 1.0;
         nrow++;
        }
     }
   for(int i = 0; i < n; i++)
     {
      if(!SLSQP_IsNaN(xu[xuo+i]))
        {
         buffer[wG_o + nrow + i*n_wG_rows] = -1.0;
         nrow++;
        }
     }

   int lsei_scratch_o = wh_o + mineq + 2*n;

   SLSQP_lsei(ld, meq, n_wG_rows, n, buffer, wA_o, buffer, wb_o, buffer, wE_o, buffer, wf_o,
              buffer, wG_o, buffer, wh_o, x, xo, buffer, lsei_scratch_o, jw, xnorm, mode);

   if(mode == 1)
     {
      for(int i = 0; i < meq; i++)
         y[yo+i] = buffer[lsei_scratch_o + i + n_wG_rows];
      for(int i = 0; i < mineq; i++)
         y[yo + meq + i] = buffer[lsei_scratch_o + i];
      double nanv = SLSQP_NaN();
      for(int i = 0; i < 2*n; i++)
         y[yo + m + i] = nanv;
     }

   for(int i = 0; i < n; i++)
     {
      if((!SLSQP_IsNaN(xl[xlo+i])) && (x[xo+i] < xl[xlo+i]))
         x[xo+i] = xl[xlo+i];
      else
         if((!SLSQP_IsNaN(xu[xuo+i])) && (x[xo+i] > xu[xuo+i]))
            x[xo+i] = xu[xuo+i];
     }
  }

//+------------------------------------------------------------------+
//| SLSQP algorithm state, persisted by the caller across calls to    |
//| SLSQPBody (this is the reverse-communication "instruction         |
//| pointer" -- see CSLSQPSolver below for the ready-to-use driver).  |
//+------------------------------------------------------------------+
struct SSLSQPVars
  {
   double            acc, alpha, f0, gs, h1, h2, h3, h4, t, t0, tol;
   int               exact, inconsistent, reset, iter, itermax, line, m, meq, mode, n;
  };
//+------------------------------------------------------------------+
//|Enumeration of states marking sequence of operations for optimizer|
//+------------------------------------------------------------------+
enum ESLSQPState
  {
   ST_MODE0, ST_RESET_BFGS, ST_ITER_START, ST_LINE_SEARCH,
   ST_MODE1, ST_LABEL255, ST_MODEM1
  };

//+------------------------------------------------------------------+
//| SLSQPBody: the core reverse-communication SLSQP step function,    |
//| a direct (goto-free) port of SciPy's __slsqp_body. On return,     |
//| S.mode tells the caller what to do next:                          |
//|    S.mode ==  1 : evaluate funx and d (constraint values) ONLY    |
//|                   at 'sol', then call again with S.mode still 1   |
//|    S.mode == -1 : evaluate funx, gradx, C (constraint Jacobian)   |
//|                   and d at 'sol', then call again with mode -1    |
//|    S.mode ==  0 : converged; 'sol' holds the solution              |
//|    otherwise    : terminated abnormally (see docs for mode codes) |
//| Before the very first call, set S.mode=0 and evaluate funx,       |
//| gradx, C, d at the starting point (this is exactly what           |
//| CSLSQPSolver::Minimize does for you).                              |
//+------------------------------------------------------------------+
void SLSQPBody(SSLSQPVars &S, double &funx, double &gradx[], double &C[], double &d[],
               double &sol[], double &mult[], double &xl[], double &xu[],
               double &buffer[], int &indices[])
  {
   int lda = (S.m > 0 ? S.m : 1);
   int j;
   double alfmin = 0.1;
   int n = S.n;
   int m = S.m;
   int n1 = n + 1;
   int n2 = n1*n/2;

   int bfgs_o       = 0;
   int x0_o         = n2;
   int mu_o         = n2 + n;
   int s_o          = n2 + n + m;
   int u_o          = n2 + n + m + n1;
   int v_o          = n2 + n + m + n1 + n1;
   int lsqbuf_o     = n2 + n + m + n1 + n1 + n1;

   bool badlin = false;
   ESLSQPState state;

   if(S.mode == 0)
      state = ST_MODE0;
   else
      if(S.mode == -1)
         state = ST_MODEM1;
      else
         if(S.mode == 1)
            state = ST_MODE1;
         else
            return;

   while(true)
     {
      switch(state)
        {
         case ST_MODE0:
           {
            S.exact = 0;
            S.acc = MathAbs(S.acc);
            S.tol = 10*S.acc;
            S.iter = 0;
            S.reset = 0;
            for(int i = 0; i < n; i++)
               buffer[s_o+i] = 0.0;
            for(int i = 0; i < m; i++)
               buffer[mu_o+i] = 0.0;
            state = ST_RESET_BFGS;
            continue;
           }

         case ST_RESET_BFGS:
           {
            S.reset++;
            if(S.reset > 5)
              {
               state = ST_LABEL255;
               continue;
              }
            for(int i = 0; i < n2; i++)
               buffer[bfgs_o+i] = 0.0;
            j = 0;
            for(int i = 0; i < n; i++)
              {
               buffer[bfgs_o+j] = 1.0;
               j += n - i;
              }
            state = ST_ITER_START;
            continue;
           }

         case ST_ITER_START:
           {
            S.mode = 9;
            if(S.iter >= S.itermax)
               return;
            S.iter++;

            for(int i = 0; i < n; i++)
              {
               buffer[u_o+i] = -sol[i] + xl[i];
               buffer[v_o+i] = -sol[i] + xu[i];
              }

            S.h4 = 1.0;
            SLSQP_lsq(m, S.meq, n, false, 0.0, buffer, bfgs_o, gradx, 0, C, 0, d, 0,
                      buffer, u_o, buffer, v_o, buffer, s_o, mult, 0, buffer, lsqbuf_o,
                      indices, S.mode);

            badlin = false;

            if((S.mode == 6) && (n == S.meq))
               S.mode = 4;

            if(S.mode == 4)
              {
               badlin = true;
               for(int i = 0; i < n; i++)
                  buffer[s_o+i] = 0.0;
               S.h3 = 0.0;
               double rho = 100.0;
               S.inconsistent = 0;
               while(true)
                 {
                  SLSQP_lsq(m, S.meq, n, true, rho, buffer, bfgs_o, gradx, 0, C, 0, d, 0,
                            buffer, u_o, buffer, v_o, buffer, s_o, mult, 0, buffer, lsqbuf_o,
                            indices, S.mode);
                  S.h4 = 1.0 - buffer[s_o+n];
                  if(S.mode == 4)
                    {
                     rho *= 10.0;
                     S.inconsistent++;
                     if(S.inconsistent > 5)
                        return;
                     continue;
                    }
                  else
                     if(S.mode != 1)
                       {
                        return;
                       }
                  break;
                 }
              }
            else
               if(S.mode != 1)
                 {
                  return;
                 }

            for(int i = 0; i < n; i++)
               buffer[v_o+i] = gradx[i];
            SLSQP_dgemv(true, m, n, -1.0, C, 0, lda, mult, 0, 1, 1.0, buffer, v_o, 1);

            S.f0 = funx;
            for(int i = 0; i < n; i++)
               buffer[x0_o+i] = sol[i];
            S.gs = SLSQP_ddot(n, gradx, 0, 1, buffer, s_o, 1);
            S.h1 = MathAbs(S.gs);
            S.h2 = 0.0;
            for(int jj = 0; jj < m; jj++)
              {
               if(jj < S.meq)
                  S.h3 = d[jj];
               else
                  S.h3 = 0.0;
               S.h2 = S.h2 + MathMax(-d[jj], S.h3);
               S.h3 = MathAbs(mult[jj]);
               buffer[mu_o+jj] = MathMax(S.h3, (buffer[mu_o+jj] + S.h3)/2.0);
               S.h1 = S.h1 + S.h3*MathAbs(d[jj]);
              }

            S.mode = 0;
            if((S.h1 < S.acc) && (S.h2 < S.acc) && (!badlin) && (!SLSQP_IsNaN(funx)))
               return;
            S.h1 = 0.0;
            for(int jj = 0; jj < m; jj++)
              {
               if(jj < S.meq)
                  S.h3 = d[jj];
               else
                  S.h3 = 0.0;
               S.h1 += buffer[mu_o+jj]*MathMax(-d[jj], S.h3);
              }
            S.t0 = funx + S.h1;
            S.h3 = S.gs - S.h1*S.h4;
            S.mode = 8;
            if(S.h3 >= 0.0)
              {
               state = ST_RESET_BFGS;
               continue;
              }

            S.line = 0;
            S.alpha = 1.0;
            state = ST_LINE_SEARCH;
            continue;
           }

         case ST_LINE_SEARCH:
           {
            S.line++;
            S.h3 = S.alpha * S.h3;
            SLSQP_dscal(n, S.alpha, buffer, s_o, 1);
            for(int i = 0; i < n; i++)
               sol[i] = buffer[x0_o+i];
            SLSQP_daxpy(n, 1.0, buffer, s_o, 1, sol, 0, 1);

            S.mode = 1;
            return;
           }

         case ST_MODE1:
           {
            S.t = funx;
            for(int jj = 0; jj < m; jj++)
              {
               if(jj < S.meq)
                  S.h1 = d[jj];
               else
                  S.h1 = 0.0;
               S.t = S.t + buffer[mu_o+jj]*MathMax(-d[jj], S.h1);
              }
            S.h1 = S.t - S.t0;

            if(!((S.h1 <= (S.h3/10.0)) || (S.line > 10)))
              {
               S.alpha = MathMax(S.h3/(2.0*(S.h3-S.h1)), alfmin);
               state = ST_LINE_SEARCH;
               continue;
              }

            S.h3 = 0.0;
            for(int jj = 0; jj < m; jj++)
              {
               if(jj < S.meq)
                  S.h1 = d[jj];
               else
                  S.h1 = 0.0;
               S.h3 = S.h3 + MathMax(-d[jj], S.h1);
              }
            if(
               ((MathAbs(funx - S.f0) < S.acc) || (SLSQP_dnrm2(n, buffer, s_o, 1) < S.acc)) &&
               (S.h3 < S.acc) &&
               (!badlin) &&
               (!SLSQP_IsNaN(funx))
            )
              {
               S.mode = 0;
               return;
              }
            else
              {
               S.mode = -1;
              }
            return;
           }

         case ST_LABEL255:
           {
            S.h3 = 0.0;
            for(int jj = 0; jj < m; jj++)
              {
               if(jj < S.meq)
                  S.h1 = d[jj];
               else
                  S.h1 = 0.0;
               S.h3 = S.h3 + MathMax(-d[jj], S.h1);
              }
            if(((MathAbs(funx - S.f0) < S.tol) || (SLSQP_dnrm2(n, buffer, s_o, 1) < S.tol)) &&
               (S.h3 < S.tol) &&
               (!badlin) &&
               (!SLSQP_IsNaN(funx))
              )
              {
               S.mode = 0;
              }
            else
              {
               S.mode = 8;
              }
            return;
           }

         case ST_MODEM1:
           {
            for(int i = 0; i < n; i++)
               buffer[u_o+i] = gradx[i];
            SLSQP_dgemv(true, m, n, -1.0, C, 0, lda, mult, 0, 1, 1.0, buffer, u_o, 1);
            for(int i = 0; i < n; i++)
               buffer[u_o+i] = buffer[u_o+i] - buffer[v_o+i];

            for(int i = 0; i < n; i++)
               buffer[v_o+i] = buffer[s_o+i];
            SLSQP_dtpmv(true, n, buffer, bfgs_o, buffer, v_o, 1);

            j = 0;
            for(int i = 0; i < n; i++)
              {
               buffer[v_o+i] = buffer[bfgs_o+j]*buffer[v_o+i];
               j += n - i;
              }

            SLSQP_dtpmv(false, n, buffer, bfgs_o, buffer, v_o, 1);

            S.h1 = SLSQP_ddot(n, buffer, s_o, 1, buffer, u_o, 1);
            S.h2 = SLSQP_ddot(n, buffer, s_o, 1, buffer, v_o, 1);
            S.h3 = 0.2*S.h2;
            if(S.h1 < S.h3)
              {
               S.h4 = (S.h2 - S.h3) / (S.h2 - S.h1);
               S.h1 = S.h3;
               double tmp_dbl = 1.0 - S.h4;
               SLSQP_dscal(n, S.h4, buffer, u_o, 1);
               SLSQP_daxpy(n, tmp_dbl, buffer, v_o, 1, buffer, u_o, 1);
              }

            if((S.h1 == 0.0) || (S.h2 == 0.0))
              {
               state = ST_RESET_BFGS;
               continue;
              }

            SLSQP_ldl_update(n, buffer, bfgs_o, buffer, u_o, 1.0/S.h1, buffer, v_o);
            SLSQP_ldl_update(n, buffer, bfgs_o, buffer, v_o, -1.0/S.h2, buffer, u_o);

            state = ST_ITER_START;
            continue;
           }
        }
     }
  }

//+------------------------------------------------------------------+
//| Scratch-buffer sizing, matching the (additive, deliberately      |
//| generous) formula documented in SciPy's slsqp.c, using n+1 in    |
//| place of n throughout so it also covers lsq's internal           |
//| "inconsistent linearization" augmented sub-problem.              |
//+------------------------------------------------------------------+
int SLSQP_BufferSize(const int n0, const int m, const int meq)
  {
   int n = n0 + 1;
   int mineq = m - meq;
   long total = 0;
   total += n*(n+1)/2 + m + 4*n + 3;                                             // SLSQP
   total += (long)(n+1)*(n+2) + (long)(n+1)*meq + m + (long)(mineq+2*n+2)*(n+1) + 3*n + 3; // LSQ
   total += mineq + 2*n + 2 + 2*meq + (n+1) + (long)(mineq+3*n+3)*(n+1-meq);      // LSEI
   total += (long)(mineq+2*n+2+2)*(n+2) + mineq + 2*n + 2;                       // LDP
   total += mineq + 2*n + 2;                                                     // NNLS
   total *= 3; // safety margin (validated generously in the C reference tests)
   return((int)MathMax(total, 64));
  }

//+------------------------------------------------------------------+
//|Optimization results                                              |
//+------------------------------------------------------------------+
// Container structure consolidating final calculation snapshots from a completed optimization run
struct OptimizeResult
  {
   int               return_code;        // Solver output code status mapping from the slsqp_result enum set
   int               nfeval;             // Total number of objective function evaluation loops executed
   int               niter;              // Total mathematical convergence iteration steps consumed
   vector            solution;           // Final calculated coordinates vector within the parameter space
   vector            objective_result;   // Final minimized objective cost function scalar result packed inside a vector
   vector            objective_gradient; // Calculated local slope vector at the final solution coordinates

   // Default initialization constructor
                     OptimizeResult(void)
     {
      return_code = WRONG_VALUE; // Sentinel flag representing a non-evaluated or initialized result state
      nfeval = niter = 0;
      solution = objective_result = objective_gradient = vector::Zeros(0);
     }

   // Parameterized instantiation constructor for explicit data binding
                     OptimizeResult(int rc,int feval,int iter,vector &x, vector& f, vector& g)
     {
      return_code = rc;
      nfeval = feval;
      niter = iter;
      solution = x;
      objective_result = f;
      objective_gradient = g;
     }

   // Copy constructor safeguarding native vector deep replication paths
                     OptimizeResult(OptimizeResult& other)
     {
      return_code = other.return_code;
      nfeval = other.nfeval;
      niter = other.niter;
      solution = other.solution;
      objective_result = other.objective_result;
      objective_gradient = other.objective_gradient;
     }

   // Assignment operator overload ensuring thread and allocation safe copying
   void              operator=(OptimizeResult& other)
     {
      return_code = other.return_code;
      nfeval = other.nfeval;
      niter = other.niter;
      solution = other.solution;
      objective_result = other.objective_result;
      objective_gradient = other.objective_gradient;
     }
  };

//+------------------------------------------------------------------+
//| OOP interface for SLSQP minimizer                                |
//+------------------------------------------------------------------+
// Main class orchestration wrapper driving Sequential Least Squares Programming optimization
// routines using clean object-oriented control handles.
class CSlsqp : public CObject
  {
protected:
   slsqp_constraint  m_eq_constraints[];  // Tracked collection of zero-equality mathematical constraints
   slsqp_constraint  m_ineq_constraints[];// Tracked collection of bounds-inequality mathematical constraints
   CFunctor*         m_obj;               // Pointer to objective;
   int               m_n, m_m, m_meq;
   SSLSQPVars        m_S;
   //-- Interface to internal solver
   void              evaluate(double &x[], double &f, double &grad[],
                              double &c[], double &jac[], const bool need_derivatives)
     {
      vector vx = vector::Zeros(m_n); // Initialize MQL5 native vector container matching decision parameters
      Copy(vx,x,m_n);
      //---
      double empty_array[];
      //---
      if(need_derivatives)
        {
         ObjReturn or = m_obj.fun_and_grad(vx);
         f = or.f;
         //---
         Copy(grad,or.g,m_n);
        }
      else
        {
         vector fv = m_obj.objective_function(vx);
         f = fv[0];
        }
      //---
      int ii = 0;
      for(uint i = 0; i<m_eq_constraints.Size(); ++i)
         if(need_derivatives)
            slsqp_eval_constraint(c,jac,m_eq_constraints[i],m_n,x,ii);
         else
            slsqp_eval_constraint(c,empty_array,m_eq_constraints[i],m_n,x,ii);
      //---
      for(uint i = 0; i<m_ineq_constraints.Size(); ++i)
         if(need_derivatives)
            slsqp_eval_constraint(c,jac,m_ineq_constraints[i],m_n,x,ii);
         else
            slsqp_eval_constraint(c,empty_array,m_ineq_constraints[i],m_n,x,ii);

     }

public:
                     CSlsqp(void)
     {
      m_n = 0;
      m_m = 0;
      m_meq = 0;
      m_S.acc = 1.e-8;
      m_S.itermax = 100;
      m_S.n = m_n;
      m_S.m = m_m;
      m_S.meq = m_m;
     }
                    ~CSlsqp(void) {}




   // Sets precision/accuracy targets for optimization subproblems
   void              SetAcc(double acc_)
     {
      m_S.acc = acc_;
     }
   // Sets the upper limit on the allowed number of function evaluation cycles
   void              SetMaxEval(int maxeval)
     {
      m_S.itermax = fabs(maxeval);
     }

   // Configures equality constraints with explicit tolerance parameters
   bool              SetEqualityConstraints(CConstraints &constraints)
     {
      if(!m_eq_constraints.Size())
         if(ArrayResize(m_eq_constraints,1)!=1)
           {
            Print(__FUNCTION__,": Error ", GetLastError());
            return false;
           }

      m_eq_constraints[0].f_data = GetPointer(constraints);
      m_eq_constraints[0].m = constraints.numConstraints();
      return true;
     }

   // Configures inequality constraints with explicit tolerance parameters
   bool              SetInequalityConstraints(CConstraints &constraints)
     {
      if(!m_ineq_constraints.Size())
         if(ArrayResize(m_ineq_constraints,1)!=1)
           {
            Print(__FUNCTION__, ": Error ", GetLastError());
            return false;
           }

      m_ineq_constraints[0].f_data = GetPointer(constraints);
      m_ineq_constraints[0].m = constraints.numConstraints();
      return true;
     }

   //  Get string representation of optimizer's exit code
   string            GetExitMode(const int mode)
     {
      string message = "";
      switch(mode)
        {
         case -2:
            message = "Forced termination by user";
            break;
         case -1:
            message = "Gradient evaluation required (g & a)";
            break;
         case  0:
            message = "Optimization terminated successfully";
            break;
         case  1:
            message = "Function evaluation required (f & c)";
            break;
         case  2:
            message = "More equality constraints than independent variables";
            break;
         case  3:
            message = "More than 3*n iterations in LSQ subproblem";
            break;
         case  4:
            message = "Inequality constraints incompatible";
            break;
         case  5:
            message = "Singular matrix E in LSQ subproblem";
            break;
         case  6:
            message = "Singular matrix C in LSQ subproblem";
            break;
         case  7:
            message = "Rank-deficient equality constraint subproblem HFTI";
            break;
         case  8:
            message = "Positive directional derivative for linesearch";
            break;
         case  9:
            message = "Iteration limit reached";
            break;
         default:
            message = NULL;
            break;
        }
      return message;
     }

   // Main processing pipeline wrapper driving mathematical optimization mechanics over target parameters
   OptimizeResult    Minimize(CFunctor &fungrad, bool display_log_info = false)
     {
      vector x0 = fungrad.initial_params(); // Query starting guess parameters from user optimization configuration class
      int n = (int)x0.Size();
      m_n = n;
      vector x_data = x0;
      //---
      vector low = fungrad.lower_bounds();
      vector up = fungrad.upper_bounds();
      // Dimension checking to ensure shape match constraints for boundary conditions are satisfied
      if((low.Size()&&low.Size()!=ulong(n)) || (up.Size()&&up.Size()!=ulong(n)))
        {
         Print(__FUNCTION__,": All vector inputs should be the same size if not empty");
         return OptimizeResult();
        }

      // Verify bounds ordering logic
      for(int i  = 0; i<n; ++i)
         if(low[i]>up[i])
           {
            Print(__FUNCTION__ ": Invalid boundary constraints");
            return OptimizeResult();
           }

      m_obj = GetPointer(fungrad);
      m_meq = slsqp_count_constraints(ArraySize(m_eq_constraints),m_eq_constraints);
      m_m = m_meq + slsqp_count_constraints(ArraySize(m_ineq_constraints),m_ineq_constraints);

      m_S.mode = 0;
      m_S.n = m_n;
      m_S.m = m_m;
      m_S.meq = m_meq;

      int m = m_m;
      double funx = 0.0;

      double gradx[];
      ArrayResize(gradx, n);
      double C[];
      ArrayResize(C, MathMax(m*n,1));
      double d[];
      ArrayResize(d, MathMax(m,1));
      double mult[];
      ArrayResize(mult, m + 2*n + 2);

      ArrayInitialize(mult, 0.0);
      double xlw[];
      ArrayResize(xlw, n+1);
      double xuw[];
      ArrayResize(xuw, n+1);
      for(int i = 0; i < n; i++)
        {
         xlw[i] = low[i];
         xuw[i] = up[i];
        }

      int bufsize = SLSQP_BufferSize(m_n, m_m, m_meq);
      double buffer[];
      ArrayResize(buffer, bufsize);
      ArrayInitialize(buffer, 0.0);
      int indices[];
      ArrayResize(indices, m + 2*n + 2 + 16);

      double sol[];
      ArrayResize(sol, n+1);
      for(int i = 0; i < n; i++)
         sol[i] = x0[i];

      evaluate(sol, funx, gradx, d, C, true);

      int iter = 0;
      int prev_iter = -1;
      int max_driver_loops = 50*(m_S.itermax + 10);
      vector grad_v(n);
      vector fvector(1);
      //---
      if(display_log_info)
         PrintFormat("%5s %5s %16s %16s","NIT", "FC", "OBJFUN", "GNORM");
      //---
      while(!IsStopped())
        {
         iter++;
         if(iter > max_driver_loops)
           {
            if(display_log_info && m_S.iter != prev_iter)
              {
               vector gnorm;
               gnorm.Assign(gradx);
               PrintFormat("%5i %5i % 16.6E % 16.6E", m_S.iter, iter, funx, gnorm.Norm(VECTOR_NORM_P));
              }
            m_S.mode = 9;
            break;
           }
         SLSQPBody(m_S, funx, gradx, C, d, sol, mult, xlw, xuw, buffer, indices);
         if(m_S.mode == 1)
           {
            evaluate(sol, funx, gradx, d, C, false);
           }
         else
            if(m_S.mode == -1)
              {
               evaluate(sol, funx, gradx, d, C, true);
              }
            else
              {
               if(display_log_info && m_S.iter != prev_iter)
                 {
                  vector gnorm;
                  gnorm.Assign(gradx);
                  PrintFormat("%5i %5i % 16.6E % 16.6E", m_S.iter, iter, funx, gnorm.Norm(VECTOR_NORM_P));
                 }
               break;
              }

         if(display_log_info && m_S.iter != prev_iter)
           {
            vector gnorm;
            gnorm.Assign(gradx);
            PrintFormat("%5i %5i % 16.6E % 16.6E", m_S.iter, iter, funx, gnorm.Norm(VECTOR_NORM_P));
           }

         prev_iter = m_S.iter;
        }
      //---
      Copy(x0,sol,n);
      Copy(grad_v,gradx,n);
      fvector[0] = funx;
      if(IsStopped())
         m_S.mode = -2;
      //---
      if(display_log_info)
        {
         Print("Exit mode: ", GetExitMode(m_S.mode));
         Print("            Current function value:", funx);
         Print("            Iterations:", m_S.iter);
         Print("            Function evaluations:", m_obj.nfev());
         Print("            Gradient evaluations:", m_obj.ngev());
        }
      // Pack successful calculations tracking profiles into result layout container blocks
      return OptimizeResult(int(m_S.mode),m_S.iter,iter,x0,fvector,grad_v);
     }
  };
//+------------------------------------------------------------------+