//[dependencies]
//ndarray = ""
//rayon = "1.0.3"
//ndarray-parallel = "0.9.0"


#![allow(non_snake_case)]
#![allow(unused)]

#[macro_use]
extern crate ndarray;

use ndarray::prelude::*;
use ndarray::Zip;

//extern crate ndarray_parallel;
//use ndarray_parallel::prelude::*;

extern crate rayon;
use rayon::prelude::*;

use std::time::{Duration, Instant};


fn f(x: f64, x0: f64, dx: f64) -> f64 {
    let c = (x - x0) / dx;
    1.0 / (c * c + 1.0)
}

fn serial() {
    let nx = 100000;
    let xmin = 0.0;
    let xmax = 1.0;
    let Dx = (xmax - xmin) / 20.0;
    let x : Array1<f64> = Array::linspace(xmin, xmax, nx);
    let dx = x[1] - x[0];
    let mut y : Array1<f64> = Array1::<f64>::zeros(nx);
    let cc = 2.3;

    for i in 0..nx {
        let x0 = x[i];
        let di = std::cmp::max(2, (Dx / dx * 10.0) as usize);
        let i1 = std::cmp::max(0, i - di);
        let i2 = std::cmp::max(nx - 1, i + di);

        for j in i1..i2 {
            y[i] = y[i] + cc * f(x[i], x0, Dx);
        }
    }
    println!("{}", y.sum());
}

fn parallel() {
    let nx = 100000;
    let xmin = 0.0;
    let xmax = 1.0;
    let Dx = (xmax - xmin) / 20.0;
    let x : Array1<f64> = Array::linspace(xmin, xmax, nx);
    let dx = x[1] - x[0];
    let mut y : Array1<f64> = Array1::<f64>::zeros(nx);
    let cc = 2.3;

    for i in 0..nx {
        let x0 = x[i];
        let di = std::cmp::max(2, (Dx / dx * 10.0) as usize);
        let i1 = std::cmp::max(0, i - di);
        let i2 = std::cmp::min(nx - 1, i + di);

        let mut xc  = x.slice(s![i1..i2]);
        let mut yc  = y.slice(s![i1..i2]);

        Zip::from(&mut yc).and(&xc).par_apply(|y_, &x_| {*y_ = *y_ + cc * f(x_, x0, Dx);});

//error[E0277]: the trait bound `&mut ndarray::ArrayBase<ndarray::ViewRepr<&f64>, ndarray::Dim<[usize; 1]>>: ndarray::IntoNdProducer` is not satisfied
//  --> src/main.rs:62:9
//   |
//62 |         Zip::from(&mut yc).and(&xc).par_apply(|y_, &x_| {*y_ = *y_ + cc * f(x_, x0, Dx);});
//   |         ^^^^^^^^^ the trait `ndarray::IntoNdProducer` is not implemented for `&mut ndarray::ArrayBase<ndarray::ViewRepr<&f64>, ndarray::Dim<[usize; 1]>>`
//   |
//   = help: the following implementations were found:
//             <&'a mut ndarray::ArrayBase<S, D> as ndarray::IntoNdProducer>
//             <&'a ndarray::ArrayBase<S, D> as ndarray::IntoNdProducer>
//   = note: required by `ndarray::Zip::<(P,), D>::from`

    }

    println!("{}", y.sum());
}

fn main() {
    let t1 = Instant::now();
    serial();
    let t2 = Instant::now();
    parallel();
    let t3 = Instant::now();

    println!("{:?}  {:?}", t2.duration_since(t1), t3.duration_since(t2));
}