//+------------------------------------------------------------------+ //| 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 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