Article-22610-Joint-Recurre.../Include/RQA/JRQAWindow.mqh
2026-05-30 13:14:34 +05:00

625 lines
24 KiB
MQL5

//+------------------------------------------------------------------+
//| JRQAWindow.mqh |
//| RQA Library for MQL5 |
//| Windowed (rolling) Joint-RQA for two time series |
//| |
//| OpenCL-accelerated: GPU computes both self-recurrence matrices |
//| and their AND in a single kernel pass. |
//| Falls back to CPU-only if OpenCL is unavailable. |
//+------------------------------------------------------------------+
#ifndef JRQAWINDOW_MQH
#define JRQAWINDOW_MQH
#include "JRQAMatrix.mqh"
#include "JRQAMetrics.mqh"
#include <OpenCL/OpenCL.mqh>
//+------------------------------------------------------------------+
//| OpenCL kernel — joint recurrence in a single pass |
//| For each (i,j,w): compute distX = ||x_i - x_j|| and |
//| distY = ||y_i - y_j||, output JR = (distX<=epsX)&&(distY<=epsY) |
//+------------------------------------------------------------------+
const string cl_jrqa_source =
"__kernel void jrqa_recurrence( \r\n"
" __global const float *seriesX, \r\n"
" __global const float *seriesY, \r\n"
" __global int *outR, \r\n"
" const int N, \r\n"
" const int embDim, \r\n"
" const int tau, \r\n"
" const int norm, \r\n"
" const float epsilonX, \r\n"
" const float epsilonY, \r\n"
" const int step, \r\n"
" const int baseWin) \r\n"
"{ \r\n"
" int i = get_global_id(0); \r\n"
" int j = get_global_id(1); \r\n"
" int w = get_global_id(2); \r\n"
" if(i >= N || j >= N) return; \r\n"
" int winStart = (baseWin + w) * step; \r\n"
" float distX = 0.0f; \r\n"
" float distY = 0.0f; \r\n"
" if(embDim == 1) { \r\n"
" float diffX = seriesX[winStart+i] - seriesX[winStart+j]; \r\n"
" float diffY = seriesY[winStart+i] - seriesY[winStart+j]; \r\n"
" if(norm == 1) { distX = diffX*diffX; distY = diffY*diffY;}\r\n"
" else { distX = fabs(diffX); distY = fabs(diffY);}\r\n"
" } else { \r\n"
" for(int d = 0; d < embDim; d++) { \r\n"
" float diffX = seriesX[winStart + i + d*tau] \r\n"
" - seriesX[winStart + j + d*tau]; \r\n"
" float diffY = seriesY[winStart + i + d*tau] \r\n"
" - seriesY[winStart + j + d*tau]; \r\n"
" if(norm == 0) { \r\n"
" distX = fmax(distX, fabs(diffX)); \r\n"
" distY = fmax(distY, fabs(diffY)); \r\n"
" } else if(norm == 1) { \r\n"
" distX += diffX*diffX; \r\n"
" distY += diffY*diffY; \r\n"
" } else { \r\n"
" distX += fabs(diffX); \r\n"
" distY += fabs(diffY); \r\n"
" } \r\n"
" } \r\n"
" } \r\n"
" float threshX = (norm==1) ? epsilonX*epsilonX : epsilonX; \r\n"
" float threshY = (norm==1) ? epsilonY*epsilonY : epsilonY; \r\n"
" outR[(w*N + i)*N + j] = \r\n"
" ((distX <= threshX) && (distY <= threshY)) ? 1 : 0; \r\n"
"} \r\n";
//+------------------------------------------------------------------+
//| Per-window result struct |
//+------------------------------------------------------------------+
struct SJRQAWindowResult
{
int barIndex;
SJRQAResult metrics;
};
//+------------------------------------------------------------------+
//| CJRQAWindow — rolling Joint-RQA with GPU/CPU |
//+------------------------------------------------------------------+
class CJRQAWindow
{
private:
int m_windowSize;
int m_step;
double m_epsilonX;
double m_epsilonY;
int m_embDim;
int m_delay;
ENUM_RQA_NORM m_norm;
int m_minDiagLine;
int m_minVertLine;
bool RunGPU(const double &seriesX[], const double &seriesY[],
int seriesLen, int numWindows,
SJRQAWindowResult &results[]);
void ComputeFusedCPU(const double &sX[], int offX,
const double &sY[], int offY,
SJRQAResult &result);
void ScanMetrics(const int &R[], int N, int baseIdx,
SJRQAResult &result);
public:
CJRQAWindow();
void SetWindow(int windowSize, int step = 1) { m_windowSize = windowSize; m_step = step; }
void SetEpsilon(double epsX, double epsY) { m_epsilonX = epsX; m_epsilonY = epsY; }
void SetEpsilon(double eps) { m_epsilonX = eps; m_epsilonY = eps; }
void SetEmbedding(int dim, int delay) { m_embDim = dim; m_delay = delay; }
void SetNorm(ENUM_RQA_NORM norm) { m_norm = norm; }
void SetMinLines(int diagMin, int vertMin) { m_minDiagLine = diagMin; m_minVertLine = vertMin; }
bool Run(const double &seriesX[], int lenX,
const double &seriesY[], int lenY,
SJRQAWindowResult &results[]);
static void ExtractJRR (const SJRQAWindowResult &r[], double &out[]);
static void ExtractJDET (const SJRQAWindowResult &r[], double &out[]);
static void ExtractJLAM (const SJRQAWindowResult &r[], double &out[]);
static void ExtractJTT (const SJRQAWindowResult &r[], double &out[]);
static void ExtractJENTR(const SJRQAWindowResult &r[], double &out[]);
static void ExtractJLmax(const SJRQAWindowResult &r[], double &out[]);
};
//+------------------------------------------------------------------+
//| Default constructor |
//+------------------------------------------------------------------+
CJRQAWindow::CJRQAWindow()
: m_windowSize(50), m_step(1), m_epsilonX(0.1), m_epsilonY(0.1),
m_embDim(1), m_delay(1), m_norm(RQA_NORM_EUCLIDEAN),
m_minDiagLine(2), m_minVertLine(2)
{
}
//+------------------------------------------------------------------+
//| Scan a flat int R[N*N] sub-array for JRQA metrics |
//+------------------------------------------------------------------+
void CJRQAWindow::ScanMetrics(const int &R[], int N, int baseIdx,
SJRQAResult &result)
{
result.Reset();
if(N < 2) return;
long NsqMinusN = (long)N * N - N;
long recCount = 0;
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
if(i != j && R[baseIdx + i * N + j] != 0)
recCount++;
result.JRR = (NsqMinusN > 0) ? (double)recCount / NsqMinusN : 0.0;
// --- Diagonal lines (excluding main diagonal) ---
int diagHist[];
ArrayResize(diagHist, N + 1);
ArrayInitialize(diagHist, 0);
for(int k = -(N - 1); k <= (N - 1); k++)
{
if(k == 0) continue;
int iS = (k < 0) ? -k : 0;
int iE = (k < 0) ? N - 1 : N - 1 - k;
int len = 0;
for(int i = iS; i <= iE; i++)
{
if(R[baseIdx + i * N + (i + k)] != 0)
len++;
else
{
if(len >= m_minDiagLine) diagHist[len]++;
len = 0;
}
}
if(len >= m_minDiagLine) diagHist[len]++;
}
long diagPoints = 0, totalDiagLines = 0;
int lmax = 0;
for(int l = m_minDiagLine; l <= N; l++)
if(diagHist[l] > 0)
{ diagPoints += (long)l * diagHist[l]; totalDiagLines += diagHist[l]; lmax = l; }
result.JLmax = (double)lmax;
result.JDIV = (lmax > 0) ? 1.0 / lmax : 0.0;
result.JDET = (recCount > 0) ? (double)diagPoints / recCount : 0.0;
result.JL = (totalDiagLines > 0) ? (double)diagPoints / totalDiagLines : 0.0;
if(totalDiagLines > 0)
{
double entr = 0.0;
for(int l = m_minDiagLine; l <= N; l++)
if(diagHist[l] > 0)
{ double p = (double)diagHist[l] / totalDiagLines; entr -= p * MathLog(p); }
result.JENTR = entr;
}
// --- Vertical lines ---
int vertHist[];
ArrayResize(vertHist, N + 1);
ArrayInitialize(vertHist, 0);
for(int j = 0; j < N; j++)
{
int len = 0;
for(int i = 0; i < N; i++)
{
if(R[baseIdx + i * N + j] != 0)
len++;
else
{ if(len >= m_minVertLine) vertHist[len]++; len = 0; }
}
if(len >= m_minVertLine) vertHist[len]++;
}
long vertPoints = 0, totalVertLines = 0;
int vmax = 0;
for(int l = m_minVertLine; l <= N; l++)
if(vertHist[l] > 0)
{ vertPoints += (long)l * vertHist[l]; totalVertLines += vertHist[l]; vmax = l; }
result.JVmax = (double)vmax;
result.JLAM = (recCount > 0) ? (double)vertPoints / recCount : 0.0;
result.JTT = (totalVertLines > 0) ? (double)vertPoints / totalVertLines : 0.0;
result.JRATIO = (result.JRR > 1e-12) ? result.JDET / result.JRR : 0.0;
result.JCOMPLEXITY = result.JRR * result.JDET;
// --- TREND ---
if(N >= 4)
{
int numDiag = N - 1;
double sumX = 0, sumYd = 0, sumXY = 0, sumX2 = 0;
for(int d = 1; d < N; d++)
{
int cnt = 0, tot = N - d;
for(int i = 0; i < tot; i++)
if(R[baseIdx + i * N + (i + d)] != 0) cnt++;
double density = (tot > 0) ? (double)cnt / tot : 0.0;
double x = (double)(d - numDiag / 2);
sumX += x;
sumYd += density;
sumXY += x * density;
sumX2 += x * x;
}
double denom = (double)numDiag * sumX2 - sumX * sumX;
if(MathAbs(denom) > 1e-12)
result.JTREND = ((double)numDiag * sumXY - sumX * sumYd) / denom;
}
}
//+------------------------------------------------------------------+
//| GPU path — using COpenCL wrapper |
//+------------------------------------------------------------------+
bool CJRQAWindow::RunGPU(const double &seriesX[], const double &seriesY[],
int seriesLen, int numWindows,
SJRQAWindowResult &results[])
{
int N = m_windowSize - (m_embDim - 1) * m_delay;
if(N < 2) return false;
float fSeriesX[], fSeriesY[];
ArrayResize(fSeriesX, seriesLen);
ArrayResize(fSeriesY, seriesLen);
for(int i = 0; i < seriesLen; i++)
{
fSeriesX[i] = (float)seriesX[i];
fSeriesY[i] = (float)seriesY[i];
}
COpenCL ocl;
if(!ocl.Initialize(cl_jrqa_source, true))
{
Print("JRQA: OpenCL init failed");
return false;
}
if(!ocl.SetKernelsCount(1) || !ocl.KernelCreate(0, "jrqa_recurrence"))
{
Print("JRQA: kernel create failed");
ocl.Shutdown();
return false;
}
if(!ocl.SetBuffersCount(3))
{
ocl.Shutdown();
return false;
}
if(!ocl.BufferFromArray(0, fSeriesX, 0, seriesLen, CL_MEM_READ_ONLY) ||
!ocl.BufferFromArray(1, fSeriesY, 0, seriesLen, CL_MEM_READ_ONLY))
{
Print("JRQA: buffer upload failed");
ocl.Shutdown();
return false;
}
ocl.SetArgumentBuffer(0, 0, 0); // seriesX
ocl.SetArgumentBuffer(0, 1, 1); // seriesY
ocl.SetArgument(0, 3, N);
ocl.SetArgument(0, 4, m_embDim);
ocl.SetArgument(0, 5, m_delay);
ocl.SetArgument(0, 6, (int)m_norm);
ocl.SetArgument(0, 7, (float)m_epsilonX);
ocl.SetArgument(0, 8, (float)m_epsilonY);
ocl.SetArgument(0, 9, m_step);
long cellsPerWin = (long)N * N;
int maxBatch = (int)MathMin((long)numWindows,
64L * 1024 * 1024 / (cellsPerWin * (long)sizeof(int)));
if(maxBatch < 1) maxBatch = 1;
ArrayResize(results, numWindows);
bool ok = true;
int processed = 0;
while(processed < numWindows && ok)
{
int batchSize = MathMin(maxBatch, numWindows - processed);
int totalCells = batchSize * (int)cellsPerWin;
ocl.BufferFree(2);
if(!ocl.BufferCreate(2, totalCells * sizeof(int), CL_MEM_WRITE_ONLY))
{ ok = false; break; }
ocl.SetArgumentBuffer(0, 2, 2); // outR
ocl.SetArgument(0, 10, processed);
uint gOff[3] = {0, 0, 0};
uint gWork[3] = {(uint)N, (uint)N, (uint)batchSize};
if(!ocl.Execute(0, 3, gOff, gWork))
{ ok = false; break; }
int rBuf[];
ArrayResize(rBuf, totalCells);
if(!ocl.BufferRead(2, rBuf, 0, 0, totalCells))
{ ok = false; break; }
for(int b = 0; b < batchSize; b++)
{
int winIdx = processed + b;
results[winIdx].barIndex = winIdx * m_step;
ScanMetrics(rBuf, N, b * (int)cellsPerWin, results[winIdx].metrics);
}
processed += batchSize;
}
ocl.Shutdown();
return ok;
}
//+------------------------------------------------------------------+
//| CPU fallback — fused single-window joint recurrence compute |
//+------------------------------------------------------------------+
void CJRQAWindow::ComputeFusedCPU(const double &sX[], int offX,
const double &sY[], int offY,
SJRQAResult &result)
{
result.Reset();
int N = m_windowSize - (m_embDim - 1) * m_delay;
if(N <= 1) return;
long NsqMinusN = (long)N * N - N;
double epsSqX = m_epsilonX * m_epsilonX;
double epsSqY = m_epsilonY * m_epsilonY;
int diagHist[], vertHist[];
ArrayResize(diagHist, N + 1);
ArrayResize(vertHist, N + 1);
ArrayInitialize(diagHist, 0);
ArrayInitialize(vertHist, 0);
long recCount = 0;
if(m_embDim == 1 && m_norm == RQA_NORM_EUCLIDEAN)
{
// --- Fast path: embDim=1, Euclidean ---
// Diagonal scan
for(int k = -(N - 1); k <= (N - 1); k++)
{
if(k == 0) continue;
int iS = (k < 0) ? -k : 0;
int iE = (k < 0) ? N - 1 : N - 1 - k;
int len = 0;
for(int i = iS; i <= iE; i++)
{
double diffX = sX[offX + i] - sX[offX + i + k];
double diffY = sY[offY + i] - sY[offY + i + k];
if(diffX * diffX <= epsSqX && diffY * diffY <= epsSqY)
{ len++; recCount++; }
else
{ if(len >= m_minDiagLine) diagHist[len]++; len = 0; }
}
if(len >= m_minDiagLine) diagHist[len]++;
}
// Vertical scan
for(int j = 0; j < N; j++)
{
int len = 0;
for(int i = 0; i < N; i++)
{
if(i == j) { if(len >= m_minVertLine) vertHist[len]++; len = 0; continue; }
double diffX = sX[offX + i] - sX[offX + j];
double diffY = sY[offY + i] - sY[offY + j];
if(diffX * diffX <= epsSqX && diffY * diffY <= epsSqY)
len++;
else
{ if(len >= m_minVertLine) vertHist[len]++; len = 0; }
}
if(len >= m_minVertLine) vertHist[len]++;
}
}
else
{
// --- General path ---
bool useEucSq = (m_norm == RQA_NORM_EUCLIDEAN);
double epsCompX = useEucSq ? epsSqX : m_epsilonX;
double epsCompY = useEucSq ? epsSqY : m_epsilonY;
// Diagonal scan
for(int k = -(N - 1); k <= (N - 1); k++)
{
if(k == 0) continue;
int iS = (k < 0) ? -k : 0;
int iE = (k < 0) ? N - 1 : N - 1 - k;
int len = 0;
for(int i = iS; i <= iE; i++)
{
double dstX = 0.0, dstY = 0.0;
for(int d = 0; d < m_embDim; d++)
{
double diffX = sX[offX + i + d * m_delay]
- sX[offX + i + k + d * m_delay];
double diffY = sY[offY + i + d * m_delay]
- sY[offY + i + k + d * m_delay];
if(useEucSq)
{ dstX += diffX * diffX; dstY += diffY * diffY; }
else if(m_norm == RQA_NORM_MAX)
{ dstX = MathMax(dstX, MathAbs(diffX)); dstY = MathMax(dstY, MathAbs(diffY)); }
else
{ dstX += MathAbs(diffX); dstY += MathAbs(diffY); }
}
if(dstX <= epsCompX && dstY <= epsCompY)
{ len++; recCount++; }
else
{ if(len >= m_minDiagLine) diagHist[len]++; len = 0; }
}
if(len >= m_minDiagLine) diagHist[len]++;
}
// Vertical scan
for(int j = 0; j < N; j++)
{
int len = 0;
for(int i = 0; i < N; i++)
{
if(i == j) { if(len >= m_minVertLine) vertHist[len]++; len = 0; continue; }
double dstX = 0.0, dstY = 0.0;
for(int d = 0; d < m_embDim; d++)
{
double diffX = sX[offX + i + d * m_delay]
- sX[offX + j + d * m_delay];
double diffY = sY[offY + i + d * m_delay]
- sY[offY + j + d * m_delay];
if(useEucSq)
{ dstX += diffX * diffX; dstY += diffY * diffY; }
else if(m_norm == RQA_NORM_MAX)
{ dstX = MathMax(dstX, MathAbs(diffX)); dstY = MathMax(dstY, MathAbs(diffY)); }
else
{ dstX += MathAbs(diffX); dstY += MathAbs(diffY); }
}
if(dstX <= epsCompX && dstY <= epsCompY)
len++;
else
{ if(len >= m_minVertLine) vertHist[len]++; len = 0; }
}
if(len >= m_minVertLine) vertHist[len]++;
}
}
// --- Assemble metrics ---
result.JRR = (NsqMinusN > 0) ? (double)recCount / NsqMinusN : 0.0;
long diagPoints = 0, totalDiagLines = 0;
int lmax = 0;
for(int l = m_minDiagLine; l <= N; l++)
if(diagHist[l] > 0)
{ diagPoints += (long)l * diagHist[l]; totalDiagLines += diagHist[l]; lmax = l; }
result.JLmax = (double)lmax;
result.JDIV = (lmax > 0) ? 1.0 / lmax : 0.0;
result.JDET = (recCount > 0) ? (double)diagPoints / recCount : 0.0;
result.JL = (totalDiagLines > 0) ? (double)diagPoints / totalDiagLines : 0.0;
if(totalDiagLines > 0)
{
double entr = 0.0;
for(int l = m_minDiagLine; l <= N; l++)
if(diagHist[l] > 0)
{ double p = (double)diagHist[l] / totalDiagLines; entr -= p * MathLog(p); }
result.JENTR = entr;
}
long vertPoints = 0, totalVertLines = 0;
int vmax = 0;
for(int l = m_minVertLine; l <= N; l++)
if(vertHist[l] > 0)
{ vertPoints += (long)l * vertHist[l]; totalVertLines += vertHist[l]; vmax = l; }
result.JVmax = (double)vmax;
result.JLAM = (recCount > 0) ? (double)vertPoints / recCount : 0.0;
result.JTT = (totalVertLines > 0) ? (double)vertPoints / totalVertLines : 0.0;
result.JRATIO = (result.JRR > 1e-12) ? result.JDET / result.JRR : 0.0;
result.JCOMPLEXITY = result.JRR * result.JDET;
// --- TREND ---
if(N >= 4)
{
int numDiag = N - 1;
double sumX = 0, sumYd = 0, sumXY = 0, sumX2 = 0;
for(int dd = 1; dd < N; dd++)
{
int cnt = 0, tot = N - dd;
for(int i = 0; i < tot; i++)
{
double diffX, diffY;
if(m_embDim == 1)
{
diffX = sX[offX + i] - sX[offX + i + dd];
diffY = sY[offY + i] - sY[offY + i + dd];
bool recX = (m_norm == RQA_NORM_EUCLIDEAN)
? (diffX * diffX <= epsSqX)
: (MathAbs(diffX) <= m_epsilonX);
bool recY = (m_norm == RQA_NORM_EUCLIDEAN)
? (diffY * diffY <= epsSqY)
: (MathAbs(diffY) <= m_epsilonY);
if(recX && recY) cnt++;
}
else
{
bool useEucSq = (m_norm == RQA_NORM_EUCLIDEAN);
double dstX = 0.0, dstY = 0.0;
for(int d = 0; d < m_embDim; d++)
{
double dfX = sX[offX + i + d * m_delay] - sX[offX + i + dd + d * m_delay];
double dfY = sY[offY + i + d * m_delay] - sY[offY + i + dd + d * m_delay];
if(useEucSq) { dstX += dfX * dfX; dstY += dfY * dfY; }
else if(m_norm == RQA_NORM_MAX)
{ dstX = MathMax(dstX, MathAbs(dfX)); dstY = MathMax(dstY, MathAbs(dfY)); }
else { dstX += MathAbs(dfX); dstY += MathAbs(dfY); }
}
double ecX = useEucSq ? epsSqX : m_epsilonX;
double ecY = useEucSq ? epsSqY : m_epsilonY;
if(dstX <= ecX && dstY <= ecY) cnt++;
}
}
double density = (tot > 0) ? (double)cnt / tot : 0.0;
double x = (double)(dd - numDiag / 2);
sumX += x;
sumYd += density;
sumXY += x * density;
sumX2 += x * x;
}
double denom = (double)numDiag * sumX2 - sumX * sumX;
if(MathAbs(denom) > 1e-12)
result.JTREND = ((double)numDiag * sumXY - sumX * sumYd) / denom;
}
}
//+------------------------------------------------------------------+
//| Main entry — tries GPU, falls back to CPU |
//+------------------------------------------------------------------+
bool CJRQAWindow::Run(const double &seriesX[], int lenX,
const double &seriesY[], int lenY,
SJRQAWindowResult &results[])
{
int minLen = MathMin(lenX, lenY);
if(minLen < m_windowSize)
{
Print("CJRQAWindow::Run — series shorter than window");
return false;
}
int numWindows = (minLen - m_windowSize) / m_step + 1;
if(RunGPU(seriesX, seriesY, minLen, numWindows, results))
return true;
Print("JRQA: OpenCL unavailable, using CPU fallback");
ArrayResize(results, numWindows);
for(int idx = 0; idx < numWindows; idx++)
{
int start = idx * m_step;
results[idx].barIndex = start;
ComputeFusedCPU(seriesX, start, seriesY, start, results[idx].metrics);
}
return true;
}
//+------------------------------------------------------------------+
//| Extract single-metric arrays from windowed JRQA results |
//+------------------------------------------------------------------+
void CJRQAWindow::ExtractJRR(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JRR; }
void CJRQAWindow::ExtractJDET(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JDET; }
void CJRQAWindow::ExtractJLAM(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JLAM; }
void CJRQAWindow::ExtractJTT(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JTT; }
void CJRQAWindow::ExtractJENTR(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JENTR; }
void CJRQAWindow::ExtractJLmax(const SJRQAWindowResult &r[], double &out[])
{ int n=ArraySize(r); ArrayResize(out,n); for(int i=0;i<n;i++) out[i]=r[i].metrics.JLmax; }
#endif // JRQAWINDOW_MQH