//

// Created by Alan Freitas on 2020-07-03.

//

#include <fstream>

#include <cmath>

#include <string>

#include <random>

#include <array>

#include <regex>

#include <set>

#include <iostream>

#include <matplot/util/common.h>

#include <matplot/util/colors.h>

#include <CImg.h>

#ifdef _WIN32

# include <windows.h>

# define PCLOSE _pclose

# define POPEN _popen

# define FILENO _fileno

#else

# define PCLOSE pclose

# define POPEN popen

# define FILENO fileno

#endif

namespace matplot {

bool iequals(const std::string& str1, const std::string& str2) {

if (str1.size() != str2.size()) {

return false;

} else {

for (size_t i = 0; i < str1.size(); ++i) {

if (toupper(str1[i]) != toupper(str2[i])) {

return false;

}

}

}

return true;

}

bool is_true(const std::string &str) {

return iequals(str, "on") || iequals(str, "true") || iequals(str, "yes");

}

bool is_false(const std::string& str) {

return iequals(str, "off") || iequals(str, "false") || iequals(str, "no");

}

std::string run_and_get_output(const std::string &cmd) {

std::unique_ptr<FILE, decltype(&PCLOSE)> pipe(POPEN(cmd.c_str(), "r"), PCLOSE);

if (!pipe) {

throw std::runtime_error("popen() failed!");

}

std::array<char, 128> buffer{};

std::string result;

while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) {

result += buffer.data();

}

return result;

}

std::string escape(const std::string &label) {

return std::regex_replace(label, std::regex("\""), "\\\"");

}

std::vector<double> linspace(double d1, double d2, size_t n) {

std::vector<double> x(n);

for (size_t i = 0; i < x.size(); ++i) {

x[i] = d1 + static_cast<double>(i) * (d2 - d1)/static_cast<double>(n - 1);

}

return x;

}

std::vector<double> linspace(double d1, double d2) {

return linspace(d1,d2,100);

}

std::vector<double> logspace(double d1, double d2, size_t n) {

return transform(linspace(d1,d2,n), [](double x) {return pow(10,x);});

}

std::vector<double> logspace(double d1, double d2) {

return logspace(d1,d2,50);

}

std::vector<double> iota(double d1, double step, double d2) {

vector_1d r;

r.reserve((d2 - d1)/step);

for (double i = d1; i <= d2; i += step) {

r.emplace_back(i);

}

return r;

}

std::vector<double> iota(double d1, double d2) {

return iota(d1,1.,d2);

}

std::vector<double> transform(const std::vector<double>& x, std::function<double(double)> fn) {

std::vector<double> y(x.size());

for (size_t i = 0; i < y.size(); ++i) {

y[i] = fn(x[i]);

}

return y;

}

double truncate(double x, double lower_bound, double upper_bound) {

if (x <= lower_bound) {

return lower_bound;

} else if (x >= upper_bound) {

return upper_bound;

} else {

return x;

}

}

std::vector<double> transform(const std::vector<double>& x, const std::vector<double>& y, std::function<double(double,double)> fn) {

std::vector<double> z(x.size());

size_t n = std::min(x.size(), y.size());

for (size_t i = 0; i < n; ++i) {

z[i] = fn(x[i], y[i]);

}

return z;

}

vector_2d transform(const vector_2d& x, std::function<double(double)> fn) {

const size_t n_rows = x.size();

const size_t n_cols = x[0].size();

vector_2d z(n_rows, vector_1d(n_cols, 0.));

for (size_t i = 0; i < n_rows; ++i) {

for (size_t j = 0; j < n_cols; ++j) {

z[i][j] = fn(x[i][j]);

}

}

return z;

}

vector_2d transform(const vector_2d& x, const vector_2d& y, std::function<double(double,double)> fn) {

const size_t n_rows = std::min(x.size(), y.size());

const size_t n_cols = std::min(x[0].size(), y[0].size());

vector_2d z(n_rows, vector_1d(n_cols, 0.));

for (size_t i = 0; i < n_rows; ++i) {

for (size_t j = 0; j < n_cols; ++j) {

z[i][j] = fn(x[i][j], y[i][j]);

}

}

return z;

}

vector_1d flatten(const vector_2d& x) {

vector_1d x_line;

x_line.reserve(x.size() * x[0].size());

for (size_t i = 0; i < x.size(); ++i) {

x_line.insert(x_line.end(), x[i].begin(), x[i].end());

}

return x_line;

}

double max(const std::vector<double>& x) {

if (x.empty()) {

return -std::numeric_limits<double>::max();

} else {

return *std::max_element(x.begin(), x.end());

}

}

double min(const std::vector<double>& x) {

if (x.empty()) {

return std::numeric_limits<double>::max();

} else {

return *std::min_element(x.begin(), x.end());

}

}

double min(const std::vector<std::vector<double>>& x) {

double m = min(x[0]);

for (size_t i = 1; i < x.size(); ++i) {

m = min(m, min(x[i]));

}

return m;

}

double max(const std::vector<std::vector<double>>& x) {

double m = max(x[0]);

for (size_t i = 1; i < x.size(); ++i) {

m = max(m, max(x[i]));

}

return m;

}

std::pair<double,double> minmax(const std::vector<double>& x) {

if (x.empty()) {

return std::make_pair(min(x), max(x));

} else {

auto [min_it, max_it] = std::minmax_element(x.begin(), x.end());

return std::make_pair(*min_it, *max_it);

}

}

double mean(const std::vector<double>& x) {

double sum = 0;

for (const double& xi: x) {

sum += xi/static_cast<double>(x.size());

}

return sum;

}

double stddev(const std::vector<double>& x) {

double m = mean(x);

double sum = 0;

for (const double& xi: x) {

sum += pow(xi - m, 2);

}

return sqrt(sum/static_cast<double>(x.size() - 1));

}

// windows.h library defines a function-like macro with the name of max()

#ifdef max

// store the predefined macro in a new one

#define _TEMP_MACRO_ max

// undefine the problamatic macro.

#undef max

#endif

void wait() {

#ifdef MATPLOT_BUILD_FOR_DOCUMENTATION_IMAGES

return;

#endif

// http://www.cplusplus.com/forum/beginner/106769/#msg578325

// Clear failed/error states of the stream if they are set

std::cin.clear() ;

// If there are any characters in the input buffer

if( std::cin.rdbuf()->in_avail() ) {

// Throw them away

std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');

}

std::cout << "Press ENTER to continue...";

std::cin.ignore( std::numeric_limits<std::streamsize>::max(), '\n' );

}

#ifdef _Temp_MACRO_

// restore the max() macro.

#define max _TEMP_MACRO_

// undefine the temporary macro.

#undef _TEMP_MACRO_

#endif

namespace detail {

static std::mt19937& default_generator() {

static std::mt19937 g((std::random_device())());

return g;

}

}

double randn(double mean, double stddev) {

std::normal_distribution<double> d(mean,stddev);

return d(detail::default_generator());

}

std::vector<double> randn(size_t n, double mean, double stddev) {

std::vector<double> x(n);

std::generate(x.begin(), x.end(), [&]{return randn(mean,stddev);});

return x;

}

double rande(double lambda, double multiplier) {

std::exponential_distribution<double> d(lambda);

return d(detail::default_generator()) * multiplier;

}

std::vector<double> rande(size_t n, double lambda, double multiplier) {

std::vector<double> x(n);

std::generate(x.begin(), x.end(), [&]{return rande(lambda,multiplier);});

return x;

}

double randp(double scale, double shape) {

std::exponential_distribution<double> d(shape);

return scale * exp(d(detail::default_generator()));

}

std::vector<double> randp(size_t n, double scale, double shape) {

std::vector<double> x(n);

std::generate(x.begin(), x.end(), [&]{return randp(scale,shape);});

return x;

}

int randi(int imin, int imax) {

std::uniform_int_distribution<int> d(imin,imax);

return d(detail::default_generator());

}

std::vector<int> randi(size_t n, int imin, int imax) {

std::vector<int> x(n);

std::generate(x.begin(), x.end(), [&]{return randi(imin,imax);});

return x;

}

double rand(double lower_bound, double upper_bound) {

std::uniform_real_distribution<double> d(lower_bound,upper_bound);

return d(detail::default_generator());

}

std::vector<double> rand(size_t n, double lower_bound, double upper_bound) {

std::vector<double> x(n);

std::generate(x.begin(), x.end(), [&]{return rand(lower_bound,upper_bound);});

return x;

}

vector_2d rand(size_t rows, size_t cols, double lower_bound, double upper_bound) {

vector_2d x(rows);

for (size_t i = 0; i < x.size(); ++i) {

x[i] = rand(cols, lower_bound, upper_bound);

}

return x;

}

std::vector<double> unique(const std::vector<double>& x) {

std::set<double> s;

for (const auto &v: x) {

s.emplace(v);

}

std::vector<double> u(s.begin(), s.end());

return u;

}

std::vector<double> reshape(const std::vector<std::vector<double>>& X) {

std::vector<double> x;

for (const auto &vec : X) {

for (const auto &item : vec) {

x.emplace_back(item);

}

}

return x;

}

std::vector<double> concat(const std::vector<double>& a, const std::vector<double>& b) {

std::vector<double> c(a.begin(), a.end());

c.insert(c.end(), b.begin(), b.end());

return c;

}

std::string fileread(const std::string& filename) {

std::ifstream t(filename);

if (!t) {

throw std::runtime_error("Cannot open the file " + filename);

}

std::string str((std::istreambuf_iterator<char>(t)),

std::istreambuf_iterator<char>());

return str;

}

double cosd(double degrees) {

// 0 -> 0

// 180 -> pi

// 360 -> 2pi

return cos(pi*degrees/180);

}

double sind(double degrees) {

// 0 -> 0

// 180 -> pi

// 360 -> 2pi

return sin(pi*degrees/180);

}

std::vector<double> deg2rad(const std::vector<double>& degrees) {

std::vector<double> rad(degrees.size());

for (size_t i = 0; i < degrees.size(); ++i) {

rad[i] = pi * degrees[i]/180;

}

return rad;

}

double deg2rad(double rho) {

return rho * 180. / pi;

}

double round_polar_max(double polar_max) {

// we need to round the polar max because gnuplot

// does not use the max value for the r-axis.

// it uses a round value. If we explicitly set the

// max value, the x-axis and y-axis get into conflict

// and things get out of scale.

// So the only possible solution is to foresee what

// value gnuplot will set for the r-axis limit.

double decimals = floor(log10(polar_max));

double base = polar_max / pow(10,decimals);

if (base > 1. && base < 1.5) {

base = 1.5;

} else {

base = ceil(base);

if (base == 7.) {

base += 1.0;

} else if (base == 9 || base == 10) {

base = 1.0;

decimals++;

}

}

return base * pow(10,decimals);

}

double to_positive_radian(double radian) {

radian = std::fmod(radian, pi * 2.);

while (radian < 0) {

// pretty sure this comparison is valid for doubles and floats

radian += pi * 2.;

}

return radian;

}

double vector_radians(double x, double y) {

return to_positive_radian(atan2(y,x));

}

double vector_magnitude(double x, double y) {

return sqrt(pow(x,2)+pow(y,2));

}

double vector_radians(std::complex<double> v) {

return vector_radians(v.real(),v.imag());

}

double vector_magnitude(std::complex<double> v) {

return vector_magnitude(v.real(),v.imag());

}

std::pair<double,double> elliptic2cart(double theta, double rho_w, double rho_h) {

return std::make_pair(rho_w * cos(theta), rho_h * sin(theta));

}

std::pair<vector_1d,vector_1d> elliptic2cart(vector_1d theta, vector_1d rho_w, vector_1d rho_h) {

std::pair<vector_1d,vector_1d> r;

r.first.resize(theta.size());

r.second.resize(theta.size());

for (size_t i = 0; i < theta.size(); ++i) {

std::tie(r.first[i], r.second[i]) = elliptic2cart(theta[i], rho_w[i], rho_h[i]);

}

return r;

}

std::pair<vector_1d,vector_1d> elliptic2cart(vector_1d theta, double rho_w, double rho_h) {

std::pair<vector_1d,vector_1d> r;

r.first.resize(theta.size());

r.second.resize(theta.size());

for (size_t i = 0; i < theta.size(); ++i) {

std::tie(r.first[i], r.second[i]) = elliptic2cart(theta[i], rho_w, rho_h);

}

return r;

}

std::pair<double,double> pol2cart(double theta, double rho) {

return std::make_pair(rho * cos(theta), rho * sin(theta));

}

std::pair<vector_1d,vector_1d> pol2cart(vector_1d theta, vector_1d rho) {

return std::make_pair(

transform(theta, rho, [](double theta, double rho) { return rho * cos(theta); }),

transform(theta, rho, [](double theta, double rho) { return rho * sin(theta); })

);

}

std::pair<vector_1d,vector_1d> pol2cart(vector_1d theta, double rho) {

return std::make_pair(

transform(theta, [&](double theta) { return rho * cos(theta); }),

transform(theta, [&](double theta) { return rho * sin(theta); })

);

}

std::pair<vector_2d,vector_2d> meshgrid(const vector_1d& x, const vector_1d& y) {

vector_2d X(y.size(), std::vector<double>(x.size(),0));

vector_2d Y(y.size(), std::vector<double>(x.size(),0));

for (size_t i = 0; i < y.size(); ++i) {

for (size_t j = 0; j < x.size(); ++j) {

X[i][j] = x[j];

Y[i][j] = y[i];

}

}

return std::make_pair(X,Y);

}

std::pair<vector_2d,vector_2d> meshgrid(const vector_1d& x_and_y) {

return meshgrid(x_and_y, x_and_y);

}

vector_1d gradient(const vector_1d& z, double spacing) {

const size_t n = z.size();

vector_1d d(n);

d[0] = (z[1] - z[0]) / spacing;

d[n-1] = (z[n-1] - z[n-2]) / spacing;

for (size_t i = 1; i < n - 1; ++i) {

d[i] = ((z[i] - z[i-1]) + (z[i+1] - z[i])) / spacing / 2.;

}

return d;

}

std::pair<vector_2d,vector_2d> gradient(const vector_2d& z, double spacing_x, double spacing_y) {

const size_t n1 = z.size();

const size_t n2 = z[0].size();

vector_2d dx(n1, vector_1d(n2, 0.));

vector_2d dy(n1, vector_1d(n2, 0.));

for (size_t row = 0; row < n1; ++row) {

dx[row][0] = (z[row][1] - z[row][0]) / spacing_x;

dx[row][n2-1] = (z[row][n2-1] - z[row][n2-2]) / spacing_x;

for (size_t col = 1; col < n2 - 1; ++col) {

dx[row][col] = ((z[row][col] - z[row][col-1]) + (z[row][col+1] - z[row][col])) / spacing_x / 2.;

}

}

for (size_t col = 0; col < n2; ++col) {

dy[0][col] = (z[1][col] - z[0][col]) / spacing_y;

dy[n1-1][col] = (z[n1-1][col] - z[n1-2][col]) / spacing_y;

for (size_t row = 1; row < n1 - 1; ++row) {

dy[row][col] = ((z[row][col] - z[row-1][col]) + (z[row+1][col] - z[row][col])) / spacing_y / 2.;

}

}

return std::make_pair(dx,dy);

}

std::pair<vector_2d,vector_2d> gradient(const vector_2d& z, double spacing) {

return gradient(z,spacing,spacing);

}

double peaks(double x, double y) {

return 3. * pow(1. - x, 2.) * exp(-pow(x, 2.) - pow(y + 1., 2.))

- 10. * (x / 5. - pow(x, 3.) - pow(y, 5.))

* exp(-pow(x, 2.) - pow(y, 2.))

- 1. / 3. * exp(-pow(x + 1., 2.) - pow(y, 2.));

}

vector_2d peaks(const vector_2d& X, const vector_2d& Y) {

const size_t n_rows = std::min(X.size(), Y.size());

const size_t n_cols = std::min(X[0].size(), Y[0].size());

vector_2d Z(n_rows, vector_1d(n_cols, 0.));

for (size_t i = 0; i < n_rows; ++i) {

for (size_t j = 0; j < n_cols; ++j) {

Z[i][j] = peaks(X[i][j],Y[i][j]);

}

}

return Z;

}

std::tuple<vector_2d,vector_2d,vector_2d> peaks(size_t N) {

std::vector<double> x = linspace(-3.,+3., N);

std::vector<double> y = linspace(-3.,+3., N);

auto [X,Y] = meshgrid(x,y);

return std::make_tuple(X,Y,peaks(X,Y));

}

vector_2d zeros(size_t rows, size_t cols) {

return vector_2d(rows, vector_1d(cols, 0.));

}

vector_2d ones(size_t rows, size_t cols) {

return vector_2d(rows, vector_1d(cols, 0.));

}

image_channels_t cimg2channels(const cimg_library::CImg<unsigned char>& cimg_image) {

const size_t n_channels = cimg_image.spectrum();

const size_t height = cimg_image.height();

const size_t width = cimg_image.width();

image_channels_t img(n_channels, image_channel_t(height, image_row_t(width, 0)));

for (size_t channel = 0; channel < n_channels; ++channel) {

for (size_t i = 0; i < cimg_image.height(); ++i) {

for (size_t j = 0; j < cimg_image.width(); ++j) {

img[channel][i][j] = cimg_image.operator()(j,i,0,channel);

}

}

}

return img;

}

cimg_library::CImg<unsigned char> channels2cimg(const image_channels_t& A) {

cimg_library::CImg<unsigned char> cimg_image(A[0].size(),A[0][0].size(),1,A.size());

const size_t n_channels = A.size();

const size_t height = A[0].size();

const size_t width = A[0][0].size();

for (size_t channel = 0; channel < n_channels; ++channel) {

for (size_t i = 0; i < cimg_image.height(); ++i) {

for (size_t j = 0; j < cimg_image.width(); ++j) {

cimg_image.operator()(j,i,0,channel) = A[channel][i][j];

}

}

}

return cimg_image;

}

image_channels_t imread(const std::string& filename) {

try {

cimg_library::CImg<unsigned char> cimg_image(filename.c_str());

return cimg2channels(cimg_image);

} catch (...) {

// return empty tuple if we can't open the file

return image_channels_t{0};

}

}

image_channels_t imresize(const image_channels_t& A, size_t height, size_t width, image_interpolation m) {

auto cimg_image = channels2cimg(A);

int interpolation_type = 0;

switch (m) {

case image_interpolation::raw:

interpolation_type = -1;

break;

case image_interpolation::additional_space:

interpolation_type = 0;

break;

case image_interpolation::nearest:

interpolation_type = 1;

break;

case image_interpolation::moving_average:

interpolation_type = 2;

break;

case image_interpolation::bilinear:

interpolation_type = 3;

break;

case image_interpolation::grid:

interpolation_type = 4;

break;

case image_interpolation::bicubic:

interpolation_type = 5;

break;

case image_interpolation::lanczos:

interpolation_type = 6;

break;

}

cimg_image.resize(width, height, 1, A.size(), interpolation_type);

return cimg2channels(cimg_image);

}

image_channels_t imresize(const image_channels_t& A, double scale, image_interpolation m) {

if (A.empty()) {

return image_channels_t{};

}

auto [h,w] = size(A[0]);

return imresize(A, h * scale, w * scale, m);

}

void imwrite(const image_channels_t& A, const std::string& filename) {

auto cimg_image = channels2cimg(A);

cimg_image.save(filename.c_str());

}

image_channel_t rgb2gray(const image_channels_t& A) {

image_channel_t img(A[0].size(), image_row_t(A[0][0].size()));

for (size_t i = 0; i < A.size(); ++i) {

for (size_t j = 0; j < A[i].size(); ++j) {

img[i][j] = (A[0][i][j]+A[1][i][j]+A[2][i][j])/3;

}

}

return img;

}

image_channels_t gray2rgb(const image_channel_t& A, const std::vector<std::vector<double>>& colormap) {

image_channels_t img(3, image_channel_t(A.size(), image_row_t(A[0].size())));

for (size_t i = 0; i < A.size(); ++i) {

for (size_t j = 0; j < A[i].size(); ++j) {

color_array c = colormap_interpolation(A[i][j],0,255,colormap);

img[0][i][j] = round(c[1] * 255);

img[1][i][j] = round(c[2] * 255);

img[2][i][j] = round(c[3] * 255);

}

}

return img;

}

image_channels_t gray2rgb(const image_channel_t& A) {

static std::vector<std::vector<double>> map = {{0, 0, 0}, {1, 1, 1}};

return gray2rgb(A, map);

}

image_channels_t gray2rgb(const image_channels_t& A, const std::vector<std::vector<double>>& colormap) {

return gray2rgb(A[0],colormap);

}

image_channels_t gray2rgb(const image_channels_t& A) {

return gray2rgb(A[0]);

}

void imwrite(const image_channel_t& A, const std::vector<std::vector<double>>& colormap, const std::string& filename) {

imwrite(gray2rgb(A,colormap), filename);

}

image_channels_t imvignette(const image_channels_t& A, double min_radius, double exponent) {

auto image = A;

if (image.empty()) {

return A;

}

if (image.size() < 2) {

image = gray2rgb(image);

}

if (image.size() < 3) {

image.emplace_back(image[0]);

}

auto [h,w] = size(image[0]);

if (image.size() < 4) {

image.emplace_back(image_channel_t(h, image_row_t(w, 255)));

}

double center_x = w/2;

double center_y = h/2;

double radius_sq = pow(min(h/2,w/2) * min_radius, 2);

double max_t1 = pow(static_cast<double>(h) - center_y, 2);

double max_t2 = pow(static_cast<double>(w) - center_x,2);

double max_t_minus_radius = (max_t1 + max_t2) - radius_sq;

auto& alpha_channel = image[3];

for (int i = 0; i < h; ++i) {

for (int j = 0; j < w; ++j) {

double t1 = pow(static_cast<double>(i) - center_y, 2);

double t2 = pow(static_cast<double>(j) - center_x,2);

if (t1 + t2 > radius_sq) {

double norm_dist_from_r = ((t1 + t2) - radius_sq) / max_t_minus_radius;

norm_dist_from_r = pow(norm_dist_from_r, exponent);

alpha_channel[i][j] = 255 * (1. - norm_dist_from_r);

}

}

}

return image;

}

vector_2d transpose(const vector_2d& z) {

vector_2d z2(z[0].size(), vector_1d(z.size()));

for (size_t i = 0; i < z.size(); ++i) {

for (size_t j = 0; j < z[0].size(); ++j) {

z2[j][i] = z[i][j];

}

}

return z2;

}

std::vector<std::string> tokenize(const std::string &text, std::string delimiters) {

std::vector<std::string> tokens;

size_t pos = 0;

while ((pos = text.find_first_not_of(delimiters, pos)) != std::string::npos) {

size_t end_pos = text.find_first_of(delimiters, pos);

tokens.emplace_back(text.substr(pos, end_pos - pos));

pos = end_pos;

}

return tokens;

}

std::pair<std::vector<std::string>,std::vector<size_t>> wordcount(const std::vector<std::string> &tokens, const std::vector<std::string> &black_list, const std::string &delimiters, size_t max_cloud_size) {

// count the frequency of each token

const bool bl_sorted = std::is_sorted(black_list.begin(), black_list.end());

std::map<std::string, size_t> word_frequency;

for (const auto &token : tokens) {

const bool is_in_black_list = bl_sorted ?

std::binary_search(black_list.begin(), black_list.end(), token) :

std::find(black_list.begin(), black_list.end(), token) != black_list.end();

if (!is_in_black_list) {

auto it = word_frequency.find(token);

if (it != word_frequency.end()) {

++word_frequency[token];

} else {

word_frequency[token] = 1;

}

}

}

// sort tokens by frequency

std::multimap<size_t, std::string, std::greater<>> frequency_word;

for (const auto &[token, count] : word_frequency) {

frequency_word.emplace(count, token);

}

// get the max_cloud_size most frequent tokens

std::vector<std::string> unique_tokens;

std::vector<size_t> token_count;

size_t i = 0;

for (const auto &[count, token] : frequency_word) {

unique_tokens.emplace_back(token);

token_count.emplace_back(count);

++i;

if (i > max_cloud_size) {

break;

}

}

// return tokens and counts

return std::make_pair(unique_tokens, token_count);

}

std::pair<std::vector<std::string>,std::vector<size_t>> wordcount(const std::string &text, const std::vector<std::string> &black_list, const std::string &delimiters, size_t max_cloud_size) {

auto tokens = tokenize(text);

return wordcount(tokens, black_list, delimiters, max_cloud_size);

}

ticks_results calcticks(double limits_min, double limits_max, bool horizontal, double text_size, bool separateExp, bool log) {

// Default limit on label length (in characters)

constexpr size_t defMaxChars = 9;

// Minimum value of the (upper) limit on label length (characters). Setting

// this value too small will cause problems.

// Note: This value does not affect the minimum length of the ticklabels.

constexpr size_t minChars = 6;

// Default font size for exponents (assumes that font units are 'points')

constexpr size_t defExpFontSize = 7;

// Upper limit on number of ticks returned by CALCTICKS

constexpr size_t initMaxTicks = 11;

// Multiplier for textSize for vertical orientation when scale is 'log', to

// account for labels using exponential notation.

constexpr double vertExpScale = 1.3;

// Maximum number of characters in label string

// Determines numerical precision displayed by the labels,

// and also affects the tick spacing for horizontal orientation.

size_t maxChars = defMaxChars;

if (maxChars < minChars) {

maxChars = minChars;

}

// Initial calculations

if (log) {

std::invalid_argument("Not implemented yet. The library does not need that as it is.");

}

// Data range

double range = limits_max - limits_min;

// Get eps values for rounding

double lEps = eps(limits_min);

double uEps = eps(limits_max);

double minEps = min(lEps,uEps);

// Vector of allowed tick counts

std::vector<double> testTickCounts = iota(2,initMaxTicks);

// Make a list of rough intervals as a starting point

std::vector<double> roughInts = transform(testTickCounts, [&](double x) { return range/(x-1); });

// Vector of 'nice' intervals

std::vector<double> niceVec = {1,2,5,10};

// Find nice intervals

// Normalize rough intervals by their scale

std::vector<double> decRoughInts = transform(roughInts, [](double x) { return floor(log10(x)); });

std::vector<double> normRoughInts = transform(roughInts, decRoughInts, [](double x, double y) { return x/pow(10.,y); });

// Get the distances to nice intervals, pick the shortest

std::vector<std::vector<double>> deltas(niceVec.size(), std::vector<double>(roughInts.size()));

for (size_t i = 0; i < deltas.size(); ++i) {

for (size_t j = 0; j < deltas[i].size(); ++j) {

deltas[i][j] = abs(normRoughInts[j] - niceVec[i]);

}

}

// idx of ticks with min distance

std::vector<size_t> idx(roughInts.size(), 0);

for (size_t i = 0; i < deltas.size(); ++i) {

for (size_t j = 0; j < deltas[i].size(); ++j) {

if (deltas[i][j] < deltas[idx[j]][j]) {

idx[j] = i;

}

}

}

// Get the nice intervals and scores

std::vector<double> niceInts(roughInts.size());

for (size_t i = 0; i < niceInts.size(); ++i) {

niceInts[i] = niceVec[idx[i]] * pow(10,decRoughInts[i]);

}

// Remove duplicates

niceInts = unique(niceInts);

// Get upper and lower limits, fixed by the list of nice intervals. Round

// out to make sure we get ticks at the original limits.

std::vector<double> lLims = transform(niceInts, [&](double x) { return floor(limits_min/x)*x; });

std::vector<double> uLims = transform(niceInts, [&](double x) { return ceil(limits_max/x)*x; });

// Get tick counts using the list of nice intervals and limits

std::vector<double> nTicks(niceInts.size());

for (size_t i = 0; i < nTicks.size(); ++i) {

nTicks[i] = floor(1 + (uLims[i] - lLims[i] + 10*minEps)/niceInts[i]);

}

// Shrink nice limits that are outside of original limits

for (size_t idx = 0; idx < lLims.size(); ++idx) {

if (lLims[idx] < limits_min - 10*eps(limits_min)) {

nTicks[idx] = nTicks[idx]-1;

lLims[idx] = lLims[idx] + niceInts[idx];

}

if (uLims[idx] > limits_max + 10*eps(limits_min)) {

nTicks[idx] = nTicks[idx] - 1;

uLims[idx] = uLims[idx] - niceInts[idx];

}

}

// Set values that are almost exactly the original limits to be the original

// limit value.

for (size_t idx = 0; idx < lLims.size(); ++idx) {

if (abs(lLims[idx] - limits_min) < 10*eps(limits_min)) {

lLims[idx] = limits_min;

}

if (abs(uLims[idx] - limits_max) < 10*eps(limits_max)) {

uLims[idx] = limits_max;

}

}

// Discard values where the limits are reversed or equal

size_t lim_idx = 0;

while (lim_idx < lLims.size()) {

if (lLims[lim_idx] >= uLims[lim_idx]) {

lLims.erase(lLims.begin() + lim_idx);

uLims.erase(uLims.begin() + lim_idx);

nTicks.erase(nTicks.begin() + lim_idx);

niceInts.erase(niceInts.begin() + lim_idx);

} else {

++lim_idx;

}

}

// Determine label size for each interval

// Get the decade span of the limits and the decade of the intervals

std::vector<double> maxAbs = transform(lLims,uLims,[](double x, double y) {

return abs(x) > abs(y) ? abs(x) : abs(y);

});

std::vector<double> decMax = transform(maxAbs, [](double x) { return floor(log10(x)); });// - nDec;

std::vector<double> decInts = transform(niceInts, [](double x) { return floor(log10(x)); }); // - nDec;

// Get the number of characters needed for tick labels for normal notation

std::vector<double> labelChars = transform(decInts, decMax, [&](double x, double y) {

return max(y + 1, 1.) + ((x < 0) ? (1 - x) : 0.);

});

for (size_t i = 0; i < labelChars.size(); ++i) {

if (labelChars[i] > maxChars - 1) {

labelChars[i] = maxChars - 1;

}

}

// Handle exponential notation

// Determine whether or not to use exponential notation

std::vector<bool> isExp(decInts.size());

if (separateExp) {

for (size_t i = 0; i < decInts.size(); ++i) {

// Large numbers:

isExp[i] = (decMax[i] > 6) || (decMax[i] == 6 && decInts[i] > 0);

// Small numbers:

isExp[i] = isExp[i] || (decMax[i] < -3) || (decMax[i] == -3 && decInts[i] < -5);

}

} else {

for (size_t i = 0; i < decInts.size(); ++i) {

// Large numbers:

isExp[i] = decMax[i] > 6 || (decMax[i] == 6 && decInts[i] > 3);

// Small numbers:

isExp[i] = isExp[i] || (decMax[i] < -3) || (decMax[i] == -3 && decInts[i] < -5);

}

}

// Get length of exponential labels depending on style

std::vector<double> expChars;

if (separateExp) {

expChars = transform(decInts, decMax, [&](double x, double y) {

return 2 + max(0., min(static_cast<double>(maxChars) - 3, y - x));

});

} else {

expChars = transform(decInts, decMax, [&](double x, double y) {

return 2 + max(0., min(max(0., static_cast<double>(maxChars) - 7), y - x)) + 4;

});

}

// Select between normal and exponential label lengths

for (size_t i = 0; i < labelChars.size(); ++i) {

if (isExp[i]) {

labelChars[i] = expChars[i];

}

}

// For consistency, always include space for a negative sign, regardless of

// the sign of the actual limits

labelChars = transform(labelChars, [](double x) { return x + 1; });

// Uncomment to not include the negative sign space if it isn't needed

// labelChars = labelChars + (lLims < 0);

// Get label size based on textSize, orientation and length of label string

std::vector<double> labelSize;

if (horizontal) {

labelSize = transform(labelChars, [&](double labelChars) { return text_size * labelChars; });

} else {

labelSize = transform(labelChars, [&](double labelChars) { return text_size; });

}

// Choose the best interval

// Maximum number of ticks without overlapping labels

std::vector<double> nMax(uLims.size());

for (size_t i = 0; i < nMax.size(); ++i) {

nMax[i] = floor((uLims[i]-lLims[i]+10*minEps)/labelSize[i]) + 1;

}

// Modify this value based on initMaxTicks

std::vector<double> nMaxScore = transform(nMax, [&](double x) { return min(x,initMaxTicks + 0.25*(x-initMaxTicks)); });

// Calculate a score based on the number of ticks relative to the maximum.

std::vector<double> nTickScore = transform(nTicks, nMaxScore, [](double x, double y) { return 1-pow((x/y - 0.7),2); });

// Severe penalty for more than nMax ticks

std::vector<double> penalty = transform(nTicks, nMax, [](double x, double y) { return 4*x/y; });

for (size_t i = 0; i < nTicks.size(); ++i) {

if (nTicks[i]/nMax[i] <= 1.) {

penalty[i] = 0;

}

}

// Penalty for more than 0.75*nMax but less than nMax

std::vector<double> penaltyScale = transform(nTicks, nMax, [](double x, double y) { return x/y; });

for (size_t i = 0; i < penaltyScale.size(); ++i) {

if (penaltyScale[i] > 0.75 && nTicks[i]/nMax[i] <= 1.) {