Gamestudio Links
Zorro Links
Newest Posts
Purchase A8 full licence version
by ukgamer. 04/29/26 18:09
Z9 getting Error 058
by k_ivan. 04/25/26 19:13
ZorroGPT
by TipmyPip. 04/25/26 16:09
Stooq now requires an API key
by jcl. 04/13/26 09:42
Strange "Alien" Skull created with >Knubber<
by NeoDumont. 04/10/26 18:58
AUM Magazine
Latest Screens
Dorifto samurai
Shadow 2
Rocker`s Revenge
Stug 3 Stormartillery
Who's Online Now
0 registered members (), 5,489 guests, and 68 spiders.
Key: Admin, Global Mod, Mod
Newest Members
ukgamer, valino, juergenwue, VladMak, Geir
19210 Registered Users
Previous Thread
Next Thread
Print Thread
Rate Thread
Page 23 of 23 1 2 21 22 23
The Cartographer of Compact Currents [Re: TipmyPip] #489355
04/08/26 17:07
04/08/26 17:07
Joined: Sep 2017
Posts: 293
TipmyPip Offline OP
Member
TipmyPip  Offline OP
Member

Joined: Sep 2017
Posts: 293
A multi asset Zorro strategy engine that watches a basket of currency pairs, turns each pair into a compact feature profile, builds a live map of how pairs are related through behavior and shared currency exposure, measures which assets sit in the most coherent and central parts of that map, and then refines those rankings through several learning layers that estimate market regime, novelty, stability, density, and uncertainty, using an optional accelerated compute path when available and a full processor fallback when not, so that instead of blindly picking the strongest raw signals it produces a cautious, diversified, regime aware shortlist of the most structurally attractive pairs at each update.

Code
// TGr06A_CompactDominant_v14.cpp - Zorro64 Strategy DLL
// Strategy A v13: Compactness-Dominant with MX06 OOP + OpenCL + Learning Controller
//
// Notes:
// - Keeps full CPU fallback.
// - OpenCL is optional: if OpenCL.dll missing / no device / kernel build fails -> CPU path.
// - OpenCL accelerates the heavy correlation matrix step by offloading pairwise correlations.
// - Correlation is computed in float on GPU; results are stored back into fvar corrMatrix.

#define _CRT_SECURE_NO_WARNINGS
#include <zorro.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <windows.h>
#include <stddef.h>

#define INF 1e30
#define EPS 1e-12
#define N_ASSETS 28
#define FEAT_N 9
#define FEAT_WINDOW 200
#define UPDATE_EVERY 5
#define TOP_K 5

#define ALPHA 0.1
#define BETA 0.2
#define GAMMA 3.0
#define LAMBDA_META 0.7

#define USE_ML 1
#define USE_UNSUP 1
#define USE_RL 1
#define USE_PCA 1
#define USE_GMM 1
#define USE_HMM 1
#define HMM_K 3
#define HMM_DIM 8
#define HMM_VAR_FLOOR 1e-4
#define HMM_SMOOTH 0.02
#define HMM_ENTROPY_TH 0.85
#define HMM_SWITCH_TH 0.35
#define HMM_MIN_RISK 0.25
#define HMM_COOLDOWN_UPDATES 2
#define HMM_ONLINE_UPDATE 1
#define USE_KMEANS 1
#define KMEANS_K 3
#define KMEANS_DIM 8
#define KMEANS_ETA 0.03
#define KMEANS_DIST_EMA 0.08
#define KMEANS_STABILITY_MIN 0.35
#define KMEANS_ONLINE_UPDATE 1
#define USE_SPECTRAL 1
#define SPECTRAL_K 4
#define USE_HCLUST 1
#define HCLUST_COARSE_K 4
#define HCLUST_FINE_K 8
#define USE_COMMUNITY 1
#define COMM_W_MIN 0.15
#define COMM_TOPM 6
#define COMM_ITERS 4
#define COMM_Q_EMA 0.20
#define COMM_Q_LOW 0.20
#define COMM_Q_HIGH 0.45
#define USE_AE 1
#define AE_INPUT_DIM 8
#define AE_LATENT_DIM 4
#define AE_NORM_ALPHA 0.02
#define AE_ERR_EMA 0.10
#define AE_Z_LOW 1.0
#define AE_Z_HIGH 2.0
#define USE_SOM 1
#define SOM_W 10
#define SOM_H 10
#define SOM_DIM 12
#define SOM_ALPHA_MAX 0.30
#define SOM_ALPHA_MIN 0.05
#define SOM_SIGMA_MAX 5.0
#define SOM_SIGMA_MIN 1.0
#define SOM_CONF_MIN 0.15
#define SOM_ONLINE_UPDATE 1
#define USE_DENSITY 1
#define DENSITY_DIM 10
#define DENSITY_HISTORY 640
#define DENSITY_EPS 2.35
#define DENSITY_MIN_SAMPLES 10
#define DENSITY_N_REGIMES 4
#define DENSITY_NOISE_PERSIST 2
#define DENSITY_CLEAR_PERSIST 2
#define DENSITY_Z_EMA 0.03
#define GMM_K 3
#define GMM_DIM 8
#define GMM_ALPHA 0.02
#define GMM_VAR_FLOOR 1e-4
#define GMM_ENTROPY_COEFF 0.45
#define GMM_MIN_RISK 0.25
#define GMM_ONLINE_UPDATE 1
#define STRATEGY_PROFILE 0
#define PCA_DIM 6
#define PCA_COMP 3
#define PCA_WINDOW 128
#define PCA_REBUILD_EVERY 4

#ifdef TIGHT_MEM
typedef float fvar;
#else
typedef double fvar;
#endif

static const char* ASSET_NAMES[] = {
  "EURUSD","GBPUSD","USDCHF","USDJPY","AUDUSD","AUDCAD","AUDCHF","AUDJPY","AUDNZD",
  "CADJPY","CADCHF","EURAUD","EURCAD","EURCHF","EURGBP","EURJPY","EURNZD","GBPAUD",
  "GBPCAD","GBPCHF","GBPJPY","GBPNZD","NZDCAD","NZDCHF","NZDJPY","NZDUSD","USDCAD"
};
static const char* CURRENCIES[] = {"EUR","GBP","USD","CHF","JPY","AUD","CAD","NZD"};
#define N_CURRENCIES 8

// ---------------------------- Exposure Table ----------------------------

struct ExposureTable {
  int exposure[N_ASSETS][N_CURRENCIES];
  double exposureDist[N_ASSETS][N_ASSETS];

  void init() {
    for(int i=0;i<N_ASSETS;i++){
      for(int c=0;c<N_CURRENCIES;c++){
        exposure[i][c] = 0;
      }
    }
    for(int i=0;i<N_ASSETS;i++){
      for(int j=0;j<N_ASSETS;j++){
        exposureDist[i][j] = 0.0;
      }
    }
  }

  inline double getDist(int i,int j) const { return exposureDist[i][j]; }
};

// ---------------------------- Slab Allocator ----------------------------

template<typename T>
class SlabAllocator {
public:
  T* data;
  int capacity;

  SlabAllocator() : data(NULL), capacity(0) {}
  ~SlabAllocator() { shutdown(); }

  void init(int size) {
    shutdown();
    capacity = size;
    data = (T*)malloc((size_t)capacity * sizeof(T));
    if(data) memset(data, 0, (size_t)capacity * sizeof(T));
  }

  void shutdown() {
    if(data) free(data);
    data = NULL;
    capacity = 0;
  }

  T& operator[](int i) { return data[i]; }
  const T& operator[](int i) const { return data[i]; }
};

// ---------------------------- Feature Buffer (SoA ring) ----------------------------

struct FeatureBufferSoA {
  SlabAllocator<fvar> buffer;
  int windowSize;
  int currentIndex;

  void init(int assets, int window) {
    windowSize = window;
    currentIndex = 0;
    buffer.init(FEAT_N * assets * window);
  }

  void shutdown() { buffer.shutdown(); }

  inline int offset(int feat,int asset,int t) const {
    return (feat * N_ASSETS + asset) * windowSize + t;
  }

  void push(int feat,int asset,fvar value) {
    buffer[offset(feat, asset, currentIndex)] = value;
    currentIndex = (currentIndex + 1) % windowSize;
  }

  // t=0 => most recent
  fvar get(int feat,int asset,int t) const {
    int idx = (currentIndex - 1 - t + windowSize) % windowSize;
    return buffer[offset(feat, asset, idx)];
  }
};

// ---------------------------- Minimal OpenCL (dynamic) ----------------------------

typedef struct _cl_platform_id*   cl_platform_id;
typedef struct _cl_device_id*     cl_device_id;
typedef struct _cl_context*       cl_context;
typedef struct _cl_command_queue* cl_command_queue;
typedef struct _cl_program*       cl_program;
typedef struct _cl_kernel*        cl_kernel;
typedef struct _cl_mem*           cl_mem;
typedef unsigned int              cl_uint;
typedef int                       cl_int;
typedef unsigned long long        cl_ulong;
typedef size_t                    cl_bool;

#define CL_SUCCESS 0
#define CL_DEVICE_TYPE_CPU (1ULL << 1)
#define CL_DEVICE_TYPE_GPU (1ULL << 2)
#define CL_MEM_READ_ONLY   (1ULL << 2)
#define CL_MEM_WRITE_ONLY  (1ULL << 1)
#define CL_MEM_READ_WRITE  (1ULL << 0)
#define CL_TRUE  1
#define CL_FALSE 0
#define CL_PROGRAM_BUILD_LOG 0x1183

class OpenCLBackend {
public:
  HMODULE hOpenCL;
  int ready;

  cl_platform_id platform;
  cl_device_id device;
  cl_context context;
  cl_command_queue queue;
  cl_program program;
  cl_kernel kCorr;

  cl_mem bufFeat;
  cl_mem bufCorr;

  int featBytes;
  int corrBytes;

  cl_int (*clGetPlatformIDs)(cl_uint, cl_platform_id*, cl_uint*);
  cl_int (*clGetDeviceIDs)(cl_platform_id, cl_ulong, cl_uint, cl_device_id*, cl_uint*);
  cl_context (*clCreateContext)(void*, cl_uint, const cl_device_id*, void*, void*, cl_int*);
  cl_command_queue (*clCreateCommandQueue)(cl_context, cl_device_id, cl_ulong, cl_int*);
  cl_program (*clCreateProgramWithSource)(cl_context, cl_uint, const char**, const size_t*, cl_int*);
  cl_int (*clBuildProgram)(cl_program, cl_uint, const cl_device_id*, const char*, void*, void*);
  cl_int (*clGetProgramBuildInfo)(cl_program, cl_device_id, cl_uint, size_t, void*, size_t*);
  cl_kernel (*clCreateKernel)(cl_program, const char*, cl_int*);
  cl_int (*clSetKernelArg)(cl_kernel, cl_uint, size_t, const void*);
  cl_mem (*clCreateBuffer)(cl_context, cl_ulong, size_t, void*, cl_int*);
  cl_int (*clEnqueueWriteBuffer)(cl_command_queue, cl_mem, cl_bool, size_t, size_t, const void*, cl_uint, const void*, void*);
  cl_int (*clEnqueueReadBuffer)(cl_command_queue, cl_mem, cl_bool, size_t, size_t, void*, cl_uint, const void*, void*);
  cl_int (*clEnqueueNDRangeKernel)(cl_command_queue, cl_kernel, cl_uint, const size_t*, const size_t*, const size_t*, cl_uint, const void*, void*);
  cl_int (*clFinish)(cl_command_queue);
  cl_int (*clReleaseMemObject)(cl_mem);
  cl_int (*clReleaseKernel)(cl_kernel);
  cl_int (*clReleaseProgram)(cl_program);
  cl_int (*clReleaseCommandQueue)(cl_command_queue);
  cl_int (*clReleaseContext)(cl_context);

  OpenCLBackend()
  : hOpenCL(NULL), ready(0),
    platform(NULL), device(NULL), context(NULL), queue(NULL), program(NULL), kCorr(NULL),
    bufFeat(NULL), bufCorr(NULL),
    featBytes(0), corrBytes(0),
    clGetPlatformIDs(NULL), clGetDeviceIDs(NULL), clCreateContext(NULL), clCreateCommandQueue(NULL),
    clCreateProgramWithSource(NULL), clBuildProgram(NULL), clGetProgramBuildInfo(NULL),
    clCreateKernel(NULL), clSetKernelArg(NULL),
    clCreateBuffer(NULL), clEnqueueWriteBuffer(NULL), clEnqueueReadBuffer(NULL),
    clEnqueueNDRangeKernel(NULL), clFinish(NULL),
    clReleaseMemObject(NULL), clReleaseKernel(NULL), clReleaseProgram(NULL),
    clReleaseCommandQueue(NULL), clReleaseContext(NULL)
  {}

  int loadSymbol(void** fp, const char* name) {
    *fp = (void*)GetProcAddress(hOpenCL, name);
    return (*fp != NULL);
  }

  const char* kernelSource() {
    return
      "__kernel void corr_pairwise(\n"
      "  __global const float* feat,\n"
      "  __global float* outCorr,\n"
      "  const int nAssets,\n"
      "  const int nFeat,\n"
      "  const int windowSize,\n"
      "  const float eps\n"
      "){\n"
      "  int a = (int)get_global_id(0);\n"
      "  int b = (int)get_global_id(1);\n"
      "  if(a >= nAssets || b >= nAssets) return;\n"
      "  if(a >= b) return;\n"
      "  float acc = 0.0f;\n"
      "  for(int f=0; f<nFeat; f++){\n"
      "    int baseA = (f*nAssets + a) * windowSize;\n"
      "    int baseB = (f*nAssets + b) * windowSize;\n"
      "    float mx = 0.0f;\n"
      "    float my = 0.0f;\n"
      "    for(int t=0; t<windowSize; t++){\n"
      "      mx += feat[baseA + t];\n"
      "      my += feat[baseB + t];\n"
      "    }\n"
      "    mx /= (float)windowSize;\n"
      "    my /= (float)windowSize;\n"
      "    float sxx = 0.0f;\n"
      "    float syy = 0.0f;\n"
      "    float sxy = 0.0f;\n"
      "    for(int t=0; t<windowSize; t++){\n"
      "      float dx = feat[baseA + t] - mx;\n"
      "      float dy = feat[baseB + t] - my;\n"
      "      sxx += dx*dx;\n"
      "      syy += dy*dy;\n"
      "      sxy += dx*dy;\n"
      "    }\n"
      "    float den = sqrt(sxx*syy + eps);\n"
      "    float corr = (den > eps) ? (sxy/den) : 0.0f;\n"
      "    acc += corr;\n"
      "  }\n"
      "  outCorr[a*nAssets + b] = acc / (float)nFeat;\n"
      "}\n";
  }

  void printBuildLog() {
    if(!clGetProgramBuildInfo || !program || !device) return;
    size_t logSize = 0;
    clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &logSize);
    if(logSize == 0) return;
    char* log = (char*)malloc(logSize + 1);
    if(!log) return;
    memset(log, 0, logSize + 1);
    clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, logSize, log, NULL);
    printf("OpenCL build log:\n%s\n", log);
    free(log);
  }

  void init() {
    ready = 0;

    hOpenCL = LoadLibraryA("OpenCL.dll");
    if(!hOpenCL) {
      printf("OpenCL: CPU (OpenCL.dll missing)\n");
      return;
    }

    if(!loadSymbol((void**)&clGetPlatformIDs,       "clGetPlatformIDs")) return;
    if(!loadSymbol((void**)&clGetDeviceIDs,         "clGetDeviceIDs")) return;
    if(!loadSymbol((void**)&clCreateContext,        "clCreateContext")) return;
    if(!loadSymbol((void**)&clCreateCommandQueue,   "clCreateCommandQueue")) return;
    if(!loadSymbol((void**)&clCreateProgramWithSource,"clCreateProgramWithSource")) return;
    if(!loadSymbol((void**)&clBuildProgram,         "clBuildProgram")) return;
    if(!loadSymbol((void**)&clGetProgramBuildInfo,  "clGetProgramBuildInfo")) return;
    if(!loadSymbol((void**)&clCreateKernel,         "clCreateKernel")) return;
    if(!loadSymbol((void**)&clSetKernelArg,         "clSetKernelArg")) return;
    if(!loadSymbol((void**)&clCreateBuffer,         "clCreateBuffer")) return;
    if(!loadSymbol((void**)&clEnqueueWriteBuffer,   "clEnqueueWriteBuffer")) return;
    if(!loadSymbol((void**)&clEnqueueReadBuffer,    "clEnqueueReadBuffer")) return;
    if(!loadSymbol((void**)&clEnqueueNDRangeKernel, "clEnqueueNDRangeKernel")) return;
    if(!loadSymbol((void**)&clFinish,               "clFinish")) return;
    if(!loadSymbol((void**)&clReleaseMemObject,     "clReleaseMemObject")) return;
    if(!loadSymbol((void**)&clReleaseKernel,        "clReleaseKernel")) return;
    if(!loadSymbol((void**)&clReleaseProgram,       "clReleaseProgram")) return;
    if(!loadSymbol((void**)&clReleaseCommandQueue,  "clReleaseCommandQueue")) return;
    if(!loadSymbol((void**)&clReleaseContext,       "clReleaseContext")) return;

    cl_uint nPlat = 0;
    if(clGetPlatformIDs(0, NULL, &nPlat) != CL_SUCCESS || nPlat == 0) {
      printf("OpenCL: CPU (no platform)\n");
      return;
    }
    clGetPlatformIDs(1, &platform, NULL);

    cl_uint nDev = 0;
    cl_int ok = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, &nDev);
    if(ok != CL_SUCCESS || nDev == 0) {
      ok = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &device, &nDev);
      if(ok != CL_SUCCESS || nDev == 0) {
        printf("OpenCL: CPU (no device)\n");
        return;
      }
    }

    cl_int err = 0;
    context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
    if(err != CL_SUCCESS || !context) {
      printf("OpenCL: CPU (context fail)\n");
      return;
    }

    queue = clCreateCommandQueue(context, device, 0, &err);
    if(err != CL_SUCCESS || !queue) {
      printf("OpenCL: CPU (queue fail)\n");
      return;
    }

    const char* src = kernelSource();
    program = clCreateProgramWithSource(context, 1, &src, NULL, &err);
    if(err != CL_SUCCESS || !program) {
      printf("OpenCL: CPU (program fail)\n");
      return;
    }

    err = clBuildProgram(program, 1, &device, "", NULL, NULL);
    if(err != CL_SUCCESS) {
      printf("OpenCL: CPU (build fail)\n");
      printBuildLog();
      return;
    }

    kCorr = clCreateKernel(program, "corr_pairwise", &err);
    if(err != CL_SUCCESS || !kCorr) {
      printf("OpenCL: CPU (kernel fail)\n");
      printBuildLog();
      return;
    }

    featBytes = FEAT_N * N_ASSETS * FEAT_WINDOW * (int)sizeof(float);
    corrBytes = N_ASSETS * N_ASSETS * (int)sizeof(float);

    bufFeat = clCreateBuffer(context, CL_MEM_READ_ONLY, (size_t)featBytes, NULL, &err);
    if(err != CL_SUCCESS || !bufFeat) {
      printf("OpenCL: CPU (bufFeat fail)\n");
      return;
    }

    bufCorr = clCreateBuffer(context, CL_MEM_WRITE_ONLY, (size_t)corrBytes, NULL, &err);
    if(err != CL_SUCCESS || !bufCorr) {
      printf("OpenCL: CPU (bufCorr fail)\n");
      return;
    }

    ready = 1;
    printf("OpenCL: READY (kernel+buffers)\n");
  }

  void shutdown() {
    if(bufCorr) { clReleaseMemObject(bufCorr); bufCorr = NULL; }
    if(bufFeat) { clReleaseMemObject(bufFeat); bufFeat = NULL; }
    if(kCorr) { clReleaseKernel(kCorr); kCorr = NULL; }
    if(program) { clReleaseProgram(program); program = NULL; }
    if(queue) { clReleaseCommandQueue(queue); queue = NULL; }
    if(context) { clReleaseContext(context); context = NULL; }
    if(hOpenCL) { FreeLibrary(hOpenCL); hOpenCL = NULL; }
    ready = 0;
  }

  int computeCorrelationMatrixCL(const float* featLinear, float* outCorr, int nAssets, int nFeat, int windowSize) {
    if(!ready) return 0;
    if(!featLinear || !outCorr) return 0;

    cl_int err = clEnqueueWriteBuffer(queue, bufFeat, CL_TRUE, 0, (size_t)featBytes, featLinear, 0, NULL, NULL);
    if(err != CL_SUCCESS) return 0;

    float eps = 1e-12f;
    err = CL_SUCCESS;
    err |= clSetKernelArg(kCorr, 0, sizeof(cl_mem), &bufFeat);
    err |= clSetKernelArg(kCorr, 1, sizeof(cl_mem), &bufCorr);
    err |= clSetKernelArg(kCorr, 2, sizeof(int), &nAssets);
    err |= clSetKernelArg(kCorr, 3, sizeof(int), &nFeat);
    err |= clSetKernelArg(kCorr, 4, sizeof(int), &windowSize);
    err |= clSetKernelArg(kCorr, 5, sizeof(float), &eps);
    if(err != CL_SUCCESS) return 0;

    size_t global[2];
    global[0] = (size_t)nAssets;
    global[1] = (size_t)nAssets;

    err = clEnqueueNDRangeKernel(queue, kCorr, 2, NULL, global, NULL, 0, NULL, NULL);
    if(err != CL_SUCCESS) return 0;

    err = clFinish(queue);
    if(err != CL_SUCCESS) return 0;

    err = clEnqueueReadBuffer(queue, bufCorr, CL_TRUE, 0, (size_t)corrBytes, outCorr, 0, NULL, NULL);
    if(err != CL_SUCCESS) return 0;

    return 1;
  }
};

// ---------------------------- Learning Layer ----------------------------

struct LearningSnapshot {
  double meanScore;
  double meanCompactness;
  double meanVol;
  double meanAbsCorr;
  double stdAbsCorr;
  double momentumMean;
  int regime;
  double regimeConfidence;
};

class UnsupervisedModel {
public:
  double centroids[3][3];
  int counts[3];
  int initialized;

  UnsupervisedModel() : initialized(0) {
    memset(centroids, 0, sizeof(centroids));
    memset(counts, 0, sizeof(counts));
  }

  void init() {
    initialized = 0;
    memset(centroids, 0, sizeof(centroids));
    memset(counts, 0, sizeof(counts));
  }

  void update(const LearningSnapshot& s, int* regimeOut, double* confOut) {
    double x[3];
    x[0] = s.meanScore;
    x[1] = s.meanCompactness;
    x[2] = s.meanVol;

    if(!initialized) {
      for(int k=0; k<3; k++) {
        centroids[k][0] = x[0] + 0.01 * (k - 1);
        centroids[k][1] = x[1] + 0.01 * (1 - k);
        centroids[k][2] = x[2] + 0.005 * (k - 1);
        counts[k] = 1;
      }
      initialized = 1;
    }

    int best = 0;
    double bestDist = INF;
    double secondDist = INF;
    for(int k=0; k<3; k++) {
      double d0 = x[0] - centroids[k][0];
      double d1 = x[1] - centroids[k][1];
      double d2 = x[2] - centroids[k][2];
      double dist = d0*d0 + d1*d1 + d2*d2;
      if(dist < bestDist) {
        secondDist = bestDist;
        bestDist = dist;
        best = k;
      } else if(dist < secondDist) {
        secondDist = dist;
      }
    }

    counts[best]++;
    double lr = 1.0 / (double)counts[best];
    centroids[best][0] += lr * (x[0] - centroids[best][0]);
    centroids[best][1] += lr * (x[1] - centroids[best][1]);
    centroids[best][2] += lr * (x[2] - centroids[best][2]);

    *regimeOut = best;
    *confOut = 1.0 / (1.0 + sqrt(fabs(secondDist - bestDist) + EPS));
  }
};

class RLAgent {
public:
  double q[4];
  int n[4];
  double epsilon;
  int lastAction;
  double lastMeanScore;

  RLAgent() : epsilon(0.10), lastAction(0), lastMeanScore(0) {
    for(int i=0;i<4;i++){ q[i]=0; n[i]=0; }
  }

  void init() {
    epsilon = 0.10;
    lastAction = 0;
    lastMeanScore = 0;
    for(int i=0;i<4;i++){ q[i]=0; n[i]=0; }
  }

  int chooseAction(int updateCount) {
    int exploratory = ((updateCount % 10) == 0);
    if(exploratory) return updateCount % 4;
    int best = 0;
    for(int i=1;i<4;i++) if(q[i] > q[best]) best = i;
    return best;
  }

  void updateReward(double newMeanScore) {
    double reward = newMeanScore - lastMeanScore;
    n[lastAction]++;
    q[lastAction] += (reward - q[lastAction]) / (double)n[lastAction];
    lastMeanScore = newMeanScore;
  }
};

class PCAModel {
public:
  double hist[PCA_WINDOW][PCA_DIM];
  double mean[PCA_DIM];
  double stdev[PCA_DIM];
  double latent[PCA_COMP];
  double explainedVar[PCA_COMP];
  int writeIdx;
  int count;
  int rebuildEvery;
  int updates;
  double dom;
  double rot;
  double prevExplained0;

  PCAModel() : writeIdx(0), count(0), rebuildEvery(PCA_REBUILD_EVERY), updates(0), dom(0), rot(0), prevExplained0(0) {
    memset(hist, 0, sizeof(hist));
    memset(mean, 0, sizeof(mean));
    memset(stdev, 0, sizeof(stdev));
    memset(latent, 0, sizeof(latent));
    memset(explainedVar, 0, sizeof(explainedVar));
  }

  void init() {
    writeIdx = 0;
    count = 0;
    updates = 0;
    dom = 0;
    rot = 0;
    prevExplained0 = 0;
    memset(hist, 0, sizeof(hist));
    memset(mean, 0, sizeof(mean));
    memset(stdev, 0, sizeof(stdev));
    memset(latent, 0, sizeof(latent));
    memset(explainedVar, 0, sizeof(explainedVar));
  }

  void pushSnapshot(const double x[PCA_DIM]) {
    for(int d=0; d<PCA_DIM; d++) hist[writeIdx][d] = x[d];
    writeIdx = (writeIdx + 1) % PCA_WINDOW;
    if(count < PCA_WINDOW) count++;
  }

  void rebuildStats() {
    if(count <= 0) return;
    for(int d=0; d<PCA_DIM; d++) {
      double m = 0;
      for(int i=0; i<count; i++) m += hist[i][d];
      m /= (double)count;
      mean[d] = m;

      double v = 0;
      for(int i=0; i<count; i++) {
        double dd = hist[i][d] - m;
        v += dd * dd;
      }
      v /= (double)count;
      stdev[d] = sqrt(v + EPS);
    }
  }

  void update(const LearningSnapshot& snap, int regime, double conf) {
    double x[PCA_DIM];
    x[0] = snap.meanScore;
    x[1] = snap.meanCompactness;
    x[2] = snap.meanVol;
    x[3] = (double)regime / 2.0;
    x[4] = conf;
    x[5] = snap.meanScore - snap.meanCompactness;

    pushSnapshot(x);
    updates++;
    if((updates % rebuildEvery) == 0 || count < 4) rebuildStats();

    double z[PCA_DIM];
    for(int d=0; d<PCA_DIM; d++) z[d] = (x[d] - mean[d]) / (stdev[d] + EPS);

    latent[0] = 0.60*z[0] + 0.30*z[1] + 0.10*z[2];
    latent[1] = 0.25*z[0] - 0.45*z[1] + 0.20*z[2] + 0.10*z[4];
    latent[2] = 0.20*z[2] + 0.50*z[3] - 0.30*z[5];

    double a0 = fabs(latent[0]);
    double a1 = fabs(latent[1]);
    double a2 = fabs(latent[2]);
    double sumA = a0 + a1 + a2 + EPS;

    explainedVar[0] = a0 / sumA;
    explainedVar[1] = a1 / sumA;
    explainedVar[2] = a2 / sumA;

    dom = explainedVar[0];
    rot = fabs(explainedVar[0] - prevExplained0);
    prevExplained0 = explainedVar[0];
  }
};

class GMMRegimeModel {
public:
  double pi[GMM_K];
  double mu[GMM_K][GMM_DIM];
  double var[GMM_K][GMM_DIM];
  double p[GMM_K];
  double entropy;
  double conf;
  int bestRegime;
  int initialized;

  GMMRegimeModel() : entropy(0), conf(0), bestRegime(0), initialized(0) {
    memset(pi, 0, sizeof(pi));
    memset(mu, 0, sizeof(mu));
    memset(var, 0, sizeof(var));
    memset(p, 0, sizeof(p));
  }

  void init() {
    initialized = 0;
    entropy = 0;
    conf = 0;
    bestRegime = 0;
    for(int k=0;k<GMM_K;k++) {
      pi[k] = 1.0 / (double)GMM_K;
      for(int d=0; d<GMM_DIM; d++) {
        mu[k][d] = 0.02 * (k - 1);
        var[k][d] = 1.0;
      }
      p[k] = 1.0 / (double)GMM_K;
    }
    initialized = 1;
  }

  static double gaussianDiag(const double* x, const double* m, const double* v) {
    double logp = 0;
    for(int d=0; d<GMM_DIM; d++) {
      double vv = v[d];
      if(vv < GMM_VAR_FLOOR) vv = GMM_VAR_FLOOR;
      double z = x[d] - m[d];
      logp += -0.5 * (z*z / vv + log(vv + EPS));
    }
    if(logp < -80.0) logp = -80.0;
    return exp(logp);
  }

  void infer(const double x[GMM_DIM]) {
    if(!initialized) init();
    double sum = 0;
    for(int k=0;k<GMM_K;k++) {
      double g = gaussianDiag(x, mu[k], var[k]);
      p[k] = pi[k] * g;
      sum += p[k];
    }
    if(sum < EPS) {
      for(int k=0;k<GMM_K;k++) p[k] = 1.0 / (double)GMM_K;
    } else {
      for(int k=0;k<GMM_K;k++) p[k] /= sum;
    }

    bestRegime = 0;
    conf = p[0];
    for(int k=1;k<GMM_K;k++) {
      if(p[k] > conf) {
        conf = p[k];
        bestRegime = k;
      }
    }

    entropy = 0;
    for(int k=0;k<GMM_K;k++) entropy -= p[k] * log(p[k] + EPS);

#if GMM_ONLINE_UPDATE
    // lightweight incremental update (EM-like with forgetting)
    for(int k=0;k<GMM_K;k++) {
      double w = GMM_ALPHA * p[k];
      pi[k] = (1.0 - GMM_ALPHA) * pi[k] + w;
      for(int d=0; d<GMM_DIM; d++) {
        double diff = x[d] - mu[k][d];
        mu[k][d] += w * diff;
        var[k][d] = (1.0 - w) * var[k][d] + w * diff * diff;
        if(var[k][d] < GMM_VAR_FLOOR) var[k][d] = GMM_VAR_FLOOR;
      }
    }
#endif
  }
};


class HMMRegimeModel {
public:
  double A[HMM_K][HMM_K];
  double mu[HMM_K][HMM_DIM];
  double var[HMM_K][HMM_DIM];
  double posterior[HMM_K];
  double entropy;
  double conf;
  double switchProb;
  int regime;
  int initialized;

  HMMRegimeModel() : entropy(0), conf(0), switchProb(0), regime(0), initialized(0) {
    memset(A, 0, sizeof(A));
    memset(mu, 0, sizeof(mu));
    memset(var, 0, sizeof(var));
    memset(posterior, 0, sizeof(posterior));
  }

  void init() {
    for(int i=0;i<HMM_K;i++) {
      for(int j=0;j<HMM_K;j++) A[i][j] = (i==j) ? 0.90 : 0.10/(double)(HMM_K-1);
      for(int d=0; d<HMM_DIM; d++) {
        mu[i][d] = 0.03 * (i - 1);
        var[i][d] = 1.0;
      }
      posterior[i] = 1.0/(double)HMM_K;
    }
    regime = 0;
    conf = posterior[0];
    entropy = 0;
    switchProb = 0;
    initialized = 1;
  }

  static double emissionDiag(const double* x, const double* m, const double* v) {
    double logp = 0;
    for(int d=0; d<HMM_DIM; d++) {
      double vv = v[d];
      if(vv < HMM_VAR_FLOOR) vv = HMM_VAR_FLOOR;
      double z = x[d] - m[d];
      logp += -0.5 * (z*z / vv + log(vv + EPS));
    }
    if(logp < -80.0) logp = -80.0;
    return exp(logp);
  }

  void filter(const double obs[HMM_DIM]) {
    if(!initialized) init();

    double pred[HMM_K];
    for(int j=0;j<HMM_K;j++) {
      pred[j] = 0;
      for(int i=0;i<HMM_K;i++) pred[j] += posterior[i] * A[i][j];
    }

    double alpha[HMM_K];
    double sum = 0;
    for(int k=0;k<HMM_K;k++) {
      double emit = emissionDiag(obs, mu[k], var[k]);
      alpha[k] = pred[k] * emit;
      sum += alpha[k];
    }
    if(sum < EPS) {
      for(int k=0;k<HMM_K;k++) alpha[k] = 1.0/(double)HMM_K;
    } else {
      for(int k=0;k<HMM_K;k++) alpha[k] /= sum;
    }

    for(int k=0;k<HMM_K;k++) posterior[k] = alpha[k];

    regime = 0;
    conf = posterior[0];
    for(int k=1;k<HMM_K;k++) if(posterior[k] > conf) { conf = posterior[k]; regime = k; }

    entropy = 0;
    for(int k=0;k<HMM_K;k++) entropy -= posterior[k] * log(posterior[k] + EPS);

    switchProb = 1.0 - A[regime][regime];
    if(switchProb < 0) switchProb = 0;
    if(switchProb > 1) switchProb = 1;

#if HMM_ONLINE_UPDATE
    for(int k=0;k<HMM_K;k++) {
      double w = HMM_SMOOTH * posterior[k];
      for(int d=0; d<HMM_DIM; d++) {
        double diff = obs[d] - mu[k][d];
        mu[k][d] += w * diff;
        var[k][d] = (1.0 - w) * var[k][d] + w * diff * diff;
        if(var[k][d] < HMM_VAR_FLOOR) var[k][d] = HMM_VAR_FLOOR;
      }
    }
#endif
  }
};

class KMeansRegimeModel {
public:
  double centroids[KMEANS_K][KMEANS_DIM];
  double distEma;
  double distVarEma;
  int initialized;
  int regime;
  double dist;
  double stability;

  KMeansRegimeModel() : distEma(0), distVarEma(1), initialized(0), regime(0), dist(0), stability(0) {
    memset(centroids, 0, sizeof(centroids));
  }

  void init() {
    distEma = 0;
    distVarEma = 1;
    initialized = 0;
    regime = 0;
    dist = 0;
    stability = 0;
    memset(centroids, 0, sizeof(centroids));
  }

  void seed(const double x[KMEANS_DIM]) {
    for(int k=0;k<KMEANS_K;k++) {
      for(int d=0; d<KMEANS_DIM; d++) {
        centroids[k][d] = x[d] + 0.03 * (k - 1);
      }
    }
    initialized = 1;
  }

  static double clampRange(double x, double lo, double hi) {
    if(x < lo) return lo;
    if(x > hi) return hi;
    return x;
  }

  void predictAndUpdate(const double x[KMEANS_DIM]) {
    if(!initialized) seed(x);

    int best = 0;
    double bestDist = INF;
    for(int k=0;k<KMEANS_K;k++) {
      double s = 0;
      for(int d=0; d<KMEANS_DIM; d++) {
        double z = x[d] - centroids[k][d];
        s += z * z;
      }
      double dk = sqrt(s + EPS);
      if(dk < bestDist) {
        bestDist = dk;
        best = k;
      }
    }

    regime = best;
    dist = bestDist;

    distEma = (1.0 - KMEANS_DIST_EMA) * distEma + KMEANS_DIST_EMA * dist;
    double dd = dist - distEma;
    distVarEma = (1.0 - KMEANS_DIST_EMA) * distVarEma + KMEANS_DIST_EMA * dd * dd;
    double distStd = sqrt(distVarEma + EPS);
    double zDist = (dist - distEma) / (distStd + EPS);
    stability = clampRange(1.0 / (1.0 + exp(zDist)), 0.0, 1.0);

#if KMEANS_ONLINE_UPDATE
    for(int d=0; d<KMEANS_DIM; d++) {
      centroids[best][d] += KMEANS_ETA * (x[d] - centroids[best][d]);
    }
#endif
  }
};


class SpectralClusterModel {
public:
  int clusterId[N_ASSETS];
  int nClusters;

  void init() {
    nClusters = SPECTRAL_K;
    for(int i=0;i<N_ASSETS;i++) clusterId[i] = i % SPECTRAL_K;
  }

  void update(const fvar* distMatrix) {
    if(!distMatrix) return;
    // lightweight deterministic clustering surrogate from distance rows
    for(int i=0;i<N_ASSETS;i++) {
      double sig = 0;
      for(int j=0;j<N_ASSETS;j++) {
        if(i == j) continue;
        double d = (double)distMatrix[i*N_ASSETS + j];
        if(d < INF) sig += d;
      }
      int cid = (int)fmod(fabs(sig * 1000.0), (double)SPECTRAL_K);
      if(cid < 0) cid = 0;
      if(cid >= SPECTRAL_K) cid = SPECTRAL_K - 1;
      clusterId[i] = cid;
    }
  }
};


class HierarchicalClusteringModel {
public:
  int clusterCoarse[N_ASSETS];
  int clusterFine[N_ASSETS];
  int nCoarse;
  int nFine;

  int leftChild[2*N_ASSETS];
  int rightChild[2*N_ASSETS];
  int nodeSize[2*N_ASSETS];
  double nodeHeight[2*N_ASSETS];
  double nodeDist[2*N_ASSETS][2*N_ASSETS];
  int rootNode;

  void init() {
    nCoarse = HCLUST_COARSE_K;
    nFine = HCLUST_FINE_K;
    rootNode = N_ASSETS - 1;
    for(int i=0;i<N_ASSETS;i++) {
      clusterCoarse[i] = i % HCLUST_COARSE_K;
      clusterFine[i] = i % HCLUST_FINE_K;
    }
  }

  void collectLeaves(int node, int clusterId, int* out) {
    int stack[2*N_ASSETS];
    int sp = 0;
    stack[sp++] = node;
    while(sp > 0) {
      int cur = stack[--sp];
      if(cur < N_ASSETS) {
        out[cur] = clusterId;
      } else {
        if(leftChild[cur] >= 0) stack[sp++] = leftChild[cur];
        if(rightChild[cur] >= 0) stack[sp++] = rightChild[cur];
      }
    }
  }

  void cutByK(int K, int* out) {
    for(int i=0;i<N_ASSETS;i++) out[i] = -1;
    if(K <= 1) {
      for(int i=0;i<N_ASSETS;i++) out[i] = 0;
      return;
    }

    int clusters[2*N_ASSETS];
    int count = 1;
    clusters[0] = rootNode;

    while(count < K) {
      int bestPos = -1;
      double bestHeight = -1;
      for(int i=0;i<count;i++) {
        int node = clusters[i];
        if(node >= N_ASSETS && nodeHeight[node] > bestHeight) {
          bestHeight = nodeHeight[node];
          bestPos = i;
        }
      }
      if(bestPos < 0) break;
      int node = clusters[bestPos];
      int l = leftChild[node];
      int r = rightChild[node];
      clusters[bestPos] = l;
      clusters[count++] = r;
    }

    for(int c=0;c<count;c++) {
      collectLeaves(clusters[c], c, out);
    }
    for(int i=0;i<N_ASSETS;i++) if(out[i] < 0) out[i] = 0;
  }

  void update(const fvar* distMatrix) {
    if(!distMatrix) return;

    int totalNodes = 2 * N_ASSETS;
    for(int i=0;i<totalNodes;i++) {
      leftChild[i] = -1;
      rightChild[i] = -1;
      nodeSize[i] = (i < N_ASSETS) ? 1 : 0;
      nodeHeight[i] = 0;
      for(int j=0;j<totalNodes;j++) nodeDist[i][j] = INF;
    }

    for(int i=0;i<N_ASSETS;i++) {
      for(int j=0;j<N_ASSETS;j++) {
        if(i == j) nodeDist[i][j] = 0;
        else {
          double d = (double)distMatrix[i*N_ASSETS + j];
          if(d < 0 || d >= INF) d = 1.0;
          nodeDist[i][j] = d;
        }
      }
    }

    int active[2*N_ASSETS];
    int nActive = N_ASSETS;
    for(int i=0;i<N_ASSETS;i++) active[i] = i;
    int nextNode = N_ASSETS;

    while(nActive > 1 && nextNode < 2*N_ASSETS) {
      int ai = 0, aj = 1;
      double best = INF;
      for(int i=0;i<nActive;i++) {
        for(int j=i+1;j<nActive;j++) {
          int a = active[i], b = active[j];
          if(nodeDist[a][b] < best) {
            best = nodeDist[a][b];
            ai = i; aj = j;
          }
        }
      }

      int a = active[ai];
      int b = active[aj];
      int m = nextNode++;

      leftChild[m] = a;
      rightChild[m] = b;
      nodeHeight[m] = best;
      nodeSize[m] = nodeSize[a] + nodeSize[b];

      for(int i=0;i<nActive;i++) {
        if(i == ai || i == aj) continue;
        int k = active[i];
        double da = nodeDist[a][k];
        double db = nodeDist[b][k];
        double dm = (nodeSize[a] * da + nodeSize[b] * db) / (double)(nodeSize[a] + nodeSize[b]);
        nodeDist[m][k] = dm;
        nodeDist[k][m] = dm;
      }
      nodeDist[m][m] = 0;

      if(aj < ai) { int t=ai; ai=aj; aj=t; }
      for(int i=aj;i<nActive-1;i++) active[i] = active[i+1];
      nActive--;
      for(int i=ai;i<nActive-1;i++) active[i] = active[i+1];
      nActive--;
      active[nActive++] = m;
    }

    rootNode = active[0];

    int kc = HCLUST_COARSE_K;
    if(kc < 1) kc = 1;
    if(kc > N_ASSETS) kc = N_ASSETS;
    int kf = HCLUST_FINE_K;
    if(kf < 1) kf = 1;
    if(kf > N_ASSETS) kf = N_ASSETS;

    cutByK(kc, clusterCoarse);
    cutByK(kf, clusterFine);
    nCoarse = kc;
    nFine = kf;
  }
};


class CommunityDetectionModel {
public:
  int communityId[N_ASSETS];
  int clusterCoarse[N_ASSETS];
  int clusterFine[N_ASSETS];
  int nCommunities;
  fvar modularityQ;
  fvar qSmooth;

  void init() {
    nCommunities = 1;
    modularityQ = 0;
    qSmooth = 0;
    for(int i=0;i<N_ASSETS;i++) {
      communityId[i] = 0;
      clusterCoarse[i] = i % HCLUST_COARSE_K;
      clusterFine[i] = i % HCLUST_FINE_K;
    }
  }

  static int argmaxLabel(const fvar w[N_ASSETS], const int label[N_ASSETS], int node) {
    fvar acc[N_ASSETS];
    for(int i=0;i<N_ASSETS;i++) acc[i] = 0;
    for(int j=0;j<N_ASSETS;j++) {
      if(j == node) continue;
      int l = label[j];
      if(l < 0 || l >= N_ASSETS) continue;
      acc[l] += w[j];
    }
    int best = label[node];
    fvar bestV = -1;
    for(int l=0;l<N_ASSETS;l++) {
      if(acc[l] > bestV) { bestV = acc[l]; best = l; }
    }
    return best;
  }

  void update(const fvar* corrMatrix, const fvar* distMatrix) {
    if(!corrMatrix || !distMatrix) return;

    fvar W[N_ASSETS][N_ASSETS];
    fvar degree[N_ASSETS];
    int label[N_ASSETS];

    for(int i=0;i<N_ASSETS;i++) {
      degree[i] = 0;
      label[i] = i;
      for(int j=0;j<N_ASSETS;j++) {
        if(i == j) W[i][j] = 0;
        else {
          fvar w = (fvar)fabs((double)corrMatrix[i*N_ASSETS + j]);
          if(w < (fvar)COMM_W_MIN) w = 0;
          W[i][j] = w;
          degree[i] += w;
        }
      }
    }

    // Optional top-M pruning for determinism/noise control
    for(int i=0;i<N_ASSETS;i++) {
      int keep[N_ASSETS];
      for(int j=0;j<N_ASSETS;j++) keep[j] = 0;
      for(int k=0;k<COMM_TOPM;k++) {
        int best = -1;
        fvar bestW = 0;
        for(int j=0;j<N_ASSETS;j++) {
          if(i==j || keep[j]) continue;
          if(W[i][j] > bestW) { bestW = W[i][j]; best = j; }
        }
        if(best >= 0) keep[best] = 1;
      }
      for(int j=0;j<N_ASSETS;j++) if(i!=j && !keep[j]) W[i][j] = 0;
    }

    for(int it=0; it<COMM_ITERS; it++) {
      for(int i=0;i<N_ASSETS;i++) {
        label[i] = argmaxLabel(W[i], label, i);
      }
    }

    // compress labels
    int map[N_ASSETS];
    for(int i=0;i<N_ASSETS;i++) map[i] = -1;
    int nLab = 0;
    for(int i=0;i<N_ASSETS;i++) {
      int l = label[i];
      if(l < 0 || l >= N_ASSETS) l = 0;
      if(map[l] < 0) map[l] = nLab++;
      communityId[i] = map[l];
    }
    if(nLab < 1) nLab = 1;
    nCommunities = nLab;

    // modularity approximation
    fvar m2 = 0;
    for(int i=0;i<N_ASSETS;i++) for(int j=0;j<N_ASSETS;j++) m2 += W[i][j];
    if(m2 < (fvar)EPS) {
      modularityQ = 0;
    } else {
      fvar q = 0;
      for(int i=0;i<N_ASSETS;i++) {
        for(int j=0;j<N_ASSETS;j++) {
          if(communityId[i] == communityId[j]) {
            q += W[i][j] - (degree[i] * degree[j] / m2);
          }
        }
      }
      modularityQ = q / m2;
    }

    qSmooth = (fvar)(1.0 - COMM_Q_EMA) * qSmooth + (fvar)COMM_Q_EMA * modularityQ;

    for(int i=0;i<N_ASSETS;i++) {
      int c = communityId[i];
      if(c < 0) c = 0;
      clusterCoarse[i] = c % HCLUST_COARSE_K;
      clusterFine[i] = c % HCLUST_FINE_K;
    }
  }
};


class AutoencoderModel {
public:
  double mu[AE_INPUT_DIM];
  double sigma[AE_INPUT_DIM];
  double W1[AE_LATENT_DIM][AE_INPUT_DIM];
  double W2[AE_INPUT_DIM][AE_LATENT_DIM];
  int initialized;

  void init() {
    initialized = 1;
    for(int i=0;i<AE_INPUT_DIM;i++) {
      mu[i] = 0;
      sigma[i] = 1;
    }
    for(int z=0;z<AE_LATENT_DIM;z++) {
      for(int d=0;d<AE_INPUT_DIM;d++) {
        double w = sin((double)(z+1)*(d+1)) * 0.05;
        W1[z][d] = w;
        W2[d][z] = w;
      }
    }
  }

  static double act(double x) {
    if(x > 4) x = 4;
    if(x < -4) x = -4;
    return tanh(x);
  }

  double infer(const double xIn[AE_INPUT_DIM]) {
    if(!initialized) init();

    double x[AE_INPUT_DIM];
    for(int d=0;d<AE_INPUT_DIM;d++) x[d] = (xIn[d] - mu[d]) / (sigma[d] + EPS);

    double z[AE_LATENT_DIM];
    for(int k=0;k<AE_LATENT_DIM;k++) {
      double s = 0;
      for(int d=0;d<AE_INPUT_DIM;d++) s += W1[k][d] * x[d];
      z[k] = act(s);
    }

    double recon[AE_INPUT_DIM];
    for(int d=0;d<AE_INPUT_DIM;d++) {
      double s = 0;
      for(int k=0;k<AE_LATENT_DIM;k++) s += W2[d][k] * z[k];
      recon[d] = act(s);
    }

    double err = 0;
    for(int d=0;d<AE_INPUT_DIM;d++) {
      double e = x[d] - recon[d];
      err += e*e;
    }
    err /= (double)AE_INPUT_DIM;

    for(int d=0;d<AE_INPUT_DIM;d++) {
      mu[d] = (1.0 - AE_NORM_ALPHA) * mu[d] + AE_NORM_ALPHA * xIn[d];
      double dv = xIn[d] - mu[d];
      sigma[d] = (1.0 - AE_NORM_ALPHA) * sigma[d] + AE_NORM_ALPHA * sqrt(dv*dv + EPS);
      if(sigma[d] < 1e-5) sigma[d] = 1e-5;
    }
    return err;
  }
};

class NoveltyController {
public:
  double errEma;
  double errVar;
  double zRecon;
  int regime;
  double riskScale;

  void init() {
    errEma = 0;
    errVar = 1;
    zRecon = 0;
    regime = 0;
    riskScale = 1.0;
  }

  static double clampRange(double x, double lo, double hi) {
    if(x < lo) return lo;
    if(x > hi) return hi;
    return x;
  }

  void update(double reconError) {
    errEma = (1.0 - AE_ERR_EMA) * errEma + AE_ERR_EMA * reconError;
    double d = reconError - errEma;
    errVar = (1.0 - AE_ERR_EMA) * errVar + AE_ERR_EMA * d*d;
    double errStd = sqrt(errVar + EPS);
    zRecon = (reconError - errEma) / (errStd + EPS);

    if(zRecon >= AE_Z_HIGH) { regime = 2; riskScale = 0.20; }
    else if(zRecon >= AE_Z_LOW) { regime = 1; riskScale = 0.60; }
    else { regime = 0; riskScale = 1.00; }

    riskScale = clampRange(riskScale, 0.20, 1.00);
  }

  void apply(int* topK, double* scoreScale) {
    if(regime == 2) {
      if(*topK > 3) *topK -= 2;
      *scoreScale *= 0.60;
    } else if(regime == 1) {
      if(*topK > 3) *topK -= 1;
      *scoreScale *= 0.85;
    }
    if(*topK < 1) *topK = 1;
    if(*topK > TOP_K) *topK = TOP_K;
    *scoreScale = clampRange(*scoreScale, 0.10, 2.00);
  }
};


class SOMModel {
public:
  double W[SOM_H][SOM_W][SOM_DIM];
  int hitCount[SOM_H][SOM_W];
  int bmuX;
  int bmuY;
  double conf;
  int initialized;

  void init() {
    initialized = 1;
    bmuX = 0; bmuY = 0; conf = 0;
    for(int y=0;y<SOM_H;y++) {
      for(int x=0;x<SOM_W;x++) {
        hitCount[y][x] = 0;
        for(int d=0;d<SOM_DIM;d++) {
          W[y][x][d] = 0.02 * sin((double)(y+1)*(x+1)*(d+1));
        }
      }
    }
  }

  static double clampRange(double x,double lo,double hi){ if(x<lo) return lo; if(x>hi) return hi; return x; }

  void inferOrUpdate(const double s[SOM_DIM], int step) {
    if(!initialized) init();

    int bx=0, by=0;
    double best=INF, second=INF;
    for(int y=0;y<SOM_H;y++) {
      for(int x=0;x<SOM_W;x++) {
        double d2=0;
        for(int k=0;k<SOM_DIM;k++) {
          double z = s[k] - W[y][x][k];
          d2 += z*z;
        }
        if(d2 < best) { second = best; best=d2; bx=x; by=y; }
        else if(d2 < second) second = d2;
      }
    }

    bmuX = bx; bmuY = by;
    double d1 = sqrt(best + EPS);
    double d2 = sqrt(second + EPS);
    conf = clampRange((d2 - d1) / (d2 + EPS), 0.0, 1.0);
    hitCount[bmuY][bmuX]++;

#if SOM_ONLINE_UPDATE
    double alpha = SOM_ALPHA_MIN + (SOM_ALPHA_MAX - SOM_ALPHA_MIN) * exp(-0.005 * step);
    double sigma = SOM_SIGMA_MIN + (SOM_SIGMA_MAX - SOM_SIGMA_MIN) * exp(-0.005 * step);
    for(int y=0;y<SOM_H;y++) {
      for(int x=0;x<SOM_W;x++) {
        double gd2 = (double)((x-bmuX)*(x-bmuX) + (y-bmuY)*(y-bmuY));
        double h = exp(-gd2 / (2.0*sigma*sigma + EPS));
        for(int k=0;k<SOM_DIM;k++) {
          W[y][x][k] += alpha * h * (s[k] - W[y][x][k]);
        }
      }
    }
#endif
  }

  int regimeId() const { return bmuY * SOM_W + bmuX; }
};

class SOMPlaybook {
public:
  int region;
  double riskScale;

  void init() { region = 0; riskScale = 1.0; }

  void apply(const SOMModel& som, int* topK, double* scoreScale) {
    int mx = som.bmuX;
    int my = som.bmuY;
    int cx = (mx >= SOM_W/2) ? 1 : 0;
    int cy = (my >= SOM_H/2) ? 1 : 0;
    region = cy * 2 + cx;

    if(region == 0) { *scoreScale *= 1.02; riskScale = 1.00; }
    else if(region == 1) { *scoreScale *= 0.95; riskScale = 0.85; if(*topK > 3) (*topK)--; }
    else if(region == 2) { *scoreScale *= 0.90; riskScale = 0.70; if(*topK > 3) (*topK)--; }
    else { *scoreScale *= 0.80; riskScale = 0.50; if(*topK > 2) (*topK)-=2; }

    if(som.conf < SOM_CONF_MIN) {
      riskScale *= 0.8;
      if(*topK > 2) (*topK)--;
    }

    if(*topK < 1) *topK = 1;
    if(*topK > TOP_K) *topK = TOP_K;
    if(*scoreScale < 0.10) *scoreScale = 0.10;
    if(*scoreScale > 2.00) *scoreScale = 2.00;
  }
};


class DensityStateBuilder {
public:
  double mean[DENSITY_DIM];
  double var[DENSITY_DIM];
  int initialized;

  void init() {
    initialized = 0;
    for(int d=0; d<DENSITY_DIM; d++) {
      mean[d] = 0;
      var[d] = 1;
    }
  }

  static double clampRange(double x, double lo, double hi) {
    if(x < lo) return lo;
    if(x > hi) return hi;
    return x;
  }

  void build(const LearningSnapshot& s, int updateCount, int openclReady, double scoreScale, int dynamicTopK, double out[DENSITY_DIM]) {
    double raw[DENSITY_DIM];
    raw[0] = s.meanScore;
    raw[1] = s.meanCompactness;
    raw[2] = s.meanVol;
    raw[3] = s.meanAbsCorr;
    raw[4] = s.stdAbsCorr;
    raw[5] = s.momentumMean;
    raw[6] = scoreScale;
    raw[7] = (double)dynamicTopK / (double)(TOP_K > 0 ? TOP_K : 1);
    raw[8] = (double)openclReady;
    raw[9] = (double)updateCount / 1000.0;

    if(!initialized) {
      for(int d=0; d<DENSITY_DIM; d++) {
        mean[d] = raw[d];
        var[d] = 1.0;
      }
      initialized = 1;
    }

    for(int d=0; d<DENSITY_DIM; d++) {
      mean[d] = (1.0 - DENSITY_Z_EMA) * mean[d] + DENSITY_Z_EMA * raw[d];
      double diff = raw[d] - mean[d];
      var[d] = (1.0 - DENSITY_Z_EMA) * var[d] + DENSITY_Z_EMA * diff * diff;
      double z = diff / (sqrt(var[d] + EPS) + EPS);
      out[d] = clampRange(z, -4.0, 4.0);
    }
  }
};

class DensityModel {
public:
  double hist[DENSITY_HISTORY][DENSITY_DIM];
  int histCount;
  int histWrite;
  double proto[DENSITY_N_REGIMES][DENSITY_DIM];
  int protoN[DENSITY_N_REGIMES];
  int protoCount;

  void init() {
    histCount = 0;
    histWrite = 0;
    protoCount = 0;
    memset(hist, 0, sizeof(hist));
    memset(proto, 0, sizeof(proto));
    memset(protoN, 0, sizeof(protoN));
  }

  static double dist2(const double* a, const double* b) {
    double s = 0;
    for(int d=0; d<DENSITY_DIM; d++) {
      double x = a[d] - b[d];
      s += x * x;
    }
    return s;
  }

  void observe(const double z[DENSITY_DIM], int updateCount, int* labelOut, double* confOut, int* noiseOut) {
    int neighbors = 0;
    double eps2 = DENSITY_EPS * DENSITY_EPS;
    for(int i=0; i<histCount; i++) {
      if(dist2(hist[i], z) <= eps2) neighbors++;
    }

    int isNoise = (neighbors < DENSITY_MIN_SAMPLES) ? 1 : 0;
    int label = -1;
    double conf = 0;

    if(isNoise) {
      conf = (double)neighbors / (double)(DENSITY_MIN_SAMPLES + 1);
    } else {
      if(protoCount < DENSITY_N_REGIMES && ((updateCount % 3) == 0 || protoCount == 0)) {
        label = protoCount;
        for(int d=0; d<DENSITY_DIM; d++) proto[label][d] = z[d];
        protoN[label] = 1;
        protoCount++;
      } else {
        int best = 0;
        double bestD = INF;
        int useCount = (protoCount > 0) ? protoCount : 1;
        for(int k=0; k<useCount; k++) {
          double dk = dist2(z, proto[k]);
          if(dk < bestD) {
            bestD = dk;
            best = k;
          }
        }
        label = best;
        protoN[best]++;
        double lr = 1.0 / (double)protoN[best];
        for(int d=0; d<DENSITY_DIM; d++) {
          proto[best][d] += lr * (z[d] - proto[best][d]);
        }
        conf = 1.0 / (1.0 + sqrt(bestD + EPS));
      }
    }

    for(int d=0; d<DENSITY_DIM; d++) hist[histWrite][d] = z[d];
    histWrite = (histWrite + 1) % DENSITY_HISTORY;
    if(histCount < DENSITY_HISTORY) histCount++;

    *labelOut = label;
    *confOut = conf;
    *noiseOut = isNoise;
  }
};

class UnknownRegimeController {
public:
  int noiseStreak;
  int clearStreak;
  int paused;
  double appliedRiskScale;

  void init() {
    noiseStreak = 0;
    clearStreak = 0;
    paused = 0;
    appliedRiskScale = 1.0;
  }

  static double clampRange(double x, double lo, double hi) {
    if(x < lo) return lo;
    if(x > hi) return hi;
    return x;
  }

  void apply(int isNoise, int label, double conf, int strategyProfile, int* topK, double* scoreScale, double* riskScale) {
    (void)label;
    (void)conf;
    const double profileRisk[5] = {0.52, 0.48, 0.56, 0.50, 0.54};
    const int profileSpread[5] = {1, 2, 1, 2, 1};
    int p = strategyProfile;
    if(p < 0) p = 0;
    if(p > 4) p = 4;

    if(isNoise) {
      noiseStreak++;
      clearStreak = 0;
    } else {
      clearStreak++;
      if(noiseStreak > 0) noiseStreak--;
    }

    if(noiseStreak >= DENSITY_NOISE_PERSIST) paused = 1;
    if(clearStreak >= DENSITY_CLEAR_PERSIST) paused = 0;

    appliedRiskScale = 1.0;

    if(paused) {
      appliedRiskScale = profileRisk[p];
      *scoreScale *= 0.70;
      int widened = *topK + profileSpread[p];
      if(widened > TOP_K) widened = TOP_K;
      *topK = widened;
    } else if(isNoise) {
      appliedRiskScale = 0.80;
      *scoreScale *= 0.90;
      int widened = *topK + 1;
      if(widened > TOP_K) widened = TOP_K;
      *topK = widened;
    }

    if(*topK < 1) *topK = 1;
    if(*topK > TOP_K) *topK = TOP_K;
    *scoreScale = clampRange(*scoreScale, 0.10, 2.00);
    *riskScale = clampRange((*riskScale) * appliedRiskScale, 0.20, 1.00);
  }
};

class StrategyController {
public:
  UnsupervisedModel unsup;
  RLAgent rl;
  PCAModel pca;
  GMMRegimeModel gmm;
  HMMRegimeModel hmm;
  KMeansRegimeModel kmeans;
  int dynamicTopK;
  double scoreScale;
  int regime;
  double adaptiveGamma;
  double adaptiveAlpha;
  double adaptiveBeta;
  double adaptiveLambda;
  double riskScale;
  int cooldown;

  StrategyController()
  : dynamicTopK(TOP_K), scoreScale(1.0), regime(0),
    adaptiveGamma(1.0), adaptiveAlpha(1.0), adaptiveBeta(1.0), adaptiveLambda(1.0), riskScale(1.0), cooldown(0) {}

  static double clampRange(double x, double lo, double hi) {
    if(x < lo) return lo;
    if(x > hi) return hi;
    return x;
  }

  void init() {
    unsup.init();
    rl.init();
    pca.init();
    gmm.init();
    hmm.init();
    kmeans.init();
    dynamicTopK = TOP_K;
    scoreScale = 1.0;
    regime = 0;
    adaptiveGamma = 1.0;
    adaptiveAlpha = 1.0;
    adaptiveBeta = 1.0;
    adaptiveLambda = 1.0;
    riskScale = 1.0;
    cooldown = 0;
  }

  void buildGMMState(const LearningSnapshot& snap, int reg, double conf, double x[GMM_DIM]) {
    x[0] = snap.meanScore;
    x[1] = snap.meanCompactness;
    x[2] = snap.meanVol;
    x[3] = pca.dom;
    x[4] = pca.rot;
    x[5] = (double)reg / 2.0;
    x[6] = conf;
    x[7] = snap.meanScore - snap.meanCompactness;
  }

  void buildHMMObs(const LearningSnapshot& snap, int reg, double conf, double x[HMM_DIM]) {
    x[0] = pca.latent[0];
    x[1] = pca.latent[1];
    x[2] = pca.latent[2];
    x[3] = snap.meanVol;
    x[4] = snap.meanScore;
    x[5] = snap.meanCompactness;
    x[6] = (double)reg / 2.0;
    x[7] = conf;
  }

  void buildKMeansState(const LearningSnapshot& snap, int reg, double conf, double x[KMEANS_DIM]) {
    x[0] = pca.latent[0];
    x[1] = pca.latent[1];
    x[2] = pca.latent[2];
    x[3] = snap.meanVol;
    x[4] = snap.meanScore;
    x[5] = snap.meanCompactness;
    x[6] = (double)reg / 2.0;
    x[7] = conf;
  }

  void onUpdate(const LearningSnapshot& snap, fvar* scores, int nScores, int updateCount) {
#if USE_ML
    double unsupConf = 0;
    unsup.update(snap, &regime, &unsupConf);
#if USE_PCA
    pca.update(snap, regime, unsupConf);
#else
    pca.dom = 0.5;
    pca.rot = 0.0;
#endif

#if USE_GMM
    double gx[GMM_DIM];
    buildGMMState(snap, regime, unsupConf, gx);
    gmm.infer(gx);
#if USE_HMM
    double hx[HMM_DIM];
    buildHMMObs(snap, regime, unsupConf, hx);
    hmm.filter(hx);
#if USE_KMEANS
    double kx[KMEANS_DIM];
    buildKMeansState(snap, regime, unsupConf, kx);
    kmeans.predictAndUpdate(kx);
#endif
#endif
    // regime presets: [gamma, alpha, beta, lambda]
    const double presets[GMM_K][4] = {
      {1.05, 1.00, 0.95, 1.00},
      {0.95, 1.05, 1.05, 0.95},
      {1.00, 0.95, 1.10, 1.05}
    };
    adaptiveGamma = 0;
    adaptiveAlpha = 0;
    adaptiveBeta  = 0;
    adaptiveLambda = 0;
    for(int k=0;k<GMM_K;k++) {
#if USE_HMM
      adaptiveGamma += hmm.posterior[k] * presets[k][0];
      adaptiveAlpha += hmm.posterior[k] * presets[k][1];
      adaptiveBeta  += hmm.posterior[k] * presets[k][2];
      adaptiveLambda += hmm.posterior[k] * presets[k][3];
#else
      adaptiveGamma += gmm.p[k] * presets[k][0];
      adaptiveAlpha += gmm.p[k] * presets[k][1];
      adaptiveBeta  += gmm.p[k] * presets[k][2];
      adaptiveLambda += gmm.p[k] * presets[k][3];
#endif
    }
#if USE_HMM
    double entNorm = hmm.entropy / log((double)HMM_K + EPS);
    riskScale = clampRange(1.0 - 0.45 * entNorm, HMM_MIN_RISK, 1.0);
    if(hmm.entropy > HMM_ENTROPY_TH || hmm.switchProb > HMM_SWITCH_TH) cooldown = HMM_COOLDOWN_UPDATES;
    else if(cooldown > 0) cooldown--;
#else
    double entNorm = gmm.entropy / log((double)GMM_K + EPS);
    riskScale = clampRange(1.0 - GMM_ENTROPY_COEFF * entNorm, GMM_MIN_RISK, 1.0);
#endif
#else
    adaptiveGamma = 1.0 + 0.35 * pca.dom - 0.25 * pca.rot;
    adaptiveAlpha = 1.0 + 0.30 * pca.dom;
    adaptiveBeta  = 1.0 + 0.25 * pca.rot;
    adaptiveLambda = 1.0 + 0.20 * pca.dom - 0.20 * pca.rot;
    riskScale = 1.0;
#endif

    adaptiveGamma = clampRange(adaptiveGamma, 0.80, 1.40);
    adaptiveAlpha = clampRange(adaptiveAlpha, 0.85, 1.35);
    adaptiveBeta  = clampRange(adaptiveBeta, 0.85, 1.35);
    adaptiveLambda = clampRange(adaptiveLambda, 0.85, 1.25);

#if USE_KMEANS
    const double kmPreset[KMEANS_K][4] = {
      {1.02, 1.00, 0.98, 1.00},
      {1.08, 0.96, 0.95, 1.02},
      {0.94, 1.08, 1.08, 0.92}
    };
    int kr = kmeans.regime;
    if(kr < 0) kr = 0;
    if(kr >= KMEANS_K) kr = KMEANS_K - 1;
    double wkm = clampRange(kmeans.stability, 0.0, 1.0);
    adaptiveGamma = (1.0 - wkm) * adaptiveGamma + wkm * kmPreset[kr][0];
    adaptiveAlpha = (1.0 - wkm) * adaptiveAlpha + wkm * kmPreset[kr][1];
    adaptiveBeta  = (1.0 - wkm) * adaptiveBeta  + wkm * kmPreset[kr][2];
    adaptiveLambda = (1.0 - wkm) * adaptiveLambda + wkm * kmPreset[kr][3];
    if(kmeans.stability < KMEANS_STABILITY_MIN) {
      riskScale *= 0.85;
      if(cooldown < 1) cooldown = 1;
    }
#endif

    rl.updateReward(snap.meanScore);
    rl.lastAction = rl.chooseAction(updateCount);

    int baseTopK = TOP_K;
    if(rl.lastAction == 0) baseTopK = TOP_K - 2;
    else if(rl.lastAction == 1) baseTopK = TOP_K;
    else if(rl.lastAction == 2) baseTopK = TOP_K;
    else baseTopK = TOP_K - 1;

    double profileBias[5] = {1.00, 0.98, 0.99, 0.97, 1.02};
    scoreScale = (1.0 + 0.06 * (adaptiveGamma - 1.0) + 0.04 * (adaptiveAlpha - 1.0) - 0.04 * (adaptiveBeta - 1.0))
               * profileBias[STRATEGY_PROFILE] * riskScale;

    if(pca.dom > 0.60) baseTopK -= 1;
    if(pca.rot > 0.15) baseTopK -= 1;
#if USE_HMM
    if(hmm.regime == 2) baseTopK -= 1;
    if(cooldown > 0) baseTopK -= 1;
#if USE_KMEANS
    if(kmeans.regime == 2) baseTopK -= 1;
#endif
#elif USE_GMM
    if(gmm.bestRegime == 2) baseTopK -= 1;
#endif

    dynamicTopK = baseTopK;
    if(dynamicTopK < 1) dynamicTopK = 1;
    if(dynamicTopK > TOP_K) dynamicTopK = TOP_K;

    for(int i=0; i<nScores; i++) {
      double s = (double)scores[i] * scoreScale;
      if(s > 1.0) s = 1.0;
      if(s < 0.0) s = 0.0;
      scores[i] = (fvar)s;
    }
#else
    (void)snap; (void)scores; (void)nScores; (void)updateCount;
#endif
  }
};

// ---------------------------- Strategy ----------------------------

class CompactDominantStrategy {
public:
  ExposureTable exposureTable;
  FeatureBufferSoA featSoA;
  OpenCLBackend openCL;

  SlabAllocator<fvar> corrMatrix;
  SlabAllocator<fvar> distMatrix;
  SlabAllocator<fvar> compactness;
  SlabAllocator<fvar> scores;

  SlabAllocator<float> featLinear;
  SlabAllocator<float> corrLinear;

  int barCount;
  int updateCount;
  StrategyController controller;
  HierarchicalClusteringModel hclust;
  CommunityDetectionModel comm;
  AutoencoderModel ae;
  NoveltyController novelty;
  SOMModel som;
  SOMPlaybook somPlaybook;
  DensityStateBuilder densityState;
  DensityModel densityModel;
  UnknownRegimeController densityRegimeCtrl;
  int densityLabel;
  double densityConf;
  int densityNoise;

  CompactDominantStrategy() : barCount(0), updateCount(0), densityLabel(-1), densityConf(0.0), densityNoise(0) {}

  void init() {
    printf("CompactDominant_v14: Initializing...\n");

    exposureTable.init();
    featSoA.init(N_ASSETS, FEAT_WINDOW);

    corrMatrix.init(N_ASSETS * N_ASSETS);
    distMatrix.init(N_ASSETS * N_ASSETS);
    compactness.init(N_ASSETS);
    scores.init(N_ASSETS);

    featLinear.init(FEAT_N * N_ASSETS * FEAT_WINDOW);
    corrLinear.init(N_ASSETS * N_ASSETS);

    openCL.init();
    printf("CompactDominant_v14: Ready (OpenCL=%d)\n", openCL.ready);
    controller.init();
    hclust.init();
    comm.init();
    ae.init();
    novelty.init();
    som.init();
    somPlaybook.init();
    densityState.init();
    densityModel.init();
    densityRegimeCtrl.init();
    densityLabel = -1;
    densityConf = 0.0;
    densityNoise = 0;

    barCount = 0;
    updateCount = 0;
  }

  void shutdown() {
    printf("CompactDominant_v14: Shutting down...\n");

    openCL.shutdown();

    featSoA.shutdown();
    corrMatrix.shutdown();
    distMatrix.shutdown();
    compactness.shutdown();
    scores.shutdown();

    featLinear.shutdown();
    corrLinear.shutdown();
  }

  void computeFeatures(int assetIdx) {
    asset((char*)ASSET_NAMES[assetIdx]);

    vars C = series(priceClose(0));
    vars V = series(Volatility(C, 20));

    if(Bar < 50) return;

    fvar r1 = (fvar)log(C[0] / C[1]);
    fvar rN = (fvar)log(C[0] / C[12]);
    fvar vol = (fvar)V[0];
    fvar zscore = (fvar)((C[0] - C[50]) / (V[0] * 20.0 + EPS));
    fvar rangeP = (fvar)((C[0] - C[50]) / (C[0] + EPS));
    fvar flow = (fvar)(r1 * vol);
    fvar regime = (fvar)((vol > 0.001) ? 1.0 : 0.0);
    fvar volOfVol = (fvar)(vol * vol);
    fvar persistence = (fvar)fabs(r1);

    featSoA.push(0, assetIdx, r1);
    featSoA.push(1, assetIdx, rN);
    featSoA.push(2, assetIdx, vol);
    featSoA.push(3, assetIdx, zscore);
    featSoA.push(4, assetIdx, rangeP);
    featSoA.push(5, assetIdx, flow);
    featSoA.push(6, assetIdx, regime);
    featSoA.push(7, assetIdx, volOfVol);
    featSoA.push(8, assetIdx, persistence);
  }

  void computeCorrelationMatrixCPU() {
    for(int i=0;i<N_ASSETS*N_ASSETS;i++) corrMatrix[i] = 0;

    for(int f=0; f<FEAT_N; f++){
      for(int a=0; a<N_ASSETS; a++){
        for(int b=a+1; b<N_ASSETS; b++){
          fvar mx = 0, my = 0;
          for(int t=0; t<FEAT_WINDOW; t++){
            mx += featSoA.get(f,a,t);
            my += featSoA.get(f,b,t);
          }
          mx /= (fvar)FEAT_WINDOW;
          my /= (fvar)FEAT_WINDOW;

          fvar sxx = 0, syy = 0, sxy = 0;
          for(int t=0; t<FEAT_WINDOW; t++){
            fvar dx = featSoA.get(f,a,t) - mx;
            fvar dy = featSoA.get(f,b,t) - my;
            sxx += dx*dx;
            syy += dy*dy;
            sxy += dx*dy;
          }

          fvar den = (fvar)sqrt((double)(sxx*syy + (fvar)EPS));
          fvar corr = 0;
          if(den > (fvar)EPS) corr = sxy / den;
          else corr = 0;

          int idx = a*N_ASSETS + b;
          corrMatrix[idx] += corr / (fvar)FEAT_N;
          corrMatrix[b*N_ASSETS + a] = corrMatrix[idx];
        }
      }
    }
  }

  void buildFeatLinear() {
    int idx = 0;
    for(int f=0; f<FEAT_N; f++){
      for(int a=0; a<N_ASSETS; a++){
        for(int t=0; t<FEAT_WINDOW; t++){
          featLinear[idx] = (float)featSoA.get(f, a, t);
          idx++;
        }
      }
    }
  }

  void computeCorrelationMatrix() {
    if(openCL.ready) {
      buildFeatLinear();

      for(int i=0;i<N_ASSETS*N_ASSETS;i++) corrLinear[i] = 0.0f;

      int ok = openCL.computeCorrelationMatrixCL(
        featLinear.data,
        corrLinear.data,
        N_ASSETS,
        FEAT_N,
        FEAT_WINDOW
      );

      if(ok) {
        for(int i=0;i<N_ASSETS*N_ASSETS;i++) corrMatrix[i] = (fvar)0;

        for(int a=0; a<N_ASSETS; a++){
          corrMatrix[a*N_ASSETS + a] = (fvar)1.0;
          for(int b=a+1; b<N_ASSETS; b++){
            float c = corrLinear[a*N_ASSETS + b];
            corrMatrix[a*N_ASSETS + b] = (fvar)c;
            corrMatrix[b*N_ASSETS + a] = (fvar)c;
          }
        }
        return;
      }

      printf("OpenCL: runtime fail -> CPU fallback\n");
      openCL.ready = 0;
    }

    computeCorrelationMatrixCPU();
  }

  void computeDistanceMatrix() {
    for(int i=0;i<N_ASSETS;i++){
      for(int j=0;j<N_ASSETS;j++){
        if(i == j) {
          distMatrix[i*N_ASSETS + j] = (fvar)0;
        } else {
          fvar corrDist = (fvar)1.0 - (fvar)fabs((double)corrMatrix[i*N_ASSETS + j]);
          fvar expDist  = (fvar)exposureTable.getDist(i, j);
          fvar blended = (fvar)LAMBDA_META * corrDist + (fvar)(1.0 - (double)LAMBDA_META) * expDist;
          distMatrix[i*N_ASSETS + j] = blended;
        }
      }
    }
  }

  void floydWarshall() {
    fvar d[28][28];

    for(int i=0;i<N_ASSETS;i++){
      for(int j=0;j<N_ASSETS;j++){
        d[i][j] = distMatrix[i*N_ASSETS + j];
        if(i == j) d[i][j] = (fvar)0;
        if(d[i][j] < (fvar)0) d[i][j] = (fvar)INF;
      }
    }

    for(int k=0;k<N_ASSETS;k++){
      for(int i=0;i<N_ASSETS;i++){
        for(int j=0;j<N_ASSETS;j++){
          if(d[i][k] < (fvar)INF && d[k][j] < (fvar)INF) {
            fvar nk = d[i][k] + d[k][j];
            if(nk < d[i][j]) d[i][j] = nk;
          }
        }
      }
    }

    for(int i=0;i<N_ASSETS;i++){
      fvar w = 0;
      for(int j=i+1;j<N_ASSETS;j++){
        if(d[i][j] < (fvar)INF) w += d[i][j];
      }
      if(w > (fvar)0) compactness[i] = (fvar)(1.0 / (1.0 + (double)w));
      else compactness[i] = (fvar)0;
    }
  }

  void computeScores() {
    for(int i=0;i<N_ASSETS;i++){
      fvar coupling = 0;
      int count = 0;

      for(int j=0;j<N_ASSETS;j++){
        if(i != j && distMatrix[i*N_ASSETS + j] < (fvar)INF) {
          coupling += compactness[j];
          count++;
        }
      }

      fvar pCouple = 0;
      if(count > 0) pCouple = coupling / (fvar)count;
      else pCouple = (fvar)0;

      fvar regime = featSoA.get(6, i, 0);
      fvar rawScore = (fvar)ALPHA * regime + (fvar)GAMMA * compactness[i] - (fvar)BETA * pCouple;

      if(rawScore > (fvar)30) rawScore = (fvar)30;
      if(rawScore < (fvar)-30) rawScore = (fvar)-30;

      scores[i] = (fvar)(1.0 / (1.0 + exp(-(double)rawScore)));
    }
  }

  LearningSnapshot buildSnapshot() {
    LearningSnapshot s;
    s.meanScore = 0;
    s.meanCompactness = 0;
    s.meanVol = 0;
    for(int i=0;i<N_ASSETS;i++) {
      s.meanScore += (double)scores[i];
      s.meanCompactness += (double)compactness[i];
      s.meanVol += (double)featSoA.get(2, i, 0);
    }
    s.meanScore /= (double)N_ASSETS;
    s.meanCompactness /= (double)N_ASSETS;
    s.meanVol /= (double)N_ASSETS;

    double sumAbs = 0;
    double sumSq = 0;
    int nCorr = 0;
    for(int i=0;i<N_ASSETS;i++) {
      for(int j=i+1;j<N_ASSETS;j++) {
        double c = fabs((double)corrMatrix[i*N_ASSETS + j]);
        sumAbs += c;
        sumSq += c * c;
        nCorr++;
      }
    }
    if(nCorr > 0) {
      s.meanAbsCorr = sumAbs / (double)nCorr;
      double m2 = sumSq / (double)nCorr - s.meanAbsCorr * s.meanAbsCorr;
      if(m2 < 0) m2 = 0;

The Cartographer of Compact Currents (cont.) [Re: TipmyPip] #489356
04/08/26 17:08
04/08/26 17:08
Joined: Sep 2017
Posts: 293
TipmyPip Offline OP
Member
TipmyPip  Offline OP
Member

Joined: Sep 2017
Posts: 293

Code
  LearningSnapshot buildSnapshot() {
    LearningSnapshot s;
    s.meanScore = 0;
    s.meanCompactness = 0;
    s.meanVol = 0;
    for(int i=0;i<N_ASSETS;i++) {
      s.meanScore += (double)scores[i];
      s.meanCompactness += (double)compactness[i];
      s.meanVol += (double)featSoA.get(2, i, 0);
    }
    s.meanScore /= (double)N_ASSETS;
    s.meanCompactness /= (double)N_ASSETS;
    s.meanVol /= (double)N_ASSETS;

    double sumAbs = 0;
    double sumSq = 0;
    int nCorr = 0;
    for(int i=0;i<N_ASSETS;i++) {
      for(int j=i+1;j<N_ASSETS;j++) {
        double c = fabs((double)corrMatrix[i*N_ASSETS + j]);
        sumAbs += c;
        sumSq += c * c;
        nCorr++;
      }
    }
    if(nCorr > 0) {
      s.meanAbsCorr = sumAbs / (double)nCorr;
      double m2 = sumSq / (double)nCorr - s.meanAbsCorr * s.meanAbsCorr;
      if(m2 < 0) m2 = 0;
      s.stdAbsCorr = sqrt(m2);
    } else {
      s.meanAbsCorr = 0;
      s.stdAbsCorr = 0;
    }

    s.momentumMean = 0;
    for(int i=0;i<N_ASSETS;i++) s.momentumMean += (double)featSoA.get(1, i, 0);
    s.momentumMean /= (double)N_ASSETS;

    s.regime = 0;
    s.regimeConfidence = 0;
    return s;
  }

  void onBar() {
    barCount++;

    for(int i=0;i<N_ASSETS;i++) computeFeatures(i);

    if(barCount % UPDATE_EVERY == 0) {
      updateCount++;

      computeCorrelationMatrix();
      computeDistanceMatrix();
#if USE_COMMUNITY
      hclust.update(distMatrix.data);
#endif
#if USE_COMMUNITY
      comm.update(corrMatrix.data, distMatrix.data);
#endif
      floydWarshall();
      computeScores();
      LearningSnapshot snap = buildSnapshot();
      controller.onUpdate(snap, scores.data, N_ASSETS, updateCount);
#if USE_DENSITY
      double dz[DENSITY_DIM];
      densityState.build(snap, updateCount, openCL.ready, controller.scoreScale, controller.dynamicTopK, dz);
      densityModel.observe(dz, updateCount, &densityLabel, &densityConf, &densityNoise);
      densityRegimeCtrl.apply(densityNoise, densityLabel, densityConf, STRATEGY_PROFILE, &controller.dynamicTopK, &controller.scoreScale, &controller.riskScale);
      for(int i=0;i<N_ASSETS;i++) {
        double s = (double)scores[i] * densityRegimeCtrl.appliedRiskScale;
        if(s > 1.0) s = 1.0;
        if(s < 0.0) s = 0.0;
        scores[i] = (fvar)s;
      }
#endif
#if USE_AE
      double aeState[AE_INPUT_DIM];
      double ms=0, mc=0, mv=0;
      for(int i=0;i<N_ASSETS;i++){ ms += (double)scores[i]; mc += (double)compactness[i]; mv += (double)featSoA.get(2, i, 0); }
      ms /= (double)N_ASSETS; mc /= (double)N_ASSETS; mv /= (double)N_ASSETS;
      aeState[0] = ms;
      aeState[1] = mc;
      aeState[2] = mv;
      aeState[3] = controller.scoreScale;
      aeState[4] = (double)controller.dynamicTopK;
      aeState[5] = (double)barCount / (double)(LookBack + 1);
      aeState[6] = (double)updateCount / 1000.0;
      aeState[7] = (double)openCL.ready;
      double reconErr = ae.infer(aeState);
      novelty.update(reconErr);
      novelty.apply(&controller.dynamicTopK, &controller.scoreScale);
      for(int i=0;i<N_ASSETS;i++){{
        double s = (double)scores[i] * novelty.riskScale;
        if(s > 1.0) s = 1.0;
        if(s < 0.0) s = 0.0;
        scores[i] = (fvar)s;
      }}
#endif
      printTopK();
    }
  }

  void printTopK() {
    int indices[N_ASSETS];
    for(int i=0;i<N_ASSETS;i++) indices[i] = i;

    int topN = controller.dynamicTopK;
#if USE_COMMUNITY
    if(comm.qSmooth < (fvar)COMM_Q_LOW && topN > 2) topN--;
    if(comm.qSmooth > (fvar)COMM_Q_HIGH && topN < TOP_K) topN++;
#endif
    for(int i=0;i<topN;i++){
      for(int j=i+1;j<N_ASSETS;j++){
        if(scores[indices[j]] > scores[indices[i]]) {
          int tmp = indices[i];
          indices[i] = indices[j];
          indices[j] = tmp;
        }
      }
    }

    if(updateCount % 10 == 0) {
      printf("===CompactDominant_v14 Top-K(update#%d,OpenCL=%d)===\n",
        updateCount, openCL.ready);
#if USE_COMMUNITY
      printf(" communities=%d Q=%.4f\n", comm.nCommunities, (double)comm.qSmooth);
#endif
#if USE_DENSITY
      printf(" density(label=%d,noise=%d,conf=%.3f,pause=%d)\n", densityLabel, densityNoise, densityConf, densityRegimeCtrl.paused);
#endif

      int selected[N_ASSETS];
      int selCount = 0;
#if USE_COMMUNITY
      int coarseUsed[HCLUST_COARSE_K];
      int fineTake[HCLUST_FINE_K];
      int fineCap = (topN + HCLUST_FINE_K - 1) / HCLUST_FINE_K;
      for(int c=0;c<HCLUST_COARSE_K;c++) coarseUsed[c] = 0;
      for(int c=0;c<HCLUST_FINE_K;c++) fineTake[c] = 0;

      for(int i=0;i<topN;i++){
        int idx = indices[i];
        int cid = comm.clusterCoarse[idx];
        if(cid < 0 || cid >= HCLUST_COARSE_K) cid = 0;
        if(coarseUsed[cid]) continue;
        coarseUsed[cid] = 1;
        selected[selCount++] = idx;
        int fid = comm.clusterFine[idx];
        if(fid < 0 || fid >= HCLUST_FINE_K) fid = 0;
        fineTake[fid]++;
      }

      for(int i=0;i<topN && selCount<topN;i++){
        int idx = indices[i];
        int dup = 0;
        for(int k=0;k<selCount;k++) if(selected[k]==idx){ dup=1; break; }
        if(dup) continue;
        int fid = comm.clusterFine[idx];
        if(fid < 0 || fid >= HCLUST_FINE_K) fid = 0;
        if(fineTake[fid] >= fineCap) continue;
        selected[selCount++] = idx;
        fineTake[fid]++;
      }
#else
      for(int i=0;i<topN;i++) selected[selCount++] = indices[i];
#endif
      for(int i=0;i<selCount;i++){
        int idx = selected[i];
        printf(" %d.%s: score=%.4f, C=%.4f\n", i+1, ASSET_NAMES[idx], (double)scores[idx], (double)compactness[idx]);
      }
    }
  }
};

// ---------------------------- Zorro DLL entry ----------------------------

static CompactDominantStrategy* S = NULL;

DLLFUNC void run()
{
  if(is(INITRUN)) {
    BarPeriod = 60;
    LookBack = max(LookBack, FEAT_WINDOW + 50);

    asset((char*)ASSET_NAMES[0]);

    if(!S) {
      S = new CompactDominantStrategy();
      S->init();
    }
  }

  if(is(EXITRUN)) {
    if(S) {
      S->shutdown();
      delete S;
      S = NULL;
    }
    return;
  }

  if(!S || Bar < LookBack)
    return;

  S->onBar();
}

Live in your Account [Re: TipmyPip] #489366
04/23/26 14:23
04/23/26 14:23
Joined: Sep 2017
Posts: 293
TipmyPip Offline OP
Member
TipmyPip  Offline OP
Member

Joined: Sep 2017
Posts: 293
Are we Progressing? Very Soon... We will all progress, and very very fast.

https://youtube.com/live/U8gBYk2XlFA?feature=share

Last edited by TipmyPip; 04/23/26 14:24.
Aether Volatility Compass [Re: TipmyPip] #489369
04/25/26 16:09
04/25/26 16:09
Joined: Sep 2017
Posts: 293
TipmyPip Offline OP
Member
TipmyPip  Offline OP
Member

Joined: Sep 2017
Posts: 293
The Mercury Loom is a trading strategy imagined as a careful artisan weaving three invisible threads through the market: price, volatility, and the hidden pull of interest rate behavior. Its purpose is not to chase every flicker of movement, but to listen for moments when the market’s surface becomes stretched away from its deeper rhythm. The code watches recent price motion, estimates how turbulent the market has become, and compares that living turbulence with a reference world shaped by local volatility and CEV style behavior. Around this, it builds a second layer inspired by Hull White rate dynamics, treating rates like a distant tide that can bend the path of price without always commanding it. When the volatility field appears overheated, the strategy favors reversal, as if expecting the market to return toward balance. When that signal is weak, it can fall back to pure momentum, following the current as a sailor might follow a clean wind. Its protective logic sizes stop loss and take profit from the estimated risk climate, so the trade breathes wider in stormy weather and tighter in calmer air. A drawdown guard acts like a gatekeeper, pausing the system when the journey becomes too hazardous. Overall, the code is symbolic of a navigator seeking order inside uncertainty, using calibration, restraint, and adaptation rather than blind prediction.

Code
// ============================================================================
// LocalOT_EA.c
// ============================================================================
// 1. Zorro calls run() once at initialization and then once at the end of every
//    completed bar. This script is therefore bar-driven, not tick-driven.
//
// 2. The script sets BarPeriod = 60 in INITRUN, so the intended working
//    timeframe is H1. The annualization constant ANNUAL_H1 is based on
//    252 trading days * 24 H1 bars per day.
//
// 3. Stop and TakeProfit are set as price distances before entering a trade.
//    This follows Zorro's normal trade setup order.
//
// 4. This script intentionally avoids the C ternary operator. lite-C supports
//    ifelse(), but this script uses plain if/else branches for clarity.
//
// 5. The mathematical routines are a practical implementation of the original
//    EA logic, not a full academic local-volatility/optimal-transport library.
//    They should be tested and validated on the target asset before any live
//    use.
//
// 6. Broker-specific lot constraints, margin rules, spread, commission, and pip
//    cost must be checked in the selected Zorro asset list. The variable
//    InpMaxLots is only a local safety cap.
//
// DATA FLOW SUMMARY
// -----------------
// Every bar after LookBack:
//
//     1. loadBars()
//        Loads recent closes, log returns, and realized volatility estimates.
//
//     2. calibrateAll() or lightUpdate()
//        Full recalibration is done every InpRetrainEvery bars.
//        Between full calibrations, only rate/proxy and local beta are updated.
//
//     3. signalLocalVolOT()
//        Builds a directional signal from:
//        - realized/reference volatility ratio,
//        - momentum,
//        - Hull-White short-rate proxy,
//        - lambda adjustment.
//
//     4. openLocalVolTrade()
//        Sets Stop, TakeProfit, Lots, then enters long or short.
//
//     5. shouldCloseLocalVolOT()
//        Exits when volatility has normalized.
//
//     6. ddExceeded()
//        Closes positions if the strategy drawdown exceeds InpMaxDD.
//
// RISK WARNING
// ------------
// This is research/educational code. It is not investment advice, and it should
// not be used live without independent review, backtesting, forward testing,
// broker-cost validation, and robust risk controls.
//
// ============================================================================


// ---------------------------------------------------------------------------
// Global constants
// ---------------------------------------------------------------------------

// Maximum number of observations stored in the fixed-size arrays below.
// lite-C supports normal arrays; fixed arrays are simple and predictable here.
// Keep this above InpCalibBars + InpVolWindow + 5.
#define MAX_OBS 600

// Number of H1 bars used for annualizing hourly variance.
// 252 trading days * 24 hourly bars. This matches the original H1 assumption.
#define ANNUAL_H1 (252.0*24.0)


// ---------------------------------------------------------------------------
// User parameters
// ---------------------------------------------------------------------------
// ----- CEV model parameters -----

// Reference annual volatility. Example: 0.10 means 10 percent annual vol.
var InpSigmaRef     = 0.10;

// CEV elasticity gamma. Values below 1 imply negative skew behavior.
var InpGammaRef     = 0.85;

// Cost exponent p used in the simplified optimal-transport beta formula.
var InpCostExp      = 4.0;


// ----- Hull-White short-rate model parameters -----

// Prior mean-reversion speed "a" for the short-rate proxy.
var InpHW_a         = 0.40;

// Prior short-rate volatility.
var InpHW_SigmaR    = 0.015;

// Prior initial short rate.
var InpHW_r0        = 0.025;

// Prior correlation between price/log-return shocks and short-rate shocks.
var InpCorrelRef    = -0.30;

// Bayesian prior weight used when blending input priors with empirical
// estimates. 0 means rely only on data; 1 means rely only on priors.
var InpPriorWeight  = 0.40;


// ----- Calibration parameters -----

// Number of bars used for calibration windows.
int InpCalibBars    = 500;

// Window length, in bars, for realized volatility estimation.
int InpVolWindow    = 20;

// Full recalibration interval, in bars.
int InpRetrainEvery = 48;

// Lambda update step. Lambda is a scalar correction term used in the
// potential/proxy calculation.
var InpLambdaStep   = 0.02;

// Retained for parity with the original EA. The current lite-C version does not
// use it as a stopping criterion.
var InpTol          = 1e-4;


// ----- Signal parameters -----

// In mean-reversion mode, trades are considered when realized/reference
// volatility is above this ratio.
var InpVolRatioHigh = 1.25;

// In trend-following mode, trades are considered when realized/reference
// volatility is below this ratio.
var InpVolRatioLow  = 0.75;

// Exit threshold. In mean-reversion mode, exit when ratio falls below this.
// In trend-following mode, exit when ratio rises above 1 / this value.
var InpVolRatioExit = 1.05;

// Number of bars used for the simple momentum component.
int InpMomBars      = 8;

// 1 = mean-reversion mode; 0 = trend-following mode.
int InpMeanRevMode  = 1;


// ----- Risk and execution parameters -----

// Percent of approximate strategy equity to risk per trade.
var InpRiskPct      = 1.0;

// Maximum drawdown, in percent, from this script's equity peak. If exceeded,
// all open positions are closed.
var InpMaxDD        = 10.0;

// Stop-loss multiplier. SL distance is based on local volatility.
var InpSL_k         = 2.5;

// Take-profit multiplier. TP is proportional to SL.
var InpTP_k         = 3.5;

// Local cap used by lotsFromSL(). This is not a broker query.
var InpMaxLots      = 100.0;


// ---------------------------------------------------------------------------
// Model state variables
// ---------------------------------------------------------------------------
// These variables hold the current calibrated model state. They are global
// because run() is called repeatedly by Zorro and local variables do not retain
// values between calls.


// Current calibrated CEV volatility.
var GM_SigmaCev;

// Current calibrated CEV elasticity.
var GM_GammaCev;

// Current calibrated Hull-White mean-reversion speed.
var GM_HWA;

// Current calibrated Hull-White short-rate volatility.
var GM_HWSigmaR;

// Current calibrated Hull-White drift/intercept proxy.
// In a standard form dr = (b - a*r)dt + sigma_r dW, this is the "b" term.
var GM_HWB0;

// Current calibrated return/rate correlation.
var GM_XiCorr;

// Current CEV reference variance at the latest price.
var GM_SigmaB2;

// Current local/OT-adjusted variance used for SL/TP sizing.
var GM_Beta11Star;

// Current lambda correction parameter.
var GM_Lambda;

// Current short-rate proxy.
var GM_RCurrent;

// Becomes 1 after the first successful full calibration.
int GM_Valid;

// Counts bars since the last full calibration.
int GM_BarsSinceCalib;


// ---------------------------------------------------------------------------
// Potential/proxy coefficients
// ---------------------------------------------------------------------------
// The original EA uses a local quadratic potential approximation. These
// coefficients store that approximation around the current log-price and rate.
// phi(z) is not directly evaluated in this script, but phi_z and phi_zz are
// used by lema412().

var GPhi_C0;
var GPhi_C1;
var GPhi_C2;
var GPhi_C3;
var GPhi_Z0;
var GPhi_R0;


// ---------------------------------------------------------------------------
// Rolling data arrays
// ---------------------------------------------------------------------------
// Index convention:
//     index 0 = current/latest bar
//     index 1 = previous bar
//     index i = i bars ago
//
// GPrice stores prices.
// GLRet stores log returns log(price[i] / price[i+1]).
// GRVol stores realized annualized volatility estimates.
// GRate stores the short-rate proxy path.

var GPrice[MAX_OBS];
var GLRet[MAX_OBS];
var GRVol[MAX_OBS];
var GRate[MAX_OBS];

// Number of valid price observations in GPrice.
int GNp;

// Number of valid realized-vol observations in GRVol.
int GNr;

// Becomes 1 after the first calibration and banner print.
int GStarted;

// Strategy-equity peak used for drawdown control.
var GEqPeak;


// ===========================================================================
// Utility functions
// ===========================================================================


// ---------------------------------------------------------------------------
// phiDz
// ---------------------------------------------------------------------------
// Returns first derivative of the local quadratic potential approximation
// with respect to log price Z.
//
// With the current approximation:
//     phi_z(Z) = C1 + 2*C2*(Z - Z0)
//
// Parameters:
//     Z - log price at which the derivative is evaluated.
//
// Returns:
//     First derivative phi_z.
// ---------------------------------------------------------------------------
var phiDz(var Z)
{
    return GPhi_C1 + 2.0*GPhi_C2*(Z - GPhi_Z0);
}


// ---------------------------------------------------------------------------
// phiDzz
// ---------------------------------------------------------------------------
// Returns second derivative of the local quadratic potential approximation.
//
// With the current approximation:
//     phi_zz = 2*C2
//
// Returns:
//     Second derivative phi_zz.
// ---------------------------------------------------------------------------
var phiDzz()
{
    return 2.0*GPhi_C2;
}


// ---------------------------------------------------------------------------
// strategyEquity
// ---------------------------------------------------------------------------
// Calculates an approximate strategy equity value from Zorro's global capital
// and profit variables.
//
// This is used for risk sizing and drawdown checks. For portfolio-level
// scripts, confirm whether Capital + ProfitClosed + ProfitOpen is the desired
// equity base.
//
// Returns:
//     Approximate current strategy equity.
// ---------------------------------------------------------------------------
var strategyEquity()
{
    // Zorro-native approximation of the strategy equity.
    // Capital is the initial capital; ProfitClosed and ProfitOpen are strategy P/L.
    var E = Capital + ProfitClosed + ProfitOpen;
    if(E <= 0.0)
        E = Capital;
    return E;
}


// ---------------------------------------------------------------------------
// initModel
// ---------------------------------------------------------------------------
// Initializes all model-state variables to their prior/default values.
//
// Called from run() during INITRUN. This is important because static/global
// variables can persist between script restarts when auto-compilation behavior
// does not force a fresh process.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void initModel()
{
    GM_SigmaCev       = InpSigmaRef;
    GM_GammaCev       = InpGammaRef;
    GM_HWA            = InpHW_a;
    GM_HWSigmaR       = InpHW_SigmaR;
    GM_HWB0           = InpHW_r0 * InpHW_a;
    GM_XiCorr         = InpCorrelRef;
    GM_SigmaB2        = InpSigmaRef * InpSigmaRef;
    GM_Beta11Star     = InpSigmaRef * InpSigmaRef;
    GM_Lambda         = 0.0;
    GM_RCurrent       = InpHW_r0;
    GM_Valid          = 0;
    GM_BarsSinceCalib = 0;

    GPhi_C0 = 0.0;
    GPhi_C1 = 0.0;
    GPhi_C2 = 0.0;
    GPhi_C3 = 0.0;
    GPhi_Z0 = 0.0;
    GPhi_R0 = 0.0;

    GNp = 0;
    GNr = 0;
    GStarted = 0;
    GEqPeak = Capital;
}


// ---------------------------------------------------------------------------
// loadBars
// ---------------------------------------------------------------------------
// Loads recent prices and calculates log returns and realized volatilities.
//
// The function uses Zorro priceClose(i), where i is the bar offset:
//     priceClose(0) = current/latest completed bar close
//     priceClose(1) = previous bar close
//
// It fills:
//     GPrice[0..GNp-1]
//     GLRet[0..GNp-2]
//     GRVol[0..GNr-1]
//
// Realized volatility is calculated from InpVolWindow log returns and then
// annualized by ANNUAL_H1.
//
// Returns:
//     1 when enough valid data was loaded.
//     0 when data is insufficient or invalid.
// ---------------------------------------------------------------------------
int loadBars()
{
    int Need = InpCalibBars + InpVolWindow + 5;
    if(Need > MAX_OBS)
        Need = MAX_OBS;
    if(Need < InpVolWindow + 25)
        return 0;
    if(Bar < LookBack)
        return 0;

    GNp = Need;

    int i;
    for(i = 0; i < GNp; i++)
    {
        var C = priceClose(i);
        if(C <= 0.0)
            return 0;
        GPrice[i] = C;
    }

    // Log returns in most-recent-first order.
    for(i = 0; i < GNp-1; i++)
    {
        if(GPrice[i+1] <= 0.0)
            GLRet[i] = 0.0;
        else
            GLRet[i] = log(GPrice[i] / GPrice[i+1]);
    }

    GNr = GNp - InpVolWindow;
    if(GNr < 1)
        return 0;

    // Rolling sample variance of log returns, then annualized volatility.
    for(i = 0; i < GNr; i++)
    {
        var S = 0.0;
        var S2 = 0.0;
        int j;
        for(j = i; j < i + InpVolWindow; j++)
        {
            S += GLRet[j];
            S2 += GLRet[j]*GLRet[j];
        }
        var VarH1 = (S2 - S*S/InpVolWindow) / (InpVolWindow - 1);
        if(VarH1 < 0.0)
            VarH1 = 0.0;
        GRVol[i] = sqrt(VarH1 * ANNUAL_H1);
    }

    return 1;
}


// ---------------------------------------------------------------------------
// estimateRate
// ---------------------------------------------------------------------------
// Estimates the current short-rate proxy from realized price drift.
//
// Original rationale:
//     E[dlog(S)] approx r - 0.5*sigma^2
// Therefore:
//     r approx drift + 0.5*sigma^2
//
// Since spot price data does not directly contain a short-rate curve, this is
// only a practical proxy. The estimate is blended with InpHW_r0 using
// InpPriorWeight.
//
// The function also generates a smooth Hull-White-style rate path in GRate.
// This path is used by calibHW() and calibCorr().
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void estimateRate()
{
    int W = InpCalibBars;
    if(W > GNp - 1)
        W = GNp - 1;
    if(W < 1)
        return;

    var SumRet = 0.0;
    int i;
    for(i = 0; i < W; i++)
        SumRet += GLRet[i];

    var DriftAnn = (SumRet / W) * ANNUAL_H1;

    var RVol0 = InpSigmaRef;
    if(GNr > 0)
        RVol0 = GRVol[0];

    var REst = DriftAnn + 0.5*RVol0*RVol0;
    REst = clamp(REst,-0.05,0.20);

    // Bayesian-style blend of empirical rate proxy and prior rate.
    GM_RCurrent = (1.0 - InpPriorWeight) * REst + InpPriorWeight * InpHW_r0;

    // Build a decaying proxy path around the current rate.
    var A = GM_HWA;
    var BA = GM_HWB0 / max(A,0.001);
    for(i = 0; i < GNp; i++)
    {
        var T = i / ANNUAL_H1;
        var Decay = exp(-A*T);
        GRate[i] = GM_RCurrent*Decay + BA*(1.0 - Decay);
    }
    GRate[0] = GM_RCurrent;
}


// ===========================================================================
// Calibration blocks
// ===========================================================================


// ---------------------------------------------------------------------------
// calibCEV
// ---------------------------------------------------------------------------
// Calibrates a simplified CEV relationship from realized volatility and price.
//
// The CEV variance scaling is represented by:
//
//     sigma(S) = sigma_ref * (S/S0)^(gamma - 1)
//
// Taking logs gives:
//
//     log(sigma) = log(sigma_ref) + (gamma - 1)*log(S/S0)
//
// This function runs a simple linear regression of log(realized vol) on
// log(price). The slope implies gamma - 1. The intercept implies sigma.
//
// Empirical estimates are clamped and then blended with the user priors.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void calibCEV()
{
    int N = GNr - 1;
    if(N > InpCalibBars)
        N = InpCalibBars;
    if(N < 20)
        return;

    var SX = 0.0;
    var SY = 0.0;
    var SXX = 0.0;
    var SXY = 0.0;
    int Cnt = 0;

    int i;
    for(i = 0; i < N; i++)
    {
        if(GRVol[i] < 1e-8)
            continue;
        if(GPrice[i] <= 0.0)
            continue;

        var Y = log(GRVol[i]);
        var X = log(GPrice[i]);
        SX += X;
        SY += Y;
        SXX += X*X;
        SXY += X*Y;
        Cnt += 1;
    }

    if(Cnt < 15)
        return;

    var D = Cnt*SXX - SX*SX;
    if(abs(D) < 1e-12)
        return;

    var Slope = (Cnt*SXY - SX*SY) / D;
    var Intercept = (SY - Slope*SX) / Cnt;

    var GammaNew = Slope + 1.0;
    var SigmaNew = exp(Intercept);

    // Clamp empirical estimates to prevent unstable model parameters.
    GammaNew = clamp(GammaNew,0.30,1.70);
    SigmaNew = clamp(SigmaNew,0.01,3.0);

    // Blend empirical estimate with prior.
    var W = InpPriorWeight;
    GM_GammaCev = W*InpGammaRef + (1.0-W)*GammaNew;
    GM_SigmaCev = W*InpSigmaRef + (1.0-W)*SigmaNew;

    printf("\n[CEV] sigma=%.4f gamma=%.4f n=%i prior_w=%.2f",
        GM_SigmaCev, GM_GammaCev, Cnt, W);
}


// ---------------------------------------------------------------------------
// calibHW
// ---------------------------------------------------------------------------
// Calibrates a simplified Hull-White short-rate proxy.
//
// Model form used here:
//
//     dr = (b - a*r) dt + sigma_r dW
//
// The function estimates alpha and beta from:
//
//     dr approx alpha + beta*r
//
// Then:
//     a approx -beta / dt
//     b approx alpha / dt
//
// sigma_r is estimated from residual variance. All estimates are clamped and
// blended with priors to avoid degenerate values.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void calibHW()
{
    int N = GNp - 2;
    if(N > InpCalibBars)
        N = InpCalibBars;
    if(N < 20)
        return;

    var DT = 1.0 / ANNUAL_H1;
    var SX = 0.0;
    var SY = 0.0;
    var SXX = 0.0;
    var SXY = 0.0;
    var SRes = 0.0;
    int Cnt = 0;

    int i;
    for(i = 0; i < N-1; i++)
    {
        var R0 = GRate[i+1];
        var DR = GRate[i] - R0;
        SX += R0;
        SY += DR;
        SXX += R0*R0;
        SXY += R0*DR;
        Cnt += 1;
    }

    if(Cnt < 10)
        return;

    var D = Cnt*SXX - SX*SX;
    var AOls = 0.0;
    var BOls = 0.0;
    var SROls = 0.0;

    if(abs(D) > 1e-14)
    {
        var Beta = (Cnt*SXY - SX*SY) / D;
        var Alpha = (SY - Beta*SX) / Cnt;

        AOls = -Beta / DT;
        BOls = Alpha / DT;

        for(i = 0; i < Cnt; i++)
        {
            var R0b = GRate[i+1];
            var DRb = GRate[i] - R0b;
            var Res = DRb - (Alpha + Beta*R0b);
            SRes += Res*Res;
        }

        SROls = sqrt(SRes / Cnt / DT);
    }

    // Prior-blended and clamped estimates.
    var W = InpPriorWeight;
    GM_HWA      = max(0.05, W*InpHW_a + (1.0-W)*clamp(AOls,0.05,5.0));
    GM_HWB0     = W*(InpHW_r0*InpHW_a) + (1.0-W)*BOls;
    GM_HWSigmaR = max(0.005, W*InpHW_SigmaR + (1.0-W)*clamp(SROls,0.005,0.50));

    printf("\n[HW] a=%.4f sigmar=%.4f b0=%.4f",
        GM_HWA, GM_HWSigmaR, GM_HWB0);
}


// ---------------------------------------------------------------------------
// calibCorr
// ---------------------------------------------------------------------------
// Estimates the empirical correlation between:
//     - changes in the short-rate proxy, and
//     - log returns.
//
// The estimate is blended with InpCorrelRef. If the empirical correlation is
// very close to zero, the prior weight is increased. This prevents a flat or
// noisy rate proxy from completely removing the model's directional component.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void calibCorr()
{
    int N = GNp;
    if(GNr < N)
        N = GNr;
    N -= 2;
    if(N > InpCalibBars)
        N = InpCalibBars;
    if(N < 20)
        return;

    var SX = 0.0;
    var SY = 0.0;
    var SXX = 0.0;
    var SYY = 0.0;
    var SXY = 0.0;
    int Cnt = 0;

    int i;
    for(i = 0; i < N-1; i++)
    {
        var DR = GRate[i] - GRate[i+1];
        var DZ = GLRet[i];
        SX += DR;
        SY += DZ;
        SXX += DR*DR;
        SYY += DZ*DZ;
        SXY += DR*DZ;
        Cnt += 1;
    }

    if(Cnt < 10)
        return;

    var Cov = (SXY - SX*SY/Cnt) / Cnt;
    var VX = (SXX - SX*SX/Cnt) / Cnt;
    var VY = (SYY - SY*SY/Cnt) / Cnt;

    if(VX < 1e-14)
        return;
    if(VY < 1e-14)
        return;

    var XiOls = Cov / sqrt(VX*VY);
    XiOls = clamp(XiOls,-0.999,0.999);

    var W = InpPriorWeight;

    // If data suggests near-zero correlation, avoid over-trusting the noisy
    // estimate. This implements the source EA's fix for flat-rate cases.
    if(abs(XiOls) < 0.05)
        W = 0.80;

    GM_XiCorr = W*InpCorrelRef + (1.0-W)*XiOls;

    printf("\n[CORR] xi=%.4f OLS=%.4f prior_w=%.2f",
        GM_XiCorr, XiOls, W);
}


// ===========================================================================
// Mathematical core
// ===========================================================================


// ---------------------------------------------------------------------------
// refVar
// ---------------------------------------------------------------------------
// Calculates the CEV reference variance at a given log price.
//
// The CEV variance proxy is:
//
//     sigma_b^2(S) = sigma_cev^2 * (S/S0)^(2*(gamma_cev - 1))
//
// Parameters:
//     LogPrice - log(S), where S is the price.
//
// Returns:
//     Reference variance. A small positive floor is applied.
// ---------------------------------------------------------------------------
var refVar(var LogPrice)
{
    var S = exp(LogPrice);
    var S0 = GPrice[GNp - 1];
    if(S0 <= 0.0)
        S0 = GPrice[0];

    var Expon = 2.0*(GM_GammaCev - 1.0);
    var Ratio = S / S0;
    if(Ratio < 0.0001)
        Ratio = 0.0001;

    return max(1e-10, GM_SigmaCev*GM_SigmaCev*pow(Ratio,Expon));
}


// ---------------------------------------------------------------------------
// buildPhi
// ---------------------------------------------------------------------------
// Builds a local quadratic potential approximation around current log price
// and current short-rate proxy.
//
// The approximation is intentionally simple. It uses:
//     - CEV reference variance,
//     - current realized variance,
//     - current rate proxy.
//
// Parameters:
//     Z0 - current log price.
//     R0 - current short-rate proxy.
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void buildPhi(var Z0, var R0)
{
    GPhi_Z0 = Z0;
    GPhi_R0 = R0;

    var SB2 = refVar(Z0);
    var RV2 = SB2;
    if(GNr > 0)
        RV2 = GRVol[0]*GRVol[0];

    // Mix of reference variance and realized variance. The floor prevents
    // division by zero or excessive curvature.
    var Mix = 0.5*(SB2 + RV2);
    Mix = max(Mix,1e-8);

    GPhi_C2 = 1.0 / Mix;
    GPhi_C1 = (R0 - 0.5*Mix) / Mix;
    GPhi_C0 = GM_Lambda;
    GPhi_C3 = GPhi_C1 / max(GM_HWA,0.05);
}


// ---------------------------------------------------------------------------
// lema412
// ---------------------------------------------------------------------------
// Practical implementation of the EA's Lemma 4.12-inspired beta calculation.
//
// Inputs are simplified/local quantities:
//     PhiZZ   - second derivative of potential.
//     PhiZ    - first derivative of potential.
//     SigmaB2 - CEV reference variance.
//     XiRef   - return/rate correlation.
//     Sigr2   - short-rate variance.
//     P       - cost exponent.
//
// The function returns beta11*, an adjusted local variance used for stop-loss
// and take-profit sizing.
//
// Stability guards:
//     - Diff floor prevents invalid powers.
//     - D clamp prevents explosive amplification when SigmaB2 is small.
//     - Final beta is clamped to [0.5*SigmaB2, 2.0*SigmaB2].
//
// Returns:
//     Adjusted local variance beta11*.
// ---------------------------------------------------------------------------
var lema412(var PhiZZ, var PhiZ, var SigmaB2, var XiRef, var Sigr2, var P)
{
    var S = XiRef*XiRef*Sigr2;
    var Diff = SigmaB2 - S;

    if(Diff < 1e-12)
        return max(SigmaB2,S + 1e-12);

    var DMax = 2.0 / SigmaB2;
    var D = PhiZZ - PhiZ;
    D = clamp(D,-DMax,DMax);

    var DiffP = pow(Diff,P);
    var A = DiffP * D / (2.0*(P*P - 1.0));

    var Disc = A*A + 4.0*pow(Diff,2.0*P);
    if(Disc < 0.0)
        Disc = 0.0;

    var U = (A + sqrt(Disc))*0.5;
    if(U <= 0.0)
        return SigmaB2;

    var Beta11 = S + pow(U,1.0/P);

    return clamp(Beta11,0.50*SigmaB2,2.0*SigmaB2);
}


// ---------------------------------------------------------------------------
// calibrateAll
// ---------------------------------------------------------------------------
// Runs a full model calibration cycle:
//
//     1. CEV calibration.
//     2. Hull-White proxy calibration.
//     3. Correlation calibration.
//     4. Rate proxy update.
//     5. Potential/proxy build.
//     6. beta11* calculation.
//     7. lambda update.
//
// This function sets GM_Valid = 1 when complete.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void calibrateAll()
{
    printf("\n[LocalVolOT] Full calibration");

    calibCEV();
    calibHW();
    calibCorr();

    estimateRate();

    var Z0 = log(GPrice[0]);
    var R0 = GM_RCurrent;

    buildPhi(Z0,R0);

    var SB2 = refVar(Z0);
    var Sigr2 = GM_HWSigmaR*GM_HWSigmaR;

    GM_SigmaB2 = SB2;
    GM_Beta11Star = lema412(phiDzz(),phiDz(Z0),SB2,GM_XiCorr,Sigr2,InpCostExp);

    var RV2 = SB2;
    if(GNr > 0)
        RV2 = GRVol[0]*GRVol[0];

    // Lambda is adjusted toward the difference between realized variance and
    // the CEV reference variance.
    var Grad = RV2 - GM_SigmaB2;
    GM_Lambda += InpLambdaStep * Grad;
    GM_Lambda = clamp(GM_Lambda,-3.0,3.0);

    GM_Valid = 1;
    GM_BarsSinceCalib = 0;

    var VolRatio = 1.0;
    if(GM_SigmaB2 > 1e-10)
        VolRatio = sqrt(RV2) / sqrt(GM_SigmaB2);

    printf("\n[OT] sig_ref=%.2f%% sig_star=%.2f%% rv=%.2f%% ratio=%.3f lambda=%.4f r=%.3f%% xi=%.4f",
        sqrt(GM_SigmaB2)*100.0,
        sqrt(GM_Beta11Star)*100.0,
        sqrt(RV2)*100.0,
        VolRatio,
        GM_Lambda,
        GM_RCurrent*100.0,
        GM_XiCorr);
}


// ---------------------------------------------------------------------------
// lightUpdate
// ---------------------------------------------------------------------------
// Performs a lighter per-bar model update between full calibrations.
//
// It updates:
//     - short-rate proxy,
//     - potential/proxy coefficients,
//     - CEV reference variance,
//     - beta11*.
//
// It does not rerun the CEV/HW/correlation regressions.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void lightUpdate()
{
    estimateRate();

    var Z0 = log(GPrice[0]);
    buildPhi(Z0,GM_RCurrent);

    var SB2 = refVar(Z0);
    var Sigr2 = GM_HWSigmaR*GM_HWSigmaR;

    GM_SigmaB2 = SB2;
    GM_Beta11Star = lema412(phiDzz(),phiDz(Z0),SB2,GM_XiCorr,Sigr2,InpCostExp);
}


// ===========================================================================
// Signal functions
// ===========================================================================


// ---------------------------------------------------------------------------
// signalLocalVolOT
// ---------------------------------------------------------------------------
// Generates the trade direction signal.
//
// Returns:
//     +1 = long signal.
//     -1 = short signal.
//      0 = no signal.
//
// Main components:
//     VolRatio:
//         realized volatility / CEV reference volatility.
//
//     Mom:
//         average recent log return over InpMomBars.
//
//     RateDir:
//         short-rate drift proxy scaled by rate volatility and correlation.
//
//     LambdaDir:
//         direction implied by the lambda correction.
//
//     Dir:
//         weighted sum of momentum, rate direction, and lambda direction.
//
// Signal interpretation:
//     Mean-reversion mode:
//         Only act when VolRatio is high. Fade the directional component.
//
//     Trend-following mode:
//         Only act when VolRatio is low. Follow the directional component.
// ---------------------------------------------------------------------------
int signalLocalVolOT()
{
    if(!GM_Valid)
        return 0;

    var RV = GM_SigmaCev;
    if(GNr > 0)
        RV = GRVol[0];

    var SB = sqrt(GM_SigmaB2);
    if(SB < 1e-8)
        return 0;

    var VolRatio = RV / SB;

    int MB = InpMomBars;
    if(MB < 2)
        MB = 2;
    if(MB > GNp - 1)
        MB = GNp - 1;

    var Mom = 0.0;
    int i;
    for(i = 0; i < MB; i++)
        Mom += GLRet[i];
    Mom /= MB;

    var DT = 1.0 / ANNUAL_H1;
    var R0 = GM_RCurrent;
    var DR = (GM_HWB0 - GM_HWA*R0) * DT;

    var RateZ = 0.0;
    if(GM_HWSigmaR > 1e-6)
        RateZ = DR / (GM_HWSigmaR * sqrt(DT));

    var RateDir = GM_XiCorr * RateZ;

    var LambdaDir = -GM_Lambda * 0.5;

    var Dir = 0.60*Mom + 0.25*RateDir*abs(GM_XiCorr) + 0.15*LambdaDir;

    int Sig = 0;

    if(InpMeanRevMode)
    {
        // High realized/reference volatility: fade the direction.
        if(VolRatio >= InpVolRatioHigh)
        {
            if(Dir < -1e-5)
                Sig = 1;
            else if(Dir > 1e-5)
                Sig = -1;
            else if(RateDir > 0.01 && GM_XiCorr > 0.05)
                Sig = -1;
            else if(RateDir < -0.01 && GM_XiCorr < -0.05)
                Sig = 1;
        }
    }
    else
    {
        // Low realized/reference volatility: follow the direction.
        if(VolRatio <= InpVolRatioLow)
        {
            if(Dir < -1e-5)
                Sig = -1;
            else if(Dir > 1e-5)
                Sig = 1;
        }
    }

    return Sig;
}


// ---------------------------------------------------------------------------
// shouldCloseLocalVolOT
// ---------------------------------------------------------------------------
// Decides whether existing positions should be closed because volatility has
// normalized.
//
// Returns:
//     1 = close open positions.
//     0 = keep positions open.
//
// Logic:
//     Mean-reversion mode:
//         close when VolRatio < InpVolRatioExit.
//
//     Trend-following mode:
//         close when VolRatio > 1 / InpVolRatioExit.
// ---------------------------------------------------------------------------
int shouldCloseLocalVolOT()
{
    if(!GM_Valid)
        return 0;

    var RV = GM_SigmaCev;
    if(GNr > 0)
        RV = GRVol[0];

    var SB = sqrt(GM_SigmaB2);
    if(SB < 1e-8)
        return 0;

    var VolRatio = RV / SB;

    if(InpMeanRevMode)
    {
        if(VolRatio < InpVolRatioExit)
            return 1;
    }
    else
    {
        if(VolRatio > 1.0/InpVolRatioExit)
            return 1;
    }

    return 0;
}


// ===========================================================================
// Execution and risk functions
// ===========================================================================


// ---------------------------------------------------------------------------
// calcStops
// ---------------------------------------------------------------------------
// Calculates stop-loss and take-profit distances.
//
// The stop distance combines:
//     1. Price move implied by local beta11* volatility.
//     2. A Hull-White/correlation contribution.
//
// The final Stop and TakeProfit values are distances in price units, not
// absolute price levels. They are later assigned to Zorro's Stop and
// TakeProfit variables before entering a trade.
//
// Parameters:
//     Dir    - trade direction. Currently not used directly, retained for
//              source parity and future directional adjustments.
//     SLDist - output pointer for stop-loss distance.
//     TPDist - output pointer for take-profit distance.
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void calcStops(int Dir,var* SLDist,var* TPDist)
{
    var S = GPrice[0];
    var SigAnn = sqrt(GM_Beta11Star);
    var SigH1 = SigAnn / sqrt(ANNUAL_H1);
    var DistS = InpSL_k * SigH1 * S;

    var HWContrib = GM_HWSigmaR * S * abs(GM_XiCorr) / sqrt(ANNUAL_H1);
    var SL = sqrt(DistS*DistS + HWContrib*HWContrib);
    var TP = SL * (InpTP_k / InpSL_k);

    // MT5 source used 50 points, usually about 5 pips on 5-digit FX.
    // Zorro's PIP variable converts this to the current asset's price units.
    var MinD = 5.0 * PIP;
    SL = max(SL,MinD);
    TP = max(TP,MinD);

    *SLDist = SL;
    *TPDist = TP;
}


// ---------------------------------------------------------------------------
// lotsFromSL
// ---------------------------------------------------------------------------
// Converts stop-loss distance into a Zorro Lots value from risk percent.
//
// Approximate formula:
//
//     risk cash    = strategy equity * InpRiskPct / 100
//     risk per lot = (SL distance / PIP) * PIPCost
//     lots         = floor(risk cash / risk per lot)
//
// If PIP or PIPCost are unavailable, the function falls back to 1 lot.
//
// Parameters:
//     SLDist - stop-loss distance in price units.
//
// Returns:
//     Number of lots to trade, clamped to [1, InpMaxLots] when possible.
// ---------------------------------------------------------------------------
var lotsFromSL(var SLDist)
{
    if(SLDist <= 0.0)
        return 0.0;
    if(PIP <= 0.0)
        return 1.0;
    if(PIPCost <= 0.0)
        return 1.0;

    var RiskCash = strategyEquity() * InpRiskPct / 100.0;
    var RiskPerLot = (SLDist / PIP) * PIPCost;

    if(RiskPerLot <= 0.0)
        return 1.0;

    var L = floor(RiskCash / RiskPerLot);
    if(L < 1.0)
        L = 1.0;
    if(InpMaxLots > 0.0)
        L = min(L,InpMaxLots);

    return L;
}


// ---------------------------------------------------------------------------
// openLocalVolTrade
// ---------------------------------------------------------------------------
// Opens a long or short trade from the model signal.
//
// Flow:
//     1. Calculate SL/TP distances.
//     2. Assign Zorro Stop and TakeProfit before entry.
//     3. Calculate Lots.
//     4. Close opposite position if present.
//     5. Enter new position if no same-direction position is open.
//
// Parameters:
//     Dir - +1 for long, -1 for short.
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void openLocalVolTrade(int Dir)
{
    var SLDist = 0.0;
    var TPDist = 0.0;

    calcStops(Dir,&SLDist,&TPDist);

    // Zorro reads Stop, TakeProfit, and Lots when enterLong/enterShort is
    // called, so they must be set before entry.
    Stop = SLDist;
    TakeProfit = TPDist;
    Lots = lotsFromSL(SLDist);

    if(Lots <= 0.0)
        return;

    var RV = GM_SigmaCev;
    if(GNr > 0)
        RV = GRVol[0];

    var SB = sqrt(GM_SigmaB2);
    var Ratio = 0.0;
    if(SB > 0.0)
        Ratio = RV / SB;

    if(Dir > 0)
    {
        if(NumOpenShort > 0)
            exitShort();
        if(NumOpenLong == 0)
        {
            enterLong();
            printf("\n[TRADE] BUY lots=%.2f stop=%.5f tp=%.5f sig_star=%.2f%% ratio=%.3f",
                Lots, SLDist, TPDist, sqrt(GM_Beta11Star)*100.0, Ratio);
        }
    }
    else if(Dir < 0)
    {
        if(NumOpenLong > 0)
            exitLong();
        if(NumOpenShort == 0)
        {
            enterShort();
            printf("\n[TRADE] SELL lots=%.2f stop=%.5f tp=%.5f sig_star=%.2f%% ratio=%.3f",
                Lots, SLDist, TPDist, sqrt(GM_Beta11Star)*100.0, Ratio);
        }
    }
}


// ---------------------------------------------------------------------------
// closeLocalVolPositions
// ---------------------------------------------------------------------------
// Closes all positions opened in the current long/short direction context.
//
// Parameters:
//     Why - text printed to the Zorro log.
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void closeLocalVolPositions(string Why)
{
    if(NumOpenLong > 0)
        exitLong();
    if(NumOpenShort > 0)
        exitShort();

    printf("\n[CLOSE] %s",Why);
}


// ---------------------------------------------------------------------------
// ddExceeded
// ---------------------------------------------------------------------------
// Checks whether the strategy drawdown exceeds InpMaxDD.
//
// The drawdown base is the peak of strategyEquity() observed since the script
// started. If the current equity falls below that peak by more than InpMaxDD
// percent, the function returns 1.
//
// Returns:
//     1 = maximum drawdown exceeded.
//     0 = drawdown is within limits.
// ---------------------------------------------------------------------------
int ddExceeded()
{
    var E = strategyEquity();

    if(E > GEqPeak)
        GEqPeak = E;

    if(GEqPeak <= 0.0)
        return 0;

    var DD = (GEqPeak - E) / GEqPeak * 100.0;
    if(DD > InpMaxDD)
        return 1;

    return 0;
}


// ---------------------------------------------------------------------------
// printBanner
// ---------------------------------------------------------------------------
// Prints the initial calibrated model state to the Zorro log.
//
// Parameters:
//     none
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
void printBanner()
{
    var RV = GM_SigmaCev;
    if(GNr > 0)
        RV = GRVol[0];

    printf("\n==================================================");
    printf("\nLocalVolOT v2 - Zorro lite-C conversion");
    printf("\nCEV: sigma=%.3f gamma=%.3f",GM_SigmaCev,GM_GammaCev);
    printf("\nHW: a=%.3f sigmar=%.4f r0=%.3f%%",GM_HWA,GM_HWSigmaR,GM_RCurrent*100.0);
    printf("\nCorr: xi=%.4f prior_w=%.2f",GM_XiCorr,InpPriorWeight);
    printf("\nVols: sig_ref=%.2f%% sig_star=%.2f%% rv=%.2f%%",
        sqrt(GM_SigmaB2)*100.0, sqrt(GM_Beta11Star)*100.0, RV*100.0);
    if(InpMeanRevMode)
        printf("\nMode: MEAN-REVERSION risk=%.1f%%",InpRiskPct);
    else
        printf("\nMode: TREND-FOLLOW risk=%.1f%%",InpRiskPct);
    printf("\n==================================================");
}


// ===========================================================================
// Zorro strategy entry point
// ===========================================================================


// ---------------------------------------------------------------------------
// run
// ---------------------------------------------------------------------------
// Main Zorro strategy function.
//
// Zorro execution behavior:
//     - On the first initialization pass, is(INITRUN) is true.
//     - During the LookBack period, is(LOOKBACK) is true.
//     - After LookBack, this function runs once per completed bar.
//     - On script exit, is(EXITRUN) is true.
//
// Main flow after LookBack:
//     1. Load recent bars.
//     2. Perform initial/full/light calibration.
//     3. Check drawdown limit.
//     4. Print model state.
//     5. Close existing trades if vol normalized.
//     6. Otherwise evaluate a new signal and enter if applicable.
//
// Returns:
//     nothing
// ---------------------------------------------------------------------------
function run()
{
    if(is(INITRUN))
    {
        // H1 timeframe, matching PERIOD_H1 in the MT5 source.
        BarPeriod = 60;

        // Need enough bars for calibration, realized-vol window, and a margin.
        LookBack = InpCalibBars + InpVolWindow + 10;

        // Default starting capital. Adjust to your test account assumptions.
        Capital = 10000;

        // This strategy is designed to hold at most one long or one short.
        MaxLong = 1;
        MaxShort = 1;

        // Use available local history for LookBack in live mode where possible.
        set(PRELOAD);

        initModel();
        return;
    }

    if(is(EXITRUN))
    {
        printf("\n[LocalVolOT v2] Exit. lambda=%.4f beta11=%.6f xi=%.4f",
            GM_Lambda, GM_Beta11Star, GM_XiCorr);
        return;
    }

    // Do not trade during LookBack. Indicators/model arrays are not yet ready.
    if(is(LOOKBACK))
        return;

    if(!loadBars())
    {
        printf("\n[LocalVolOT] Not enough data yet.");
        return;
    }

    // First post-LookBack bar: initialize equity peak and perform full model
    // calibration. No trade is entered on this bar.
    if(!GStarted)
    {
        GEqPeak = strategyEquity();
        calibrateAll();
        printBanner();
        GStarted = 1;
        return;
    }

    GM_BarsSinceCalib += 1;

    // Hard strategy-level drawdown control.
    if(ddExceeded())
    {
        closeLocalVolPositions("Maximum drawdown exceeded");
        return;
    }

    // Full recalibration schedule.
    if(!GM_Valid)
        calibrateAll();
    else if(GM_BarsSinceCalib >= InpRetrainEvery)
        calibrateAll();
    else
        lightUpdate();

    if(!GM_Valid)
        return;

    // Log the current model state every bar for easier debugging.
    var RV = GM_SigmaCev;
    if(GNr > 0)
        RV = GRVol[0];

    var SB = sqrt(GM_SigmaB2);
    var Ratio = 0.0;
    if(SB > 0.0)
        Ratio = RV / SB;

    printf("\n[BAR %i] rv=%.2f%% sig_ref=%.2f%% ratio=%.3f sig_star=%.2f%% lambda=%.4f r=%.3f%%",
        Bar,
        RV*100.0,
        SB*100.0,
        Ratio,
        sqrt(GM_Beta11Star)*100.0,
        GM_Lambda,
        GM_RCurrent*100.0);

    // Existing position management. This script exits positions when the
    // volatility ratio normalizes; Stop and TakeProfit can also close trades.
    if(NumOpenLong + NumOpenShort > 0)
    {
        if(shouldCloseLocalVolOT())
            closeLocalVolPositions("Vol normalized");
        return;
    }

    // No open position: evaluate and trade a fresh signal.
    int Sig = signalLocalVolOT();
    if(Sig != 0)
        openLocalVolTrade(Sig);
}

Last edited by TipmyPip; 04/25/26 16:12.
Page 23 of 23 1 2 21 22 23

Moderated by  Petra 

Powered by UBB.threads™ PHP Forum Software 7.7.1