ufranco/Article-23989-APARCH-Volatility-Process
Watch
1
0
Fork
You've already forked Article-23989-APARCH-Volatility-Process
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-23989-APARCH-Volati.../slsqp.mqh

2276 lines
76 KiB
MQL5
Raw Permalink Normal View History

First Commit
2026-08-10 23:47:33 +02:00
//+------------------------------------------------------------------+
//| 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);
}
};
//+------------------------------------------------------------------+
Reference in a new issue Copy permalink