feat: add and fix cuda files

Author: dancixx

Cargo.toml

@@ -20,15 +20,16 @@ candle-datasets = "0.8.1"
 candle-nn = "0.8.1"
 candle-transformers = "0.8.1"
 chrono = "0.4.38"
-cudarc = { version = "0.13.9", optional = true, features = ["cuda-12080"] }
+cudarc = { version = "0.13.9", optional = true, features = ["cuda-12080", "cuda-version-from-build-system"] }
 flate2 = "1.0.34"
 gauss-quad = "0.2.1"
 impl-new-derive = "0.1.2"
 implied-vol = "1.0.0"
 indicatif = "0.17.8"
-itransformer = "1.0.1"
+# itransformer = "1.0.1"
 kendalls = "0.2.2"
 levenberg-marquardt = "0.14.0"
+libloading = { version = "0.8.6", optional = true }
 linreg = "0.2.0"
 mimalloc = { version = "0.1.43", optional = true }
 nalgebra = "0.33.2"
@@ -37,7 +38,7 @@ ndarray = { version = "0.16.1", features = [
     "matrixmultiply-threading",
     "blas",
 ] }
-ndarray-linalg = { version = "0.17.0", features = ["openblas-static"] }
+ndarray-linalg = { version = "0.17.0"}
 ndarray-npy = "0.9.1"
 ndarray-rand = "0.15.0"
 ndarray-stats = "0.6.0"
@@ -68,8 +69,8 @@ yahoo_finance_api = { version = "2.3.0", optional = true }
 [dev-dependencies]
 
 [features]
-cuda = ["dep:cudarc"]
-default = ["jemalloc"]
+cuda = ["dep:cudarc", "dep:libloading"]
+default = ["cuda"]
 jemalloc = ["dep:tikv-jemallocator"]
 malliavin = []
 mimalloc = ["dep:mimalloc"]
@@ -86,5 +87,17 @@ debug = false
 codegen-units = 1
 lto = true
 
+[target.'cfg(target_os = "macos")'.dependencies]
+ndarray-linalg = { version = "0.17.0", features = ["openblas-static"] }
+
+[target.'cfg(target_os = "macos")'.features]
+default = ["jemalloc"]
+
+[target.'cfg(target_os = "linux")'.features]
+default = ["jemalloc"]
+
+[target.'cfg(target_os = "windows")'.features]
+default = ["mimalloc"]
+
 # [package.metadata.docs.rs]
 # all-features = true

src/ai.rs

@@ -1,5 +1,5 @@
 use candle_core::Tensor;
-pub use itransformer::ITransformer;
+// pub use itransformer::ITransformer;
 
 pub mod fou;
 pub mod utils;

src/stochastic/cuda/fgn.cu

@@ -0,0 +1,106 @@
+#include <stdio.h>
+#include <cuda_runtime.h>
+#include <curand_kernel.h>
+#include <cufft.h>
+#include <cuComplex.h>
+#include <math.h>
+
+#ifdef _WIN32
+#define EXPORT __declspec(dllexport)
+#else
+#define EXPORT
+#endif
+
+__global__ void fill_random_with_eigs(
+    cuComplex *d_data,
+    const cuComplex *d_sqrt_eigs,
+    int traj_size,
+    int m,
+    unsigned long seed)
+{
+    int traj_id = blockIdx.x;
+    if (traj_id >= m)
+        return;
+
+    int idx = threadIdx.x;
+    if (idx >= traj_size)
+        return;
+
+    int data_idx = traj_id * traj_size + idx;
+
+    curandState state;
+    curand_init(seed + traj_id, idx, 0, &state);
+
+    float re = curand_normal(&state);
+    float im = curand_normal(&state);
+    cuComplex noise = make_cuComplex(re, im);
+
+    d_data[data_idx] = cuCmulf(noise, d_sqrt_eigs[idx]);
+}
+
+__global__ void scale_and_copy_to_output(
+    const cuComplex *d_data,
+    float *d_output,
+    int n,
+    int m,
+    int offset,
+    float hurst,
+    float t)
+{
+    int traj_id = blockIdx.x;
+    if (traj_id >= m)
+        return;
+
+    int idx = threadIdx.x;
+    int out_size = n - offset;
+    if (idx >= out_size)
+        return;
+
+    int data_idx = traj_id * (2 * n) + (idx + 1);
+    float scale = powf((float)n, -hurst) * powf(t, hurst);
+
+    int out_idx = traj_id * out_size + idx;
+    d_output[out_idx] = d_data[data_idx].x * scale;
+}
+
+extern "C" EXPORT void fgn_kernel(
+    const cuComplex *d_sqrt_eigs,
+    float *d_output,
+    int n,
+    int m,
+    int offset,
+    float hurst,
+    float t,
+    unsigned long seed)
+{
+    int traj_size = 2 * n;
+
+    cuComplex *d_data = nullptr;
+    cudaMalloc(&d_data, (size_t)m * traj_size * sizeof(cuComplex));
+
+    {
+        dim3 gridDim(m);
+        dim3 blockDim(traj_size);
+        fill_random_with_eigs<<<gridDim, blockDim>>>(
+            d_data, d_sqrt_eigs, traj_size, m, seed);
+        cudaDeviceSynchronize();
+    }
+
+    {
+        cufftHandle plan;
+        cufftPlan1d(&plan, traj_size, CUFFT_C2C, m);
+        cufftExecC2C(plan, d_data, d_data, CUFFT_FORWARD);
+        cudaDeviceSynchronize();
+        cufftDestroy(plan);
+    }
+
+    {
+        dim3 gridDim(m);
+        dim3 blockDim(n);
+        scale_and_copy_to_output<<<gridDim, blockDim>>>(
+            d_data, d_output, n, m, offset, hurst, t);
+        cudaDeviceSynchronize();
+    }
+
+    cudaFree(d_data);
+}