use core::f64::consts::PI;
use gauss_quad::GaussLegendre;
use std::io;
use std::sync::Arc;
use rayon::prelude::*;
// use scc::DashMap;

use dashmap::DashMap;

use std::{fs::File, io::BufReader};

use crate::gq_storage;

pub type Complex = num::complex::Complex<f64>;

const M_ROH: f64 = 0.7736;
const M_PI: f64 = 0.13957;
pub const S0: f64 = 4.0 * M_PI * M_PI;
const M_K: f64 = 0.496;
const B0: f64 = 1.055;
const B1: f64 = 0.15;
const LAMBDA1: f64 = 1.57;
const LAMBDA2: f64 = -1.96;
const LAMBDA_HIGH: f64 = 10.0;
pub const EPSILON: f64 = 1e-15;
const S_MID_CUTOFF: f64 = 1.42 * 1.42;
const INF: f64 = 1e6;

#[derive(Default)]
pub struct Cache {
    pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>,
    pub delta_lut: DashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
    pub omnes_lut: DashMap<
        u32,
        DashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
        nohash_hasher::BuildNoHashHasher<u32>,
    >, // First key is n, second is s
}

// pub struct Cache {
//     pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), RandomState>,
//     pub delta_lut: DashMap<u64, f64, RandomState>,
//     pub omnes_lut: DashMap<
//         u32,
//         DashMap<u64, Complex, RandomState>,
//         RandomState,
//     >, // First key is n, second is s
// }

// pub struct Cache {
//     pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>)>,
//     pub delta_lut: DashMap<u64, f64>,
//     pub omnes_lut: DashMap<u32, DashMap<u64, Complex>> // First key is n, second is s
// }

impl Cache {
    pub fn from_file(filepath: &str) -> io::Result<Cache> {
        // println!("Cache from_file");
        let file = File::open(filepath)?;
        let mut file = BufReader::new(file);

        let gql = DashMap::with_hasher(Default::default());
        for ent in gq_storage::deserialize(&mut file, |_| true)?
                .into_iter()
                .map(|(roots, weights)| {
                    (
                        u32::try_from(roots.len()).unwrap(),
                        (Arc::from(roots), Arc::from(weights)),
                    )
                })
        {
            let _ = gql.insert(ent.0, ent.1);
        }
        let me = Cache {
            gauss_quad_lut: gql,
            delta_lut: DashMap::with_hasher(Default::default()),
            omnes_lut: DashMap::with_hasher(Default::default()),
        };
        Ok(me)
    }

    #[allow(unused)]
    pub fn from_slice(data: &[u8]) -> io::Result<Cache> {
        let mut cursor = std::io::Cursor::new(data);
        let mut gq = DashMap::with_hasher(Default::default());
        for ent in gq_storage::deserialize(&mut cursor, |_| true)?
                .into_iter()
                .map(|(roots, weights)| {
                    (
                        u32::try_from(roots.len()).unwrap(),
                        (Arc::from(roots), Arc::from(weights)),
                    )
                })
        {
            gq.insert(ent.0, ent.1);
        }
        let me = Cache {
            gauss_quad_lut: gq,
            delta_lut: DashMap::with_hasher(Default::default()),
            omnes_lut: DashMap::with_hasher(Default::default()),
        };
        Ok(me)
    }
}

// derivated values

fn lambda_0() -> f64 {
    delta((2.0 * M_K).powi(2))
}

fn delta_mid_cutoff() -> f64 {
    delta(S_MID_CUTOFF)
}

fn atan_shift(x: f64) -> f64 {
    let atan_std: f64 = x.atan();
    if atan_std < 0.0 {
        atan_std + PI
    } else {
        atan_std
    }
}

fn integrate<G: Fn(f64) -> f64 + Send + Sync>(
    roots: &[f64],
    weights: &[f64],
    a: f64,
    b: f64,
    f: G,
    print_prog: bool,
) -> f64 {
    if roots.len() != weights.len() {
        panic!("Error: roots and weights are of different length");
    }
    println!("Calculating single integration: ");
    let prog = indicatif::ProgressBar::new(roots.len() as u64);
    let prog_ = prog.clone();
    let sum:f64 = (weights.par_iter(), roots.par_iter()).into_par_iter().map(move |(weight,root)| {
        let prog = &prog_;
        // println!("Orello");
        let res = weight * f(root * (b - a) / 2.0 + (a + b) / 2.0);
        if print_prog {
            prog.inc(1);
        }
        res
        // 0.0
    }).sum();
    if print_prog {
        prog.finish();
    }
    sum * (b - a) / 2.0
}

fn integrate_complex<G: FnMut(f64) -> Complex>(
    roots: &[f64],
    weights: &[f64],
    a: f64,
    b: f64,
    mut f: G,
) -> Complex {
    // println!("integrate_complex");
    let mut sum: Complex = Complex::new(0.0, 0.0);
    let mut i: usize = 0;
    if roots.len() != weights.len() {
        panic!("Error: roots and weights are of different length");
    }
    while i < roots.len() {
        sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
        i += 1;
    }
    sum * (b - a) / 2.0
}

fn omega(s: f64) -> f64 {
    let a = s.sqrt();
    let b = (1.05f64.powi(2) - s).sqrt();
    (a - b) / (a + b)
}

fn k(s: f64) -> f64 {
    (s / 4.0 - M_PI.powi(2)).sqrt()
}

pub fn delta_with_lut(cache: &Cache, s: f64) -> f64 {
    // println!("delta_with_lut");
    match cache.delta_lut.get(&s.to_bits()) {
        Some(val) => *val,
        None => {
            let val = delta(s);
            let _ = cache.delta_lut.insert(s.to_bits(), val);
            val
        }
    }
}

pub fn delta(s: f64) -> f64 {
    if s <= 2.0 * M_K {
        atan_shift(
            (s.sqrt() / (2.0 * k(s).powi(3))
                * (M_ROH.powi(2) - s)
                * ((2.0 * M_PI.powi(3)) / (M_ROH.powi(2) * s.sqrt()) + B0 + B1 * omega(s)))
            .powi(-1),
        )
    } else if s <= S_MID_CUTOFF {
        let par = 0.5 * s.sqrt() / M_K - 1.0;
        lambda_0() + LAMBDA1 * par + LAMBDA2 * par.powi(2)
    } else {
        PI + (delta_mid_cutoff() - PI) * (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
            / (LAMBDA_HIGH * LAMBDA_HIGH + s)
    }
}

pub fn omnes_integrand_tan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
    let sub = s_tick.tan();
    (delta_with_lut(cache, sub) - delta_with_lut(cache, s))
        / (sub * (sub - s + Complex::new(0.0, EPSILON)) * s_tick.cos().powi(2))
}

pub fn omnes_integrand_notan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
    // println!("omnes_integrand_notan");
    (delta_with_lut(cache, s_tick) - delta_with_lut(cache, s))
        / (s_tick * (s_tick - s + Complex::new(0.0, EPSILON)))
}

pub fn omnes(cache: &Cache, s: f64, n: u32, use_tan: bool) -> Complex {
    // println!("Omnes");
    if let Some(inner_lut) = cache.omnes_lut.get(&n) {
        if let Some(res) = inner_lut.get(&s.to_bits()) {
            return *res;
        }
    }
    if !cache.omnes_lut.contains_key(&n) {
        let _ = cache
            .omnes_lut
            .insert(n, DashMap::with_hasher(Default::default()));
    }

    // println!("Retrieving roots and weights omnes");
    let roots: Arc<[f64]>;
    let weights: Arc<[f64]>;
    match cache.gauss_quad_lut.get(&n) {
        Some(gq_values) => {
            // println!("Somobomo");
            roots = Arc::clone(&gq_values.0);
            weights = Arc::clone(&gq_values.1);
        }
        None => {
            // println!("Nononono");
            let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
            roots = Arc::from(gq_values.0.as_slice());
            weights = Arc::from(gq_values.1.as_slice());
        }
    }

    // println!("will maybe calc n_omnes");
    if !cache.gauss_quad_lut.contains_key(&n) {
        println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n);
        let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
        let _ = cache
            .gauss_quad_lut
            .insert(n, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
    }

    // println!("Will calculate integrate complex");
    let intgrl = if use_tan {
        integrate_complex(&roots, &weights, S0.atan(), PI / 2.0, |s_tick| {
            omnes_integrand_tan(cache, s_tick, s)
        })
    } else {
        integrate_complex(&roots, &weights, S0, INF, |s_tick| {
            omnes_integrand_notan(cache, s_tick, s)
        })
    };
    let omnes_res = (s / PI * (intgrl)).exp()
        * (S0 / Complex::new(S0 - s, -EPSILON)).powf(delta_with_lut(cache, s) / PI);
    let _ = cache
        .omnes_lut
        .get(&n)
        .unwrap()
        .insert(s.to_bits(), omnes_res);

    omnes_res
}

pub fn phi0_integrand(
    cache: &Cache,
    s_tick: f64,
    s: f64,
    n_omnes: u32,
    use_tan_omnes: bool,
    use_tan_phi0: bool,
    use_xsub: bool,
    _use_reduced_integrand: bool,
    use_cutoff: bool,
    cutoff_power: f64,
    cutoff_s: f64,
    cutoff_factor: Complex,
) -> f64 {
    // println!("phi0_integrand");
    let sub = match (use_tan_phi0, use_xsub) {
        (true, true) => s_tick.tan().powi(2) + S0,
        (true, false) => s_tick.tan(),
        (false, true) => s_tick.powi(2) + S0,
        (false, false) => s_tick,
    };
    // let omnes_s = if use_cutoff && s >= cutoff_s {
    //     cutoff_factor * s.powf(cutoff_power)
    // } else {
    //     omnes(cache, s, n_omnes, use_tan_omnes)
    // };
    let omnes_sub = if use_cutoff && sub >= cutoff_s {
        cutoff_factor * sub.powf(cutoff_power)
        // Complex::new((0.00002f64).powf(cutoff_power), 0.0)
    } else {
        omnes(cache, sub, n_omnes, use_tan_omnes)
    };
    let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes);
    let mut res = (omnes_sub / omnes_s).norm().ln() / (sub * (sub - s));
    res *= match (use_tan_phi0, use_xsub) {
        (true, true) => 2.0 / s_tick.cos().powi(2),
        (true, false) => 1.0 / ((sub - S0).sqrt() * s_tick.cos().powi(2)),
        (false, true) => 2.0,
        (false, false) => 1.0 / (sub - S0).sqrt(),
    };
    res
}

pub fn phi0(
    cache: &Cache,
    s: f64,
    n_omnes: u32,
    n_phi0: u32,
    use_tan_omnes: bool,
    use_tan_phi0: bool,
    use_xsub: bool,
    use_reduced_integrand: bool,
    print_prog: bool,
    use_cutoff: bool,
    cutoff_power: f64,
    cutoff_s: f64
) -> f64 {
    // println!("phi0");
    if !cache.gauss_quad_lut.contains_key(&n_phi0) {
        let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap());
        let _ = cache
            .gauss_quad_lut
            .insert(n_phi0, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
        println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n_phi0);
    }
    let (roots, weights) = {
        let ent = cache.gauss_quad_lut.get(&n_phi0).unwrap();
        (Arc::clone(&ent.0), Arc::clone(&ent.1))
    };

    // println!("Will run omnes");
    let analyt = if use_reduced_integrand {
        omnes(cache, s, n_omnes, use_tan_omnes).norm().ln() / S0.sqrt()
    } else {
        (omnes(cache, s, n_omnes, use_tan_omnes) / omnes(cache, S0, n_omnes, use_tan_omnes))
            .norm()
            .ln()
            / (s * (s - S0).powf(1.5))
    };

    let cutoff_factor = omnes(cache, cutoff_s, n_omnes, use_tan_omnes) * cutoff_s.powf(-cutoff_power);

    let intgrl = integrate(
        &*roots,
        &*weights,
        if use_xsub {
            0.0
        } else if use_tan_phi0 {
            S0.atan()
        } else {
            S0
        },
        if use_tan_phi0 { PI / 2.0 } else { INF },
        |x| {
            phi0_integrand(
                cache,
                x,
                s,
                n_omnes,
                use_tan_omnes,
                use_tan_phi0,
                use_xsub,
                use_reduced_integrand,
                use_cutoff,
                cutoff_power,
                cutoff_s,
                cutoff_factor,
            ) * s
                / PI
        },
        print_prog,
    );
    -(s - S0).sqrt() * (intgrl - analyt)
}

#[allow(unused)]
pub fn bit_pattern_randomness_tester(a: f64, b: f64, n_points: u64) -> Vec<f64> {
    let mut modulos = vec![0.0; n_points as usize];
    let x_values: Vec<u64> = (0..modulos.len() as u32)
        .map(|x| (f64::from(x) * ((b - a) / f64::from(modulos.len() as u32)) + a).to_bits())
        .collect();
    for i in 0..x_values.len() {
        let val = (x_values[i] % n_points) as usize;
        modulos[val] += 1.0;
    }
    modulos
}
