axis.cpp

//**************************************************************************************************
//*  FRIDA: fast reliable interactive data analysis
//*  (C) Joachim Wuttke 1990-, v2(C++) 2001-
//*  http://apps.jcns.fz-juelich.de/frida
//**************************************************************************************************

//! \file  axis.cpp
//! \brief CAxis: One axis of a coordinate frame.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

#include "../readplus/macro.hpp"
#include "../trivia/math.hpp"
#include "../trivia/string_convs.hpp"

#include "axis.hpp"

using std::string;
using std::vector;
using std::cout;

#define S(a) triv::strg((a))


//! Sets infinite limits so that next plot call will recompute them.

void CAxis::set_auto()
{
    inf = -INFINITY;
    sup = +INFINITY;
}


//! Sets limits to arguments.

void CAxis::set_limits(double _inf, double _sup)
{
    if (_inf >= _sup)
        throw "Invalid limits " + S(_inf) + " .. " + S(_sup) + " for axis " + name;
    if (logflag && _inf <= 0)
        throw "Invalid nonpositive limit " + S(_inf) + " for axis " + name;
    inf = _inf;
    sup = _sup;
}


//! Sets log flag to argument, resets limits to infinite defaults.

void CAxis::set_log(bool _log)
{
    logflag = _log;
    set_auto();
}


//! Asks for new limits.

void CAxis::ask_and_set(const string& quest)
{
    for (;;) {
        double inf_in, sup_in;
        string resp1 = NMacro::readwd(quest + " [" + str() + "] ? ");
        if (resp1 == "\\q") {
            throw S("user escape to main menu");
        } else if (resp1 == "") {
            return; // don't change current value
        } else if (resp1 == "*") {
            set_auto();
            return;
        } else if (resp1 == "\\EOL") {
            inf_in = inf;
        } else if (triv::any2dbl(resp1, &inf_in)) {
            // ok
        } else {
            cout << "required input: plot axis range\n"
                    "examples:\n"
                    "    0.1 1000    two floating-point numbers\n"
                    "    \\   1000    let lower bond unchanged\n"
                    "    0.1  \\      let upper bound unchanged\n"
                    "    *           redetermine bounds from data\n"
                    "just press ENTER to let present range unchanged ["
                 << str() << "]\n";
            continue;
        }
        if (logflag && inf_in <= 0) {
            cout << "log axis requires lower bound>0\n";
            continue;
        }
        string resp2 = NMacro::readwd(".. and upper limit [" + S(sup) + "] ? ");
        if (resp2 == "\\q") {
            throw S("user escape to main menu");
        } else if (resp2 == "" || resp2 == "\\EOL") {
            sup_in = sup;
        } else if (triv::any2dbl(resp2, &sup_in)) {
            // ok
        } else {
            cout << "invalid upper bound; use '?' for help\n";
            continue;
        }
        if (sup_in <= inf_in) {
            cout << "invalid input, does not fulfill lower<upper\n";
            continue;
        }
        if (logflag && (sup_in > 1e199 || inf_in < 1e-199)) {
            cout << "log range currently restricted to 1E-199..1E199\n";
            continue;
        }
        if (!logflag && (sup_in > 3e38 || inf_in < -3e38)) {
            cout << "lin range currently restricted to -3e38..3e38\n";
            continue;
        }
        inf = inf_in;
        sup = sup_in;
        return;
    }
}


//! Sets limits to rounded values of arguments.

void CAxis::set_rounded_limits(double _inf, double _sup)
{
    static double relmargin = 0.05;
    static double digits = 2;

    if (_inf == -INFINITY || _sup == +INFINITY)
        throw "BUG: rounding " + name + " limits: unexpected inf limits";
    if (_inf > _sup)
        throw "BUG: rounding " + name + " limits: invalid plot range " + S(_inf) + " > " + S(_sup);
    if (logflag && _inf <= 0)
        throw "BUG: rounding " + name + " limits: nonpositive log inf " + S(_inf);

    inf = _inf;
    sup = _sup;

    if (inf == sup) {
        if (logflag) {
            inf /= 10;
            sup *= 10;
        } else if (inf > 0) {
            inf *= 0.9;
            sup *= 1.1;
        } else if (sup < 0) {
            inf *= 1.1;
            sup *= 0.9;
        } else {
            inf = -1;
            sup = +1;
        }
        return;
    }

    double inf2, sup2;
    double mydigits = digits;
    if (!logflag) { // lin limits
        double margin = relmargin * (sup - inf);
        do {
            sup2 = triv::round_decimal(sup + margin, mydigits);
            inf2 = triv::round_decimal(inf - margin, mydigits);
            if (sup < 0 && sup2 > 0)
                sup2 = 0;
            if (inf2 < 0 && inf > 0)
                inf2 = 0;
            mydigits += 0.1;
            // printf("DEB [%g,%g] -> [%g,%g] (%g)\n", inf, sup, inf2, sup2, mydigits);
            if (mydigits > 17) {
                std::cerr << std::setprecision(17)
                          << "PLEASE REPORT THIS BUG: escaping from endless loop in Axis::"
                             "set_rounded_limits("
                          << inf << "," << sup << ")\n";
                inf2 = inf - margin;
                sup2 = sup - margin;
                break;
            }
        } while (inf2 > inf || sup2 < sup
                 || (sup2 - inf2) > (1.33 * (sup - inf) + 1e-14 * (sup + inf)));
    } else { // log limits
        double ratio = sup / inf;
        double margin = exp(relmargin * log(ratio));
        do {
            sup2 = triv::round_decimal(sup * margin, mydigits);
            inf2 = triv::round_decimal(inf / margin, mydigits);
            mydigits += 0.5;
        } while (!((inf2 < inf) && (sup < sup2)) && mydigits < 18);
    }
    set_limits(inf2, sup2);
}


//! Are both limits finite ?

bool CAxis::finite() const { return inf != -INFINITY && sup != +INFINITY; }


//! Is argument val contained in present plot range?

bool CAxis::contains(double val) const { return inf <= val && val <= sup; }


//! Do arguments v1 and v2 agree within tolerance tol, relative to plot range?

bool CAxis::close_enough(double v1, double v2, double tol) const
{
    if (!finite())
        throw S("undefined plot range");
    if (logflag) {
        if (v1 <= 0 || v2 <= 0)
            return false;
        return fabs(log(v1 / v2)) <= tol * log(sup / inf);
    } else {
        return fabs(v1 - v2) <= tol * (sup - inf);
    }
}


//! Returns plot coordinate (between 0 and 1) of argument.

double CAxis::value2plotcoord(double v) const
{
    if (!finite())
        throw S("undefined plot range");
    if (logflag) {
        if (inf < 0 || v < 0)
            throw S("negative value in log range");
        return log(v / inf) / log(sup / inf);
    } else {
        return (v - inf) / (sup - inf);
    }
}


//! Returns half length of error bar (in plot coordinates 0..1).

double CAxis::value2ploterror(double v, double dv) const
{
    if (!finite())
        throw S("undefined plot range");
    if (dv < 0)
        throw S("negative error");
    if (dv == 0)
        return 0;
    if (logflag) {
        if (inf < 0 || v < 0)
            throw S("negative value in log range");
        return dv / v / log(sup / inf);
    } else {
        return dv / (sup - inf);
    }
}


//! Sets vector x such that its plot coordinates form an equidistent grid.

void CAxis::set_xgrid(vector<double>& x, int n) const
{
    if (!finite())
        throw S("undefined plot range");
    if (n < 2)
        throw S("call to set_xgrid with less than two points");
    if (inf >= sup)
        throw S("call to set_xgrid with invalid limits");
    x.resize(n);
    x[0] = inf;
    x[n - 1] = sup;
    for (int i = 1; i < n - 1; ++i) {
        if (logflag) {
            double c = ((double)(i)) / ((double)(n - 1));
            x[i] = pow(inf, 1 - c) * pow(sup, c);
        } else {
            x[i] = ((n - 1 - i) * inf + i * sup) / (n - 1);
        }
    }
}


//! Returns string representation of plot range.

string CAxis::str() const
{
    if (inf == -INFINITY || sup == +INFINITY)
        return "*";
    else
        return S(inf) + " " + S(sup);
}


//! Returns summary (log flag, plot limits) for use in table of plot frames.

string CAxis::info() const
{
    string ret;
    ret = logflag ? "log" : "lin";
    ret += " " + S(inf);
    ret += " " + S(sup);
    return ret;
}


//! Returns argument converted to my PostScript units (0..10).

double CAxis::pc(double v) const { return 10 * value2plotcoord(v); }


//! Returns half length of error bar (in my Postscript coordinates 0..10).

double CAxis::pcerr(double v, double dv) const { return 10 * value2ploterror(v, dv); }


//! Sets Tacks, ntpt, nt4t, *ticklim to a nice division of the plot range.

void CAxis::calc_ticks(vector<double>& Tacks, int& ntpt, int& nt4t, double* ticklim) const
{
    if (logflag)
        calc_logticks(Tacks, ntpt, nt4t, ticklim);
    else
        calc_linticks(Tacks, ntpt, nt4t, ticklim);
}


//! Sets Tacks, ntpt, nt4t, *ticklim to a nice division of the plot range.

void CAxis::calc_linticks(vector<double>& Tacks, int& ntpt, int& nt4t, double* ticklim) const
{
    double dtack;
    calc_lintacks(Tacks, ntpt, dtack);
    ticklim[0] = Tacks.front() - dtack;
    ticklim[1] = Tacks.back() + dtack;
    nt4t = Tacks.size() + 2;
}

//! Sets Tacks, *ntpt, and dtack

void CAxis::calc_lintacks(vector<double>& Tacks, int& ntpt, double& dtack) const
{
    double range = sup-inf;
    if (range <= 0)
        throw "BUG detected by calc_linticks: inf=" + S(inf) + " >= sup=" + S(sup);
    double r = log10(range);
    int ir = (int)((r < 0) ? r - 1 : r);
    double rd = r - ir; // fractional part of r
    double rdr = pow(10., rd);
    if (rdr > 9.99999) {
        dtack = 2.5 * pow(10., ir); // spacing between tacks
        ntpt = 5; // ticks per tack
    } else if (rdr > 5.00001) {
        dtack = 2 * pow(10., ir);
        ntpt = 4;
    } else if (rdr > 2.500005) {
        dtack = 1 * pow(10., ir);
        ntpt = 5;
    } else if (rdr > 1.600003) {
        dtack = .5 * pow(10., ir);
        ntpt = 5;
    } else if (rdr > 1.250002) {
        dtack = .4 * pow(10., ir);
        ntpt = 4;
    } else {
        dtack = .25 * pow(10., ir);
        ntpt = 5;
    }
    Tacks.clear();
    printf( "DEBUG (%g,%g) => (%g,%g) step %g\n",
            inf, sup, dtack * (floor(inf/dtack)-1), dtack * (ceil(sup/dtack)+1), dtack);
    for (double d = dtack * (floor(inf/dtack)-1); d <= dtack * (ceil(sup/dtack)+1); d += dtack)
        if (inf-1e-4*range <= d && d <= sup+1e-4*range)
            Tacks.push_back(d);
}


//! Sets Tacks, ntpt, nt4t, *ticklim to a nice division of the plot range.

void CAxis::calc_logticks(vector<double>& Tacks, int& ntpt, int& nt4t, double* ticklim) const
{
    if (inf >= sup)
        throw "BUG detected by calc_logticks: inf=" + S(inf) + " >= sup=" + S(sup);
    if (inf <= 0)
        throw S("BUG detected by calc_logticks: negative log argument");

    static double eins = 1 + 1e-6 - 1e-12;

    double rlgmin = log10(inf);
    double rlgmax = log10(sup);
    double rlgrel = rlgmax - rlgmin;
    if (!std::isfinite(rlgrel) || rlgrel > 100)
        throw "Excessive log range " + S(inf) + " .. " + S(sup);

    int incr; // -> decades per tack
    if (rlgrel <= 10) {
        incr = 1;
    } else {
        incr = (int)(rlgrel / 7 + 1);
        if (incr > 3)
            incr = 3 * (incr / 3); // preferred values 3, 6, ...
    }

    // first and last tack
    int minexp = (int)(rlgmin > 0 ? rlgmin + 1e-6 : rlgmin - 1e-6 - 1);
    int maxexp = (int)(rlgmax > 0 ? rlgmax + 1e-6 : rlgmax - 1e-6 - 1);
    if (incr > 1) {
        minexp -= minexp % incr;
        maxexp -= maxexp % incr;
    }

    // set Tacks
    double d1 = inf / eins;
    double d2 = sup * eins;
    Tacks.clear();
    for (int i = minexp; i <= maxexp; i += incr) {
        double r0 = pow(10., i);
        if (d1 <= r0 && r0 <= d2)
            Tacks.push_back(r0);
    }
    nt4t = Tacks.size() + 2;

    // set tick limits (may exceed actual tick range)
    int ntack = Tacks.size();
    if (ntack >= 1) {
        ticklim[0] = Tacks.front() / pow(10., incr);
        ticklim[1] = Tacks.back()  * pow(10., incr);
    } else {
        // untested. just to allow plots when no tacks in window
        ticklim[0] = pow(10., minexp);
        ticklim[1] = pow(10., maxexp + 1);
    }

    // determine #ticks / tack
    if (incr == 1 && rlgrel <= 3) {
        ntpt = 9;
    } else if (incr == 1 && rlgrel <= 7) {
        ntpt = 3; // 1-2-5-10 - Schritte
    } else if (incr == 1) {
        ntpt = 1;
    } else {
        ntpt = -incr; // this one has special semantics: it means 'incr' ticks per tack
    }

    // for narrow ranges, overwrite tacks
    if (ntack<2) {
        double dtack;
        int ntpt_dummy;
        calc_lintacks(Tacks, ntpt_dummy, dtack);
    }
}