strafe-plot 0.1.1

Statistical plotting for Rust statistics based on R
// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use std::error::Error;

use polars::export::num::Float;

use crate::drawing_coords::DrawingCoords;

fn get_xy_sort_indices(x: &[f64], y: &[f64]) -> Vec<usize> {
    let mut indices = (0..x.len()).collect::<Vec<_>>();
    indices.sort_by(|&i1, &i2| {
        let x_cmp = x[i1].partial_cmp(&x[i2]).unwrap();
        if x_cmp.is_eq() {
            y[i1].partial_cmp(&y[i2]).unwrap()
        } else {
            x_cmp
        }
    });
    indices
}

fn sort_by_indices(x: &mut Vec<f64>, indices: &[usize]) {
    let mut x_ordered = vec![0.0; x.len()];
    for (old_index, new_index) in indices.iter().cloned().enumerate() {
        x_ordered[old_index] = x[new_index];
    }
    *x = x_ordered;
}

pub fn draw_points(
    x: &[f64],
    y: &[f64],
    drawing_coords: &DrawingCoords,
) -> Result<Vec<(f64, f64)>, Box<dyn Error>> {
    let mut x = x.to_vec();
    for x in &mut x {
        if x.is_infinite() {
            if x.is_sign_positive() {
                *x = drawing_coords.x_max + drawing_coords.x_space;
            } else {
                *x = drawing_coords.x_min - drawing_coords.x_space;
            }
        }
    }
    let mut y = y.to_vec();
    for y in &mut y {
        if y.is_infinite() {
            if y.is_sign_positive() {
                *y = drawing_coords.y_max + drawing_coords.y_space;
            } else {
                *y = drawing_coords.y_min - drawing_coords.y_space;
            }
        }
    }

    for i in (0..y.len()).rev() {
        if x[i] < drawing_coords.x_min - drawing_coords.x_space
            || x[i] > drawing_coords.x_max + drawing_coords.x_space
            || y[i] < drawing_coords.y_min - drawing_coords.y_space
            || y[i] > drawing_coords.y_max + drawing_coords.y_space
        {
            x.remove(i);
            y.remove(i);
        }
    }

    let indices = get_xy_sort_indices(&x, &y);
    sort_by_indices(&mut x, &indices);
    sort_by_indices(&mut y, &indices);

    Ok(x.into_iter().zip(y.into_iter()).collect())
}

pub fn draw_line(
    x: &[f64],
    y: &[f64],
    drawing_coords: &DrawingCoords,
    dash: bool,
    dash_length: f64,
    dash_iterations: usize,
) -> Result<Vec<Vec<(f64, f64)>>, Box<dyn Error>> {
    let mut x = x.to_vec();
    for x in &mut x {
        if x.is_infinite() {
            if x.is_sign_positive() {
                *x = drawing_coords.x_max + drawing_coords.x_space;
            } else {
                *x = drawing_coords.x_min - drawing_coords.x_space;
            }
        }
    }
    let mut y = y.to_vec();
    for y in &mut y {
        if y.is_infinite() {
            if y.is_sign_positive() {
                *y = drawing_coords.y_max + drawing_coords.y_space;
            } else {
                *y = drawing_coords.y_min - drawing_coords.y_space;
            }
        }
    }

    let indices = get_xy_sort_indices(&x, &y);
    sort_by_indices(&mut x, &indices);
    sort_by_indices(&mut y, &indices);

    for i in (0..y.len()).rev() {
        if !x[i].is_finite() || !y[i].is_finite() {
            x.remove(i);
            y.remove(i);
        }
    }

    let x_boundary_min = drawing_coords.x_min - drawing_coords.x_space;
    let x_boundary_max = drawing_coords.x_max + drawing_coords.x_space;
    let y_boundary_min = drawing_coords.y_min - drawing_coords.y_space;
    let y_boundary_max = drawing_coords.y_max + drawing_coords.y_space;

    let create_boundary_coord = |coord1: (f64, f64), coord2: (f64, f64)| {
        let slope = (coord2.1 - coord1.1) / (coord2.0 - coord1.0);
        let intercept = coord1.1 - (slope * coord1.0);

        let left_coord = (x_boundary_min, (slope * x_boundary_min) + intercept);
        let right_coord = (x_boundary_max, (slope * x_boundary_max) + intercept);
        let top_coord = ((y_boundary_max - intercept) / slope, y_boundary_max);
        let bottom_coord = ((y_boundary_min - intercept) / slope, y_boundary_min);

        [left_coord, right_coord, top_coord, bottom_coord]
            .into_iter()
            .filter(|&coord| {
                coord.0.is_finite()
                    && coord.1.is_finite()
                    && ((coord1.0 <= coord.0 && coord.0 <= coord2.0)
                        || (coord1.0 >= coord.0 && coord.0 >= coord2.0))
                    && ((coord1.1 <= coord.1 && coord.1 <= coord2.1)
                        || (coord1.1 >= coord.1 && coord.1 >= coord2.1))
            })
            .collect::<Vec<_>>()
    };

    let mut additional_x = Vec::new();
    let mut additional_y = Vec::new();

    for is in (0..y.len()).collect::<Vec<_>>().windows(2) {
        let i1 = is[0];
        let i2 = is[1];

        let coord1 = (x[i1], y[i1]);
        let coord2 = (x[i2], y[i2]);

        for boundary_coord in create_boundary_coord(coord1, coord2) {
            additional_x.push(boundary_coord.0);
            additional_y.push(boundary_coord.1);
        }
    }

    x.append(&mut additional_x);
    y.append(&mut additional_y);

    let mut line_segments_x = vec![Vec::new()];
    let mut line_segments_y = vec![Vec::new()];

    let indices = get_xy_sort_indices(&x, &y);
    sort_by_indices(&mut x, &indices);
    sort_by_indices(&mut y, &indices);

    // Remove outside points
    for i in (0..y.len()).rev() {
        if x[i] < x_boundary_min
            || x[i] > x_boundary_max
            || y[i] < y_boundary_min
            || y[i] > y_boundary_max
        {
            line_segments_x.push(Vec::new());
            line_segments_y.push(Vec::new());
            x.remove(i);
            y.remove(i);
        } else {
            let l = line_segments_x.len() - 1;
            line_segments_x[l].push(x[i]);
            line_segments_y[l].push(y[i]);
        }
    }

    line_segments_x = line_segments_x
        .into_iter()
        .filter(|l| !l.is_empty())
        .collect::<Vec<_>>();
    line_segments_y = line_segments_y
        .into_iter()
        .filter(|l| !l.is_empty())
        .collect::<Vec<_>>();

    let mut ret = Vec::new();
    for (mut x, mut y) in line_segments_x.iter_mut().zip(line_segments_y.iter_mut()) {
        // Remove duplicate points
        let mut i = 0;
        while i + 1 < x.len() {
            while i + 1 < x.len() && x[i] == x[i + 1] && y[i] == y[i + 1] {
                x.remove(i + 1);
                y.remove(i + 1);
            }
            i += 1;
        }

        if dash {
            let x_center = (drawing_coords.x_max + drawing_coords.x_min) / 2.0;
            let mut x_range = drawing_coords.x_max - drawing_coords.x_min;
            if x_range == 0.0 {
                x_range = f64::infinity();
            }
            let mut standard_x = x
                .iter()
                .map(|x_i| (x_i - x_center) / x_range)
                .collect::<Vec<_>>();

            let y_center = (drawing_coords.y_max + drawing_coords.y_min) / 2.0;
            let mut y_range = drawing_coords.y_max - drawing_coords.y_min;
            if y_range == 0.0 {
                y_range = f64::infinity();
            }
            let mut standard_y = y
                .iter()
                .map(|y_i| (y_i - y_center) / y_range)
                .collect::<Vec<_>>();

            for _ in 0..dash_iterations {
                for is in (0..x.len()).collect::<Vec<_>>().windows(2) {
                    if ((standard_x[is[1]] - standard_x[is[0]]).powi(2)
                        + ((standard_y[is[1]] - standard_y[is[0]])
                            * (drawing_coords.chart_height / drawing_coords.chart_width))
                            .powi(2))
                    .sqrt()
                        >= dash_length
                    {
                        x.push((x[is[0]] + x[is[1]]) / 2.0);
                        standard_x.push((standard_x[is[0]] + standard_x[is[1]]) / 2.0);
                        y.push((y[is[0]] + y[is[1]]) / 2.0);
                        standard_y.push((standard_y[is[0]] + standard_y[is[1]]) / 2.0);
                    }
                }

                let indices = get_xy_sort_indices(&x, &y);
                sort_by_indices(&mut x, &indices);
                sort_by_indices(&mut standard_x, &indices);
                sort_by_indices(&mut y, &indices);
                sort_by_indices(&mut standard_y, &indices);
            }

            // Remove close points
            let mut i = 0;
            while i + 1 < x.len() {
                while i + 1 < x.len()
                    && ((standard_x[i + 1] - standard_x[i]).powi(2)
                        + ((standard_y[i + 1] - standard_y[i])
                            * (drawing_coords.chart_height / drawing_coords.chart_width))
                            .powi(2))
                    .sqrt()
                        < dash_length
                {
                    x.remove(i + 1);
                    standard_x.remove(i + 1);
                    y.remove(i + 1);
                    standard_y.remove(i + 1);
                }
                i += 1;
            }

            for is in (0..x.len())
                .collect::<Vec<_>>()
                .windows(2)
                .enumerate()
                .filter_map(|(i, is)| if i % 2 == 0 { Some(is) } else { None })
            {
                ret.push(vec![(x[is[0]], y[is[0]]), (x[is[1]], y[is[1]])]);
            }
        } else {
            ret.push(x.iter().cloned().zip(y.iter().cloned()).collect::<Vec<_>>());
        }
    }

    Ok(ret)
}

pub fn draw_cis(
    x: &[f64],
    ci_lower: &[f64],
    ci_upper: &[f64],
    drawing_coords: &DrawingCoords,
) -> Result<Vec<(f64, f64)>, Box<dyn Error>> {
    let mut x = x.to_vec();
    for x in &mut x {
        if x.is_infinite() {
            if x.is_sign_positive() {
                *x = drawing_coords.x_max + drawing_coords.x_space;
            } else {
                *x = drawing_coords.x_min - drawing_coords.x_space;
            }
        }
    }
    let mut ci_lower = ci_lower.to_vec();
    for ci in &mut ci_lower {
        if ci.is_infinite() {
            if ci.is_sign_positive() {
                *ci = drawing_coords.y_max + drawing_coords.y_space;
            } else {
                *ci = drawing_coords.y_min - drawing_coords.y_space;
            }
        }
    }
    let mut ci_upper = ci_upper.to_vec();
    for ci in &mut ci_upper {
        if ci.is_infinite() {
            if ci.is_sign_positive() {
                *ci = drawing_coords.y_max + drawing_coords.y_space;
            } else {
                *ci = drawing_coords.y_min - drawing_coords.y_space;
            }
        }
    }

    let indices = get_xy_sort_indices(&x, &ci_lower);
    sort_by_indices(&mut x, &indices);
    sort_by_indices(&mut ci_lower, &indices);
    sort_by_indices(&mut ci_upper, &indices);

    let ci_v = x
        .iter()
        .cloned()
        .zip(ci_lower.iter().cloned())
        .chain(x.iter().cloned().zip(ci_upper.iter().cloned()).rev())
        .collect::<Vec<_>>();

    Ok(ci_v)
}