// export POISSON_INTENSITY=0.000100 export JUMP_SIZE=0.002 export P_00=0.595000 export P_11=0.3 export GBM_DRIFT=-0.0545 export GBM_SIGMA=0.038 && bash monte_carlo.sh
// POISSON_INTENSITY=1 POISSON_STD_DEV=1 NUM_STEPS=86400 GBM_DRIFT=0.00125 GBM_SIGMA=0.05 OU_SIGMA=0.05 OU_MEAN=1600 OU_SPEED=1 OU=2 INIT_X=1000 INIT_Y=1597000 BASISPOINT=100 BLOCKTIME=12; while true; do ./sim $POISSON_INTENSITY $POISSON_STD_DEV $NUM_STEPS $GBM_DRIFT $GBM_SIGMA $OU_SPEED $OU_MEAN $OU_SIGMA $OU $INIT_X $INIT_Y $BASISPOINT $BLOCKTIME; sleep 1; done 
// ########################################################################################
// ####  Simulation of Automated Market Makers                                #############
// ########################################################################################
// #### Developed by Bence Ladoczki <ladoczki.bence@vik.bme.hu> 2023-2024  ################
// ####               All rights reserved                           #######################
// ########################################################################################
// # install:  sudo apt-get install libgsl-dev && gcc -DLOGLEVEL=3 -o sim sim.c power_iteration.c -lgsl -lgslcblas -lm -fopenmp
// # run:  gcc -DLOGLEVEL=1 -o sim sim.c -lm && ./sim POISSON_INTENSITY POISSON_STD_DEV NUM_STEPS GBM_DRIFT GBM_SIGMA <OU_SPEED> <OU_MEAN> <OU_SIGMA> <Ornstein-Uhlenbeck> <INIT_X> <INIT_Y> <BASISPOINT> <BLOCKTIME>
// # run:  gcc -DLOGLEVEL=3 -o sim sim.c -lm && ./sim 0 0 100000 0.0 0.02 1
// # run:  gcc -DLOGLEVEL=3 -o sim sim.c -lm && ./sim 0 0 100000 0.0 0.02 1
// # run:  gcc -DLOGLEVEL=3 -o sim sim.c -lm && ./sim 0 0 86400 0.00125 0.05 11 1 1 1 0 1000 1597000 100 12
// ########################################################################################
// ########################################################################################
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>  // vagy C++-ban: #include <cstdint>
#include <string.h>
#include <math.h>
#include <time.h>
#include <stdbool.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "power_iteration.h"
#include <omp.h>  // OpenMP fejléc

#ifndef LOGLEVEL
#define LOGLEVEL 1
#endif
#define MAX_STEP_SIZE 10.0
#define EPSILON 1e+4
#define SECONDS_IN_A_DAY 86400
//#define INIT_X 1000 // ETH
//#define INIT_Y 1597000 // USDT

#define EPSILON_2 1e-9  // A numerikus pontosság
#define EPSILON_3 1e-3
#define MAX_ITER 1000  // Maximális iterációk száma

// Célfüggvény: f(theta) = P_{CEX} - (Delta y / Delta x)
double objective_buy(double theta, double L, double x1, double x, double P_CEX) {
    double y1 = pow(L / pow(x1, theta), 1.0 / (1.0 - theta));
    double y2 = pow(L / pow(x1+x, theta), 1.0 / (1.0 - theta));
    double delta_y = y1 - y2;

    return (delta_y / x) - P_CEX;
}
double objective_sell(double theta, double L, double y1, double y, double P_CEX) {
    double x1 = pow(L / pow(y1, 1.0 - theta), 1.0 / theta);
    double x2 = pow(L / pow(y1+y, 1.0 - theta), 1.0 / theta);
    double delta_x = x1 - x2;
    return (y / delta_x) - P_CEX;
}

// Véletlenszerű előjel beállítása
double apply_random_sign(double value) {
    double sign = (double)rand() / RAND_MAX; // Véletlen szám [0, 1]
    if (sign < 0.5) {
        return  value;  // Pozitív előjel
    } else {
        return -value; // Negatív előjel
    }
}
// Power law mintavételezés függvény
double power_law_sample(double xmin, double alpha) {
 double retval = 0.0;
// while (retval < 0.1){
 double u = (double)rand() / (RAND_MAX + 1.0); // Egyenletes eloszlás [0, 1) -> elkerüljük, hogy u = 1
 retval =  xmin * pow(1 - u, -1.0 / (alpha - 1.0));
// }
 return retval;
}
// Bisection módszer
double find_theta(double L, double x1, double x, double P_CEX, int operand) {
    double low = 0.79;  // Az alsó határ (közel 0, de nem 0)
    double high = 0.1; // A felső határ (közel 1, de nem 1)
	if(operand == 0){
	
    low = 0.1;  // Az alsó határ (közel 0, de nem 0)
    high = 0.79; // A felső határ (közel 1, de nem 1)
	}
    double theta = 0.5; // Kezdeti becslés
    int iter = 0;
    while (fabs(high - low) > EPSILON_2 && iter < MAX_ITER) {
        theta = (low + high) / 2.0;
	double f_theta = 0.0;
	if(operand == 0) f_theta = objective_buy(theta, L, x1, x, P_CEX);
	if(operand == 1) f_theta = objective_sell(theta, L, x1, x, P_CEX);
        if (f_theta > 0) {
            high = theta;  // Szűkítsd az intervallumot
        } else {
            low = theta;
        }
        iter++;
    }
    if (iter == MAX_ITER) {
        printf("Figyelem: A bisection módszer nem konvergált a maximális iteráció alatt!\n");
    }
    return theta;
}

int poisson_random(double lambda) {
    double L = exp(-lambda);
    double p = 1.0;
    int k = 0;

    do {
        k++;
        p *= (double)rand() / RAND_MAX;
    } while (p > L);

    return k - 1;
}

int generate_jump_indicator(double lambda) {
    int jumps = poisson_random(lambda);
    return (jumps > 0) ? 1 : 0; 
}
void calculate_daily_returns(double prices[], size_t num_prices, double daily_returns[]) {
    for (size_t i = 1; i < num_prices; ++i) {
        daily_returns[i - 1] = (prices[i] - prices[i - 1]) / prices[i - 1];
    }
}
double calculate_volatility(double daily_returns[], size_t num_returns) {
    double sum = 0.0;
    for (size_t i = 0; i < num_returns; ++i) {
        sum += daily_returns[i];
    }

    double mean = sum / num_returns;

    double sum_squared_diff = 0.0;
    for (size_t i = 0; i < num_returns; ++i) {
        double diff = daily_returns[i] - mean;
        sum_squared_diff += diff * diff;
    }

    double variance = sum_squared_diff / num_returns;
    double volatility = sqrt(variance);

    return volatility;
}
double randomStep(double std_dev) {
    static int initialized = 0;
    if (!initialized) {
     FILE *fp = fopen("/dev/urandom", "r");
     unsigned int seed;
     fread(&seed, sizeof(seed), 1, fp);
     fclose(fp);
     srand(seed);
        initialized = 1;
    }
    double u1 = (double)rand() / RAND_MAX;
    double u2 = (double)rand() / RAND_MAX;
    double z = sqrt(-2.0 * log(u1)) * cos(2.0 *  3.14159265358979323846  * u2);
    return z*std_dev;
}
int areEqual(double a, double b) {
    return fabs(a - b) < EPSILON;
}
typedef struct {
// The DEX with Weighted Geometric Mean Market Maker
    double asset_x;
    double asset_y;
    double l_value;
    double theta;
    double (*bfunc)(); 
    double (*gqx)(); 
    double (*gqy)(); 
    double (*gev)(); 
    int (*executetrade_y_to_x)(); 
    int (*executetrade_x_to_y)(); 
    int (*ix)(); 
    int (*iy)(); 
    double (*gpr)(); 
    double (*gpv)(); 
    double prc;
} wgmmm_dex;
typedef struct {
// The DEX with Constant Product Market Maker
    double asset_x;
    double asset_y;
    double l_value;
    double (*bfunc)(); 
    double (*gqx)(); 
    double (*gqy)(); 
    double (*gev)(); 
    int (*executetrade_y_to_x)(); 
    int (*executetrade_x_to_y)(); 
    int (*ix)(); 
    int (*iy)(); 
    double (*gpr)(); 
    double (*gpv)(); 
    double prc;
} cpmmm_dex;
double wgmmm(wgmmm_dex* dex){
// Weighted Geometric Mean Market Maker
    double x = dex->asset_x;
    double y = dex->asset_y;
    double theta = dex->theta;
    return pow(x, theta) * pow(y, 1 - theta);
}
double cpmmm(cpmmm_dex* dex){
// Constant Product Market Maker
    double x = dex->asset_x;
    double y = dex->asset_y;
    return x*y;
}
int wgmmm_executetrade_y_to_x(wgmmm_dex* dex,double x, double y){
//   printf(" wgmmm_executetrade_y_to_x: %.4f,%.4f,%.4f \n", y/x,x,y);
    double l1 = pow(dex->asset_x,dex->theta)*pow(dex->asset_y,1.0-dex->theta);
    double l2 = pow(dex->asset_x - x,dex->theta)*pow(dex->asset_y+y,1.0-dex->theta);


if (fabs(l1 - l2) > EPSILON) {
    printf("Error: x != deltax in wgmmm_executetrade_y_to_x\n");
    printf("l1: %f, l2: %f\n", l1, l2);
    printf("dex->asset_x: %f, dex->asset_y: %f, dex->theta: %f\n", dex->asset_x, dex->asset_y, dex->theta);
    printf("x: %f, y: %f\n", x, y);
    exit(1);
}
    dex->asset_x = dex->asset_x - x; 
    dex->asset_y = dex->asset_y + y; 
//    if(LOGLEVEL == 1){
//     printf("Exchanged %f x for %f y \n", x, y);
//     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
//    }
    return 0;
}
int wgmmm_executetrade_x_to_y(wgmmm_dex* dex,double x, double y){
// Transfer x and receive y
//
//
//   printf(" wgmmm_executetrade_x_to_y: %.4f,%.4f,%.4f \n", y/x,x,y);
    double l1 = pow(dex->asset_x,dex->theta)*pow(dex->asset_y,1.0-dex->theta);
    double l2 = pow(dex->asset_x - x,dex->theta)*pow(dex->asset_y+y,1.0-dex->theta);


if (fabs(l1 - l2) > EPSILON) {
    printf("Error: x != deltax in wgmmm_executetrade_x_to_y\n");
    printf("l1: %f, l2: %f\n", l1, l2);
    printf("dex->asset_x: %f, dex->asset_y: %f, dex->theta: %f\n", dex->asset_x, dex->asset_y, dex->theta);
    printf("x: %f, y: %f\n", x, y);
    exit(1);
}
    dex->asset_x = dex->asset_x + x; 
    dex->asset_y = dex->asset_y - y; 
//    if(LOGLEVEL == 1){
//     printf("Exchanged %f x for %f y \n", x, y);
//     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
//    }
    return 0;
}
int cpmmm_executetrade_y_to_x(cpmmm_dex* dex,double x, double y){
// Transfer y and receive x
    double newreservex = dex->asset_x - x;
    double newreservey = dex->l_value / newreservex;
    double deltay = newreservey - dex->asset_y;
    if(y != deltay){
     printf(" y != deltay in cpmmm_executetrade_y_to_x %f %f \n",y, deltay);
     exit(1);
    }
    dex->asset_x = dex->asset_x - x; 
    dex->asset_y = dex->asset_y + y; 
//    if(LOGLEVEL == 1){
//     printf("Exchanged %f x for %f y \n", x, y);
//     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
//    }
    return 0;
}
int cpmmm_executetrade_x_to_y(cpmmm_dex* dex,double x, double y){
// Transfer x and receive y
    double newreservey = dex->asset_y - y;
    double newreservex = (double) dex->l_value / newreservey;
    double deltax = newreservex -  dex->asset_x;
    if(x != deltax){
     printf(" x != deltax in cpmmm_executetrade_x_to_y %f %f \n",x, deltax);
     exit(1);
    }
    dex->asset_x = dex->asset_x + x; 
    dex->asset_y = dex->asset_y - y; 
//    if(LOGLEVEL == 1){
//     printf("Exchanged %f x for %f y \n", x, y);
//     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
//    }
    return 0;
}
int cpmmm_inject_x(cpmmm_dex* dex,double x){
// Inject x into the pool
    dex->asset_x = dex->asset_x + x; 
    if(LOGLEVEL == 1){
     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
    }
    return 0;
}
int cpmmm_inject_y(cpmmm_dex* dex,double y){
// Inject y into the pool
    dex->asset_y = dex->asset_y + y; 
    if(LOGLEVEL == 1){
     printf("New state %f x, %f y, price: %f \n", dex->asset_x, dex->asset_y, (double) dex->asset_y / (double) dex->asset_x);
    }
    return 0;
}
double wgmmm_get_price(wgmmm_dex* dex){
// Get the price of asset x denominated in asset y
//  printf(" wgmmm_get_price: %.2f \n", dex->gqx(dex,0.001) / 0.001);
  return dex->gqx(dex,0.001) / 0.001;
}
double cpmmm_get_price(cpmmm_dex* dex){
// Get the price of asset x denominated in asset y
    return (double) dex->asset_y / (double) dex->asset_x;
}
double wgmmm_get_pool_value(wgmmm_dex* dex,double price){
// Get the V_(P)
    return (double) dex->asset_y + (double) dex->asset_x*price;
}
double cpmmm_get_pool_value(cpmmm_dex* dex,double price){
// Get the V_(P)
    return (double) dex->asset_y + (double) dex->asset_x*price;
}

double wgmmm_get_extractable_value(wgmmm_dex* dex,double p1,double p2,int yorx){
// x1^{theta}*y1^{1-theta} = L
// x2^{theta}*y2^{1-theta} = L
// y1/x1 = P1
// y2/x2 = P2
// l/(x2**2) = P2
// x2 = sqrt(l/)/p2^{1-theta}
// y2 = sqrt(l)*p2^{theta}
//  printf("\n %d \n", yorx);
  double retval = 0.0;
	if (yorx == 0) {
         retval =  dex->asset_x-sqrt(dex->asset_y*dex->asset_x)/pow(p2,1.0-dex->theta);
	}
	if (yorx == 1) {
         retval = dex->asset_y-sqrt(dex->asset_y*dex->asset_x)*pow(p2,dex->theta);
	}
//  printf("\n %.2f \n", retval);
	return retval;
}
double cpmmm_get_extractable_value(cpmmm_dex* dex,double p1,double p2,int yorx){
// x1*y1 = L
// x2*y2 = L
// y1/x1 = P1
// y2/x2 = P2
// l/(x2**2) = P2
// x2 = sqrt(l/p2)
// y2 = sqrt(p2*l)
	if (yorx == 0) {
         return dex->asset_x-sqrt(dex->asset_y*dex->asset_x/p2);
	}
	if (yorx == 1) {
         return dex->asset_y-sqrt(p2*(dex->asset_y*dex->asset_x));
	}
	return 0;
}
double cpmmm_getquote_forx(cpmmm_dex* dex,double x){
    double newreservex = dex->asset_x - x;
    double newreservey = dex->l_value / newreservex;
    double deltay = newreservey - dex->asset_y;
    return deltay;
}
double cpmmm_getquote_fory(cpmmm_dex* dex,double y){
    double newreservey = dex->asset_y - y;
    double newreservex = dex->l_value / newreservey;
    double deltax = newreservex -  dex->asset_x;
    return deltax;
}
double wgmmm_getquote_forx(wgmmm_dex* dex,double x){
    double y1 =	pow(dex->l_value / pow(dex->asset_x, dex->theta), 1.0 / (1.0 - dex->theta));
    double y2 =	pow(dex->l_value / pow(dex->asset_x+x, dex->theta), 1.0 / (1.0 - dex->theta));
    double deltay = y1-y2;
//     printf(" deltay: %.2f\n", deltay);
    return deltay;
}
double wgmmm_getquote_fory(wgmmm_dex* dex,double y){
	double x1 = pow(dex->l_value / pow(dex->asset_y, 1.0 - dex->theta), 1.0 / dex->theta);
	double x2 = pow(dex->l_value / pow(dex->asset_y+y, 1.0 - dex->theta), 1.0 / dex->theta);
    double deltax =  x1 - x2;
//     printf(" deltax: %.2f\n", deltax);
    return deltax;
}
int poissonProcess(double lambda) {
    int count = 0;
    double L = exp(-lambda);
    double p = 1.0;

    do {
        count++;
        p *= (double)rand() / RAND_MAX;
    } while (p > L);

    return count - 1;
}
int main(int argc, char *argv[]) {

    double y;
    double x; 
    double exponent;
    if (argc < 13) {
        fprintf(stderr, "Hasznalat: %s sim POISSON_INTENSITY POISSON_STD_DEV NUM_STEPS GBM_DRIFT GBM_SIGMA <OU_SPEED> <OU_MEAN>   <OU_SIGMA> <Ornstein-Uhlenbeck> <INIT_X> <INIT_Y> <BASISPOINT>\n", argv[0]);
        return 1;  // Hibakod
    }
//    printf("Entering sim...");
    double POISSON_INTENSITY = atof(argv[1]);
    double POISSON_STD_DEV   = atof(argv[2]);
    long long int NUM_STEPS = strtoll(argv[3], NULL, 10);
    double GBM_DRIFT         = atof(argv[4]);
    double GBM_SIGMA         = atof(argv[5]);
    double OU_SPEED          = atof(argv[6]);
    double OU_MEAN           = atof(argv[7]);
    double OU_SIGMA          = atof(argv[8]);
    int OU                   = atoi(argv[9]);
    int INIT_X               = atoi(argv[10]);
//    int INIT_Y               = atoi(argv[11]);
    uint64_t INIT_Y = (uint64_t)INIT_X * 4572;
    int BASISPOINT           = atoi(argv[12]);
    int BLOCKTIME            = atoi(argv[13]);
    int NUM_STEPS_OUTER      = atoi(argv[14]);
    double PRED_SIGMA           = atof(argv[15]);
    double JUMP_SIZE           = atof(argv[16]);
    double P_00           = atof(argv[17]);
    double P_11           = atof(argv[18]);
    double DRIFT_FACTOR   = atof(argv[19]);
    double SIGMA_FACTOR   = atof(argv[20]);

//    int matrix_size = 40000;

//    int **matrix = (int **)malloc(matrix_size * sizeof(int *));
//    if (matrix == NULL) {
//        printf("Memóriafoglalási hiba!\n");
//        return 1;
//    }
//    for (int i = 0; i < matrix_size; i++) {
//     matrix[i] = (int *)malloc(matrix_size * sizeof(int));
//     if (matrix[i] == NULL) {
 ////     printf("Memóriafoglalási hiba!\n");
 //     return 1;
 //    }
//    }
    // Minden egyes sor nullázása
 //   for (int i = 0; i < matrix_size; i++) {
//        memset(matrix[i], 0, matrix_size * sizeof(int));
//    }
//    double **markov_matrix = (double **)malloc(matrix_size * sizeof(double *));
//    for (int i = 0; i < matrix_size; i++) {
//     markov_matrix[i] = (double *)malloc(matrix_size * sizeof(double));
//    }
 //   size_t total_size_bytes = matrix_size * matrix_size * sizeof(int);
 //   double total_size_mb = total_size_bytes / (1024.0 * 1024.0);
 //   printf("A mátrix teljes mérete: %.2f MB\n", total_size_mb);

    GBM_SIGMA = GBM_SIGMA / sqrt((double) NUM_STEPS);
    PRED_SIGMA = PRED_SIGMA / sqrt((double) NUM_STEPS);
    GBM_DRIFT = GBM_DRIFT / (double) NUM_STEPS ;
//    GBM_DRIFT = GBM_SIGMA*GBM_SIGMA/2;
//    GBM_DRIFT = 0.0;
    if(LOGLEVEL == 3){
     printf("Now simulating with LVR in an AMM with \n");
     printf("POISSON_INTENSITY of the jump process: %.1f\n", POISSON_INTENSITY);
     printf("STD_DEV of the jump amplitudes: %.1f\n", POISSON_STD_DEV);
     printf("Number of Monte Carlo steps: %lld\n", NUM_STEPS);
     printf("Number of outer Monte Carlo steps: %d\n", NUM_STEPS_OUTER);
     printf("GBM Drift: %.15f\n", GBM_DRIFT);
     printf("Volatility: %f\n", GBM_SIGMA);
     printf("Pred Volatility: %f\n", PRED_SIGMA);
     printf("OU_SPEED: %f\n", OU_SPEED);
     printf("OU_MEAN: %f\n", OU_MEAN);
     printf("OU_SIGMA: %f\n", OU_SIGMA);
     printf("Ornstein-Uhlenbeck: %d\n", OU);
     printf("INIT_X: %d\n", INIT_X);
     printf("INIT_Y: %llu\n", (unsigned long long) INIT_Y);
     printf("BASISPOINT: %d\n", BASISPOINT);
     printf("BLOCKTIME: %d\n", BLOCKTIME);
    }

    wgmmm_dex *myDex = (wgmmm_dex *)malloc(sizeof(wgmmm_dex));
    myDex->asset_x = INIT_X;
    myDex->asset_y = INIT_Y;
    myDex->theta = 0.5;
    myDex->bfunc = wgmmm;
    myDex->l_value = myDex->bfunc(myDex);
    myDex->gqx = wgmmm_getquote_forx;
    myDex->gqy = wgmmm_getquote_fory;
    myDex->gev = wgmmm_get_extractable_value;
    myDex->executetrade_y_to_x = wgmmm_executetrade_y_to_x;
    myDex->executetrade_x_to_y = wgmmm_executetrade_x_to_y;
    myDex->gpr = wgmmm_get_price;
    myDex->gpv = wgmmm_get_pool_value;
//    printf("%f\n", myDex->gqx(myDex,1));

    cpmmm_dex *myDex2 = (cpmmm_dex *)malloc(sizeof(cpmmm_dex));
    
    myDex2->asset_x = INIT_X;
    myDex2->asset_y = INIT_Y;
    myDex2->l_value = INIT_X*INIT_Y;
    myDex2->bfunc = cpmmm;
    myDex2->gqx = cpmmm_getquote_forx;
    myDex2->gqy = cpmmm_getquote_fory;
    myDex2->gev = cpmmm_get_extractable_value;
    myDex2->executetrade_y_to_x = cpmmm_executetrade_y_to_x;
    myDex2->executetrade_x_to_y = cpmmm_executetrade_x_to_y;
    myDex2->ix = cpmmm_inject_x;
    myDex2->iy = cpmmm_inject_y;
    myDex2->gpr = cpmmm_get_price;
    myDex2->gpv = cpmmm_get_pool_value;
//    printf("%f\n", myDex2->gqx(myDex2,3));
//    double x = 3;
//    double y = 5000;
//    for(int i=1; i <= 100; i++) myDex2->executetrade_y_to_x(myDex2,x, myDex2->gqx(myDex2,x) );
//    for(int i=1; i <= 100; i++) myDex2->executetrade_x_to_y(myDex2,myDex2->gqy(myDex2,y) , y );


    if(LOGLEVEL == 3){
    printf("Timestamp,Price_at_pga_end,Final_price,Pool_price,Pool_action,LVR_action,Exchanged_ETH,Pool_Value,LVR_Value,LVR_balance_vs_pool_balance,Number_of_jumps,Realized_price_difference,Expected_price_difference,Inst_profit,Total_profit_ideal,Total_profit_realistic,Difference\n");
    }
    srand(time(NULL)); // Random szám generátor inicializálása

    double vp_lvr;
    double vp_pool;
    double alpha = 0.5;  // Power law paraméter
    double xmin = 5.0;  // Alsó határ

    double spotprice = (double) INIT_Y / (double) INIT_X;

    double wiener_step = 0.0;
    double lvr_eth = INIT_X;
    double lvr_usdt = INIT_Y;
//    double* steps2 = (double*)malloc(sizeof(double)*NUM_STEPS);
//    double* steps = (double*)malloc(sizeof(double)*NUM_STEPS);
    double timestep = 1.0;
    int number_of_trades = 0;
    remove("data/missprice.csv");
    remove("data/theta.csv");
    remove("data/cex_price.csv");
    remove("data/dex_price.csv");
    remove("data/learning_data_file.csv");
    remove("data/profits.csv");
    remove("data/fees.csv");
    FILE *theta_file = fopen("data/theta.csv", "a");
    FILE *missprice_file = fopen("data/missprice.csv", "a");
    FILE *cex_price_file = fopen("data/cex_price.csv", "a");
    FILE *dex_price_file = fopen("data/dex_price.csv", "a");
    FILE *learning_data_file = fopen("data/learning_data_file.csv", "a");
    FILE *profits_file = fopen("data/profits.csv", "a");
    FILE *fees_file = fopen("data/fees.csv", "a");
    int previous_missprice = 0;

    fprintf(learning_data_file, "X,theta_old,usdt,eth,trade_volume,theta_delta,theta_new\n");
    fprintf(theta_file, "X,Theta,Spotprice,DEX_price,Pool_value,LVR,ETH,Estimation_error,USDT,amount\n");
    clock_t start_time = clock();

    double profit = 0.0;
    for (int step2 = 1; step2 <= NUM_STEPS_OUTER; step2++) {
//       if (step2 % (NUM_STEPS_OUTER / 10) == 0) {
//            clock_t current_time = clock();
//            double elapsed_time = (double)(current_time - start_time) / CLOCKS_PER_SEC;
//            printf("Státusz: %d%% kész, idő: %.2f másodperc\n", 
//              (step2 * 100) / NUM_STEPS_OUTER, elapsed_time);
//            fflush(stdout);
//            start_time = current_time;
//       }
	    int cs = 0;
            int ds = 1;
            double jd = GBM_DRIFT *  DRIFT_FACTOR;
            double jump_sigma = GBM_SIGMA * SIGMA_FACTOR;
            double P[2][2] = {{P_00, 1.0-P_00}, {1.0-P_11, P_11}};
            for (long long int step = 1; step <= NUM_STEPS; step++) {
                double w = ((double) rand() / RAND_MAX) * 2.0 - 1.0;
                double ex;
                if (cs == 0) {
                    ex = (GBM_DRIFT - 0.5 * GBM_SIGMA * GBM_SIGMA) * timestep + GBM_SIGMA * sqrt(timestep) * w;
                } else {
                    double r = (double) rand() / RAND_MAX;
                    int ji = (r < POISSON_INTENSITY) ? 1 : 0;
                    int js = ds;
//                    int js = (rand() % 2 == 0) ? ds : -ds;
                    const double jsz = JUMP_SIZE;
                    ex = (jd * ds - 0.5 * jump_sigma * jump_sigma) * timestep + jump_sigma * sqrt(timestep) * w + js * jsz * ji;
                }
                spotprice *= exp(ex);
                double r = (double) rand() / RAND_MAX;
                int ps = cs;
                if (cs == 0) {
                    cs = (r < P[0][0]) ? 0 : 1;
                } else {
                    cs = (r < P[1][0]) ? 0 : 1;
                }
                if (ps != cs) ds = -ds;

      double basispoint = BASISPOINT;
      double percentage = basispoint*0.01;
      double dex_price = myDex->gpr(myDex);
      double wiener_step2 = randomStep(1.0);
      dex_price = myDex->gpr(myDex);
      double theta_old = myDex->theta;
//      for (int i = -500000; i <= 500000; i++) {
//        myDex->theta = theta_old + i / 1000000000.0;
//        if (fabs(myDex->gpr(myDex) - dex_price) < EPSILON_3) {
////	    printf("%.6f,%.6f,%.6f\n",myDex->gpr(myDex),myDex->theta,dex_price);
//            break;
//        }
//      }
//    printf("Final absolute difference: %.6f,%.6f,%.6f\n", myDex->gpr(myDex) , dex_price, myDex->theta);
    double missprice = log(spotprice/dex_price);
    fprintf(missprice_file, "%lld,%f\n",step, missprice);
    fprintf(cex_price_file, "%lld,%f\n",step, spotprice);
    fprintf(dex_price_file, "%lld,%f\n",step, dex_price);
    if(fabs(missprice) > basispoint/10000 && step % BLOCKTIME == 0){
     if(step % BLOCKTIME == 0){
       if(LOGLEVEL == 3){ printf("%lld,%.2f,%.2f,", step, spotprice, dex_price); }
       if( spotprice < myDex->gpr(myDex)){
// If the price is lower on a CEX, buy the asset on it and sell it on the DEX
        y = myDex->gev(myDex,dex_price,spotprice*1.003004505,1);
//        printf(" mydex y: %.2f \n",  y );
        x = myDex->gqy(myDex,y); 
        y = x*spotprice;
//	 printf(" effective_price: %.2f,%.2f,%.2f \n", effective_price, x,myDex->gqx(myDex,x));
//        printf(" spotprice: %.2f,%.2f \n", spotprice, myDex->gpr(myDex) );
        lvr_eth = lvr_eth + x;
        lvr_usdt = lvr_usdt - y;
        vp_lvr = lvr_usdt + lvr_eth*spotprice;
        vp_pool = myDex->gpv(myDex,spotprice);
        if(y > 0.001){
        myDex->executetrade_x_to_y(myDex,myDex->gqy(myDex,y) , y );
        }
        if(LOGLEVEL == 3){
         printf("sell,buy,%.2f,%.3f,%.3f,%.0f,%.0f,%3.f,%.5f",
         		y/x,
         		myDex->gpr(myDex),
         		myDex->gqy(myDex,y),
         		vp_pool,
         		vp_lvr,vp_lvr-vp_pool,myDex->theta);
        }
         number_of_trades++;
	 double increment = fabs(spotprice*1.003004505-dex_price)*fabs(x);
//	 printf("incrementing profit 1: %f\n", increment);
	 profit += increment;
         fprintf(profits_file, "%lld,%f,%f,%f,%f,%f\n",step, fabs(spotprice-dex_price)*fabs(x),x,y,spotprice,dex_price );
         fprintf(fees_file, "%lld,%f\n",step, fabs(y)*percentage/100.0 );
       }
       else { 
        x = myDex->gev(myDex,dex_price,spotprice/1.003004505,0);
        y = myDex->gqx(myDex,x);
        lvr_eth = lvr_eth - x ;
        lvr_usdt = lvr_usdt + x * spotprice;
        vp_lvr = lvr_usdt + lvr_eth*spotprice;
        vp_pool = myDex->gpv(myDex,spotprice);
//        printf(" x: %.2f \n", x  );
//        printf(" spotprice: %.2f,%.2f \n", spotprice, myDex->gpr(myDex) );
        if(x > 0.001){ myDex->executetrade_y_to_x(myDex,x,myDex->gqx(myDex,x)); }
        if(LOGLEVEL == 3){
         printf("buy,sell,%.2f,%.3f,%.3f,%.0f,%.0f,%.3f,%.5f",
         		y/x,
         		myDex->gpr(myDex),
         		myDex->gqy(myDex,y),
         		vp_pool,
         		vp_lvr, vp_lvr-vp_pool,myDex->theta);
         }
         number_of_trades++;
	 double increment = fabs(spotprice/1.003004505-dex_price)*fabs(x);
//	 printf("incrementing profit 2: %f\n", increment);
	 profit += increment;
         fprintf(profits_file, "%lld,%f,%f,%f,%f,%f\n",step, fabs(spotprice-dex_price)*fabs(x),x,y,spotprice,dex_price );
         fprintf(fees_file, "%lld,%f\n",step, fabs(y)*percentage/100.0 );
       }
        if(LOGLEVEL == 3){ printf("\n"); }
      }
      dex_price = myDex->gpr(myDex);
      missprice = log(spotprice/dex_price);
     }
    }
    }

//    printf("%.3f\n",(double) BLOCKTIME*100.0 * (double) number_of_trades / (double) NUM_STEPS);
    printf("number_of_trades: %d, profit: %.3f\n", number_of_trades, profit );
    fclose(missprice_file);
    fclose(theta_file);
    fclose(learning_data_file);
    fclose(profits_file);
    fclose(fees_file);
    fclose(cex_price_file);
    fclose(dex_price_file);
    free(myDex);
    free(myDex2);
//    free(steps);
//    free(steps2);

//    FILE *fp = fopen("data/transition_matrix.csv", "w");
//    if (fp == NULL) {
//     printf("Nem sikerült megnyitni a fájlt!\n");
//     return 1;
//    }
//    int matrix_limit = 1980;
//    for (int i = matrix_limit; i < matrix_size-matrix_limit; i++) {
//     for (int j = matrix_limit; j < matrix_size-matrix_limit; j++) {
//      fprintf(fp, "%d", matrix[i][j]);
//      if (j < matrix_size - 1 - matrix_limit) {
//       fprintf(fp, ","); 
//      }
//     }
//     fprintf(fp, "\n");  
//    }
//    fclose(fp);

//      int max_i = -1, max_j = -1;
//      int min_i = matrix_size, min_j = matrix_size;
      
//    long long non_zero_count = 0;

//    #pragma omp parallel for reduction(max:max_i,max_j) reduction(min:min_i,min_j) reduction(+:non_zero_count)
//    for (int i = 1; i < matrix_size; i++) {
//    double row_sum = 0.0;

 //   int nonzero = 0;
 //   for (int j = 1; j < matrix_size; j++) {
 //       row_sum += matrix[i][j];
 //       if(matrix[i][j] > 0) nonzero++;
 //   }

    // Normalizáljuk a sort
 //   for (int j = 0; j < matrix_size; j++) {
 //       if (row_sum > 0 && nonzero > matrix_size/1000) {  // Elkerüljük az osztást nullával
 //           markov_matrix[i][j] = (double)matrix[i][j] / row_sum;
 ////       } else {
 //           markov_matrix[i][j] = 0.0;  // Ha a sorban nincs átmenet, állítsuk nullára
 //       }

  //      // Ellenőrizzük, hogy nem nulla az érték
  //      if (markov_matrix[i][j] > 0.00001) {
//// printf("markov_matrix[i][j]: %d,%d,%.9f \n",i,j,markov_matrix[i][j]);
  //          non_zero_count++; // Nem nulla elem számlálása

  //          // Kritikus szakasz nélkül, redukcióval kezeljük a min/max értékeket
  //          if (i > max_i) max_i = i;
  //          if (j > max_j) max_j = j;
  //          if (i < min_i) min_i = i;
  //          if (j < min_j) min_j = j;
  //      }
  //  }

  //  // Az adott sor memóriájának felszabadítása
  //  // Minden szál csak a saját iterációja végén próbálja meg felszabadítani a sort
  //  #pragma omp critical
  //  {
  //      free(matrix[i]);
  //  }
  //  }
////    for (int i = 1; i < matrix_size; i++) {
////        double prob_i = 0;
////        for (int j = 1; j < matrix_size; j++) {
////            prob_i += markov_matrix[i][j];
////        }
////        printf("prob_i: %d,%.5f \n",i,prob_i);
////    }

  //  free(matrix);

////    // Sűrűség kiszámítása
  //  double density = (double)non_zero_count / (matrix_size * matrix_size);
  //  printf("A mátrix sűrűsége: %.6f\n", density);
//    printf("Nem nulla elemek száma: %lld\n", non_zero_count);
//      
//      printf("Max i index: %d\n", max_i);
//      printf("Max j index: %d\n", max_j);
//      printf("Min i index: %d\n", min_i);
//      printf("Min j index: %d\n", min_j);
//
//     int truncated_dimension = (max_i-min_i+1 > max_j-min_j+1) ? max_i-min_i+1 : max_j-min_j+1;
//     printf("Truncated dimension: %d\n", truncated_dimension);
//
//     double **trunacted_matrix = (double **)malloc(truncated_dimension * sizeof(double *));
//     if (trunacted_matrix == NULL) {
//         printf("Memóriafoglalási hiba!\n");
//         return 1;
//     }
//     for (int i = 0; i < truncated_dimension; i++) {
//      trunacted_matrix[i] = (double *)malloc(truncated_dimension * sizeof(double));
//      if (trunacted_matrix[i] == NULL) {
//       printf("Memóriafoglalási hiba!\n");
//       return 1;
//      }
//     }
//     // Minden egyes sor nullázása
//     for (int i = 0; i < truncated_dimension; i++) {
//         memset(trunacted_matrix[i], 0, truncated_dimension * sizeof(double));
//     }
//    non_zero_count = 0;
//      for (int i = 0; i < truncated_dimension; i++) {
//          for (int j = 0; j < truncated_dimension; j++) {
//           trunacted_matrix[i][j] = markov_matrix[i+min_i-1][j+min_j-1];
//            if (trunacted_matrix[i][j] > 0.00001) { non_zero_count++; }
////            printf(" %d,%d,%.7f \n",i,j,trunacted_matrix[i][j]);
//          }
//      }
//    density = (double)non_zero_count / (truncated_dimension * truncated_dimension);
//    printf("A mátrix sűrűsége: %.6f\n", density);
//    printf("Nem nulla elemek száma: %lld\n", non_zero_count);
//// GSL mátrix létrehozása és feltöltése
//    gsl_matrix *A = gsl_matrix_alloc(matrix_size, matrix_size);
//    for (int i = 0; i < matrix_size; i++) {
//        for (int j = 0; j < matrix_size; j++) {
//            gsl_matrix_set(A, i, j, markov_matrix[i][j]);
//        }
//    }
//
//    // EigenSolver létrehozása
//    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(matrix_size);
//    gsl_vector *eval = gsl_vector_alloc(matrix_size); // Sajátértékek
//    gsl_matrix *evec = gsl_matrix_alloc(matrix_size, matrix_size); // Sajátvektorok
//
//    // Sajátértékek számítása
//    gsl_eigen_symmv(A, eval, evec, w);
//
//    // Nem nulla sajátértékek kiírása
//    printf("Nem nulla sajátértékek:\n");
//    for (int i = 0; i < matrix_size; i++) {
//        double eigenvalue = gsl_vector_get(eval, i);
//        if (eigenvalue != 0) {
//            printf("Sajátérték %d: %f\n", i + 1, eigenvalue);
//        }
//    }
//
//    // Memória felszabadítása
//    gsl_matrix_free(A);
//    gsl_matrix_free(evec);
//    gsl_vector_free(eval);
//    gsl_eigen_symmv_free(w);

    // Kezdeti eloszlás (pl. egyenletes eloszlás)
//    double *stationary_distribution = (double *)malloc(matrix_size * sizeof(double));
//    #pragma omp parallel for
//    for (int i = 0; i < matrix_size; i++) {
//        stationary_distribution[i] = 1.0 / matrix_size;
//    }
//    power_iteration(markov_matrix, stationary_distribution, matrix_size);
//    FILE* fp = fopen("data/stationary_distribution.csv", "w");
//    for (int i = 0; i < matrix_size; i++) {
//	    if(fabs(stationary_distribution[i]) > 0.00001) fprintf(fp, "%d,%f\n ",i, stationary_distribution[i]);
//    }
//    fclose(fp);
//    system("Rscript plot_stationary_distribution.R");
    return 0;

//    remove("create_missprice_histograms.R");
//    FILE *histogram_file = fopen("create_missprice_histograms.R", "a");
//    fprintf(histogram_file, "library(ggplot2)\n");
//    fprintf(histogram_file, "library(gridExtra)\n");
//    fprintf(histogram_file, "data_cex <- read.csv('./data/cex_price.csv', header = FALSE, col.names = c('Index', 'Price'))\n");
//    fprintf(histogram_file, "plot_cex <- ggplot(data_cex, aes(x = Index, y = Price)) +\n");
//    fprintf(histogram_file, "  geom_line(color = 'blue') +\n");
//    fprintf(histogram_file, "  geom_point(color = 'blue') +\n");
//    fprintf(histogram_file, "  labs(title = 'CEX Price Over Time', x = 'Index', y = 'Price') +\n");
//    fprintf(histogram_file, "  theme_minimal()\n");
//    fprintf(histogram_file, "data_missprice <- read.csv('data/missprice.csv', header = FALSE, col.names = c('Index', 'Missprice'))\n");
//    fprintf(histogram_file, "histogram_missprice <- ggplot(data_missprice, aes(x = Missprice)) +\n");
//    fprintf(histogram_file, "geom_histogram(binwidth = 0.0001, fill = 'blue', color = 'black', alpha = 0.7) +\n");
//    char plot_title[256];
//    if(!OU){
//     sprintf(plot_title, "Histogram of Missprice GBM sigma=%f", GBM_SIGMA);
//    }
//    else if(OU == 1){
//     sprintf(plot_title, "Histogram of Missprice Ornstein-Uhlenbeck speed =%f", OU_SPEED);
//    }
//    else{
//     sprintf(plot_title, "Histogram of Missprice POISSON_INTENSITY=%f", POISSON_INTENSITY);
//    }
//    fprintf(histogram_file, "  labs(title = '%s', x = 'Missprice', y = 'Frequency') +\n", plot_title);
//    fprintf(histogram_file, "  theme_minimal()\n");
//    char filename[256];
//    if(!OU){
//     sprintf(filename, "./data/missprice_histogram_gbm_sigma_%f.pdf", GBM_SIGMA);
//    }
//    else if(OU == 1){
//     sprintf(filename, "./data/missprice_histogram_ou_speed_%f.pdf", OU_SPEED);
//    }
//    else{
//     sprintf(filename, "./data/missprice_histogram_poission_intensity_%f.pdf", POISSON_INTENSITY);
//    }
//    fprintf(histogram_file, "pdf('%s', width = 12, height = 6)\n",filename);
//    fprintf(histogram_file, "grid.arrange(plot_cex, histogram_missprice, nrow = 1)\n");
//    fprintf(histogram_file, "dev.off()\n");
//    fclose(histogram_file);
//    int result = system("Rscript create_missprice_histograms.R");
//    if (result == -1) {
//     printf("Error executing system command.\n");
//     return 1;
//    }
//    result = system("Rscript plot_cex_price.R");
//    if (result == -1) {
//     printf("Error executing system command.\n");
//     return 1;
//    }

//    result = system("Rscript plot_transition_matrix.R");
//    for (int i = 0; i < matrix_size; i++) {
//        free(matrix[i]);
//        free(markov_matrix[i]);
//    }
//    free(matrix);
    return 0;
}