new function plotter respects prange ... but still does not bisect at

limits.

pub/itest/curveplot.f2t

@@ -67,4 +67,19 @@ op0 26
 gya
 p
 gfa /tmp/fridatestplot.ps
+cca p0+p1*t
+crp .2<t && t<.730812
+gx 0 1
+gya
+p
+gfa /tmp/fridatestplot.ps
+gy -3 0
+p
+gfa /tmp/fridatestplot.ps
+cca sin(t)
+crp abs(f(t))<.3
+gx 0 40
+gy -1.1 .9
+p
+gfa /tmp/fridatestplot.ps
 qui

pub/lib/commands.cpp

@@ -152,7 +152,9 @@ bool frida_command(string cmd)
                 "  cm   modify function\n"
                 "  cd   set data file to fit\n"
                 "  cv   set file to convolute with (cv-: no convolution)\n"
-                "  cr   set fit range restriction\n"
+                "  cr   set range restrictions for fit and plot\n"
+                "  crf  - for fit only\n"
+                "  crp  - for plot only\n"
                 "  cf   fit\n"
                 "  cfs  fit, allow slow convolution\n"
                 "  cp   print parameters\n"
@@ -187,8 +189,8 @@ bool frida_command(string cmd)
     } else if (cmd == "cd" || cmd == "cv" || cmd == "cv-") {
         NCurveFile::set_file_reference(cmd.substr(1));
 
-    } else if (cmd == "cr") {
-        NCurveFile::change_range();
+    } else if (cmd == "cr" || cmd == "crf" || cmd == "crp") {
+        NCurveFile::change_range(cmd.substr(2));
 
     } else if (cmd == "cp") {
         NCurveFile::show_pars();

pub/lib/curve.cpp

@@ -205,20 +205,33 @@ void NCurveFile::set_par_attr(char attr)
 
 //! Queries and sets fit range restrictions.
 
-void NCurveFile::change_range()
+void NCurveFile::change_range(const string& subcmd)
 {
     FileIterator fiter(SFSel::instance()->selC());
 
-    string expr = sask("Expression restricting the fit range");
+    string which = "";
+    if (subcmd=="")
+        which = "plot and fit";
+    else if (subcmd=="f")
+        which = "fit";
+    else if (subcmd=="p")
+        which = "plot";
+    else
+        throw "BUG: invalid command";
+
+    string query = "Expression restricting the " + which + " range";
+    string expr = sask(query); // TODO: default if only one file
 
     while (COlc* fc = fiter.nextC()) {
-        fc->range_expr = expr;
-        fc->log_action("cr " + expr);
-        if (expr == "") {
-            fc->range_T = nullptr;
-            continue;
+        RNode T = expr!="" ? user_xaxparse(expr.c_str()) : nullptr;
+        fc->log_action("cr" + subcmd + " " + expr);
+        if (subcmd=="" || subcmd=="f") {
+            fc->frange_expr = expr;
+            fc->frange_T = T;
+        } else if (subcmd=="" || subcmd=="p") {
+            fc->prange_expr = expr;
+            fc->prange_T = T;
         }
-        fc->range_T = user_xaxparse(expr.c_str());
     }
 }
 

pub/lib/curve.hpp

@@ -20,7 +20,7 @@ namespace NCurveFile
 void create_fitcurve();
 void create_freecurve();
 void change_expr();
-void change_range();
+void change_range(const string& subcmd);
 void set_file_reference(const string& which);
 
 void set_properties(string which);

pub/lib/fit.cpp

@@ -522,7 +522,7 @@ void NCurveFit::fit(bool _allow_slow_conv)
         if ((fc->weighing == COlc::_VAR || fc->weighing == COlc::_VARD) && !fd->VS(0)->has_dy())
             throw S("Data have no error bars, cannot apply chosen fit weights");
 
-        if (fc->range_T) {
+        if (fc->frange_T) {
             throw S("TODO: restore range evaluation");
             /*
             // overwrite fd with data that obye the range restriction:
@@ -532,7 +532,7 @@ void NCurveFit::fit(bool _allow_slow_conv)
                 PSpec sout = PSpec( new CSpec );
                 sout->z = sin->z;
                 vector<double> range( sin->size() );
-                fc->range_T->tree_vec_val( &range, 0, fc->kd, j );
+                fc->frange_T->tree_vec_val( &range, nullptr, fc->kd, j );
                 bool with_dy = sin->has_dy();
                 for ( int i=0; i<sin->size(); i++ ) {
                     if ( range[i] ) {

pub/lib/obj.hpp

@@ -172,7 +172,7 @@ class CObjVecInt : public CObjVecNum
 {
 public:
     vector<int> v; //!< The data.
-    CObjVecInt(int n = 0) : CObjVecNum(), v(n, -1) {}
+    CObjVecInt(int n=0, int val=0) : CObjVecNum(), v(n, val) {}
     CObjVecInt(const vector<int>& _v) : CObjVecNum(), v(_v) {}
     inline int size() const { return v.size(); }
     // PObj clone() const { return PObjVecInt( new CObjVecInt( *this ) ); }

pub/lib/olf.cpp

@@ -127,8 +127,10 @@ void COlc::copy_meta_C_from_other(const COlc* other)
     kd = other->kd;
     weighing = other->weighing;
     plot_to_grid = other->plot_to_grid;
-    range_expr = other->range_expr;
-    range_T = other->range_T;
+    frange_expr = other->frange_expr;
+    frange_T = other->frange_T;
+    prange_expr = other->prange_expr;
+    prange_T = other->prange_T;
     chi2 = other->chi2;
 }
 
@@ -637,14 +639,35 @@ RObjVec COlc::eval_curve(const vector<double>& vt, int k, int j, bool want_error
 
 RObj COlc::eval_curve_scalar(double arg, int k, int j, bool want_error) const
 {
-    vector<double> vt(1, arg);
-    RObjVec cu = eval_curve(vt, k, j, want_error);
+    RObjVec cu = eval_curve({arg}, k, j, want_error);
     if (cu->size() != 1)
         throw "BUG: vector has " + S(cu->size()) + " elements while exactly one was expected";
     return cu->to_obj(0);
 }
 
 
+//! Vectorial evaluation of a curve's prange given by an expression tree.
+
+RObjVecInt COlc::eval_prange(const vector<double>& vt, int k, int j) const
+{
+    if (!prange_T)
+        return RObjVecInt( new CObjVecInt(vt.size(), 1) );
+    CContext ctx(k, j);
+    ctx.want_error = false;
+    RObj ret;
+    if (!prange_T->has_dummy()) { // curve does not depend on t
+        ret = prange_T->tree_val(ctx);
+    } else {
+        ctx.request_VT(&vt);
+        ret = prange_T->tree_val(ctx);
+    }
+    RObjVecInt ret2 = PCAST<const CObjVecInt>(ret);
+    if (!ret2)
+        throw S("BUG in eval_prange_expr: result is not an integer vector");
+    return ret2;
+}
+
+
 //**************************************************************************************************
 //*  class COlc - inspect
 //**************************************************************************************************
@@ -667,6 +690,10 @@ vector<string> COlc::info_settings() const
         txt += "conv file: " + S(kconv) + ", ";
     txt += "weighing: " + weight_str();
     ret.push_back(txt);
+    if (frange_expr!="")
+        ret.push_back("fit range: " + frange_expr);
+    if (prange_expr!="")
+        ret.push_back("plot range: " + prange_expr);
     if (chi2 > 0)
         ret.push_back("global fit result: " + S(chi2));
     return ret;

pub/lib/olf.hpp

@@ -118,8 +118,10 @@ public:
     static string wgtNames[];
     static TWgt name2wgt(string name);
     bool plot_to_grid; ///< Use grid from file 'kd' ?
-    string range_expr; ///< Restricts points to be fitted.
-    RNode range_T; ///< Parsed range_expr.
+    string frange_expr; ///< Restricts points to be fitted.
+    string prange_expr; ///< Restricts points to be plotted.
+    RNode frange_T; ///< Parsed range_expr.
+    RNode prange_T; ///< Parsed range_expr.
     double chi2; ///< Global chi^2
 
     COlc()
@@ -130,7 +132,8 @@ public:
         , kd(-1)
         , weighing(_VAR)
         , plot_to_grid(false)
-        , range_expr("")
+        , frange_expr("")
+        , prange_expr("")
         , chi2(0)
     {
     }
@@ -161,6 +164,7 @@ public:
     void curve_query(const string& quest);
     RObjVec eval_curve(const vector<double>& vt, int k, int j, bool want_error) const;
     RObj eval_curve_scalar(double arg, int k, int j, bool want_error) const;
+    RObjVecInt eval_prange(const vector<double>& vt, int k, int j) const;
     virtual string type() const { return "curve"; }
 };
 

pub/lib/plot.cpp

@@ -264,149 +264,56 @@ int plot_curve_equidist(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
     return xp.size();
 }
 
+class CPt {
+public:
+    double xd; // data point
+    double yd;
+    bool ok;
+    double xp; // plot point
+    double yp;
+};
+
 //! Plots slice j of non-convolved curve file fc, refining the grid where necessary;
 //! return number of points plotted.
 
 int plot_curve_refine(CPlot* plot, const COlc* fc, int k, int j, int cstyle)
 {
-    // refinement cares to:
-    // - stop/start line when leaving/reentering/crossing(TODO) the frame
-    // - TODO interpolate when crossing the horizontal border
-    // - interpolate at discontinuities
+    auto x2p = [&](const vector<double>& x) -> vector<CPt> {
+        vector<double> y = PCAST<const CObjVecDbl>(fc->eval_curve(x, k, j, false))->v;
+        vector<int> inRange = fc->eval_prange(x, k, j)->v;
+        vector<CPt> ret(x.size());
+        for (size_t i=0; i<x.size(); ++i) {
+            if (!plot->Y.contains(y[i]) || !inRange[i])
+                ret[i] = {x[i], 0., false, 0., 0.};
+            else
+                ret[i] = {x[i], y[i], true, plot->X.pc(x[i]), plot->Y.pc(y[i])};
+        }
+        return ret;
+    };
 
-    vector<double> xc, yc, xn, yn;
-    vector<double> novec;
+    vector<double> xc, yc;
+    vector<CPt> pc;
 
     // start with equidistant grid:
     int npts = plot->equipoints;
     plot->X.set_xgrid(xc, npts);
+    pc = x2p(xc);
 
-    // refinement loop:
-    for (int iref = 0;; ++iref) {
-        if (!xc.size()) {
-            if (!iref)
-                throw S("BUG: plot_curve_refine called with xc.size()=0");
-            break; // regular exit
-        }
-        // evaluate curve for grid xc:
-        RObjVecDbl cu = PCAST<const CObjVecDbl>(fc->eval_curve(xc, k, j, false));
-        for (int i = 0; i < xc.size(); ++i)
-            yc = cu->v;
-
-        // merge xc,yc and xn,yn:
-        if (iref == 0) {
-            xn = xc;
-            yn = yc;
-        } else {
-            vector<double> xa, ya;
-            int ic = 0, in = 0;
-            for (; ic < xc.size() && in < xn.size();) {
-                if (plot->X.close_enough(xc[ic], xn[ic], 1e-7)) {
-                    ++in; // skip new point
-                } else if (xc[ic] < xn[in]) {
-                    xa.push_back(xc[ic]);
-                    ya.push_back(yc[ic]);
-                    ++ic;
-                } else {
-                    xa.push_back(xn[in]);
-                    ya.push_back(yn[in]);
-                    ++in;
-                }
-            }
-            for (; ic < xc.size(); ++ic) {
-                xa.push_back(xc[ic]);
-                ya.push_back(yc[ic]);
-            }
-            for (; in < xn.size(); ++in) {
-                xa.push_back(xn[in]);
-                ya.push_back(yn[in]);
-            }
-            xn = xa;
-            yn = ya;
-        }
-        if (iref >= 10 || xn.size() >= plot->maxpoints)
-            break; // alternate exit: too many refinements or too many points
-
-        // eliminate points outside y range (except border points):
-        vector<double> xa, ya;
-        xa.push_back(xn[0]);
-        ya.push_back(yn[0]);
-        for (int i = 1; i < xn.size() - 1; ++i) {
-        new_segment:
-            xa.push_back(xn[i]);
-            ya.push_back(yn[i]);
-            if (!plot->Y.contains(yn[i])) {
-                for (int ii = i + 1; ii < xn.size() - 1; ++ii) {
-                    if (plot->Y.contains(yn[ii])) {
-                        if (ii - 1 > i) {
-                            xa.push_back(xn[ii - 1]);
-                            ya.push_back(yn[ii - 1]);
-                        }
-                        i = ii;
-                        goto new_segment;
-                    }
-                }
-            }
-        }
-        xa.push_back(xn.back());
-        ya.push_back(yn.back());
-        xn = xa;
-        yn = ya;
-
-        // set additional base points:
-        bool insert_next = false;
-        xc.clear();
-        for (int i = 1; i < xn.size() - 1; ++i) {
-            if (i < xn.size() - 2) {
-                double yi = ((xn[i + 1] - xn[i]) * yn[i - 1] + (xn[i] - xn[i - 1]) * yn[i + 1])
-                    / (xn[i + 1] - xn[i - 1]);
-                if (!plot->Y.close_enough(yn[i], yi, 0.005)) {
-                    xc.push_back((xn[i - 1] + xn[i]) / 2);
-                    insert_next = true;
-                    continue;
-                }
-            }
-            if (insert_next) {
-                xc.push_back((xn[i - 1] + xn[i]) / 2);
-                insert_next = false;
-                continue;
-            }
-            if ((plot->Y.contains(yn[i]) && !plot->Y.contains(yn[i - 1]))
-                || (!plot->Y.contains(yn[i]) && plot->Y.contains(yn[i - 1]))) {
-                xc.push_back((xn[i - 1] + xn[i]) / 2);
-                continue;
-            }
+    vector<double> xf, yf;
+    for (CPt p: pc) {
+        if (p.ok) {
+            xf.push_back(p.xd);
+            yf.push_back(p.yd);
         }
-        // temporary check
-        if (xc.size())
-            for (int i = 0; i < xc.size() - 1; ++i)
-                if (xc[i + 1] <= xc[i])
-                    throw "BUG: new base points not sorted at " + S(i) + ": " + S(xc[i]) + " "
-                        + S(xc[i + 1]);
     }
 
-    // divide into segments and plot:
-    vector<double> xa, ya;
     bool first_seg = true;
-    npts = 0;
-    for (int i = 0; i < xn.size(); ++i) {
-        if (plot->Y.contains(yn[i])) {
-            xa.push_back(xn[i]);
-            ya.push_back(yn[i]);
-        }
-        if (!plot->Y.contains(yn[i]) || i == xn.size() - 1) {
-            npts += xa.size();
-            if (xa.size()) {
-                plot->add_spec(
-                    true, first_seg, cstyle, xa, ya, novec, fc->V[j]->z_str(), fc->xco.str_std(),
-                    fc->yco.str_std(), "curve file " + S(k) + " spectrum " + S(j));
-                xa.clear();
-                ya.clear();
-                first_seg = false;
-            }
-        }
-    }
-    return npts;
+    if (xf.size())
+        plot->add_spec(
+            true, first_seg, cstyle, xf, yf, vector<double>(), fc->V[j]->z_str(), fc->xco.str_std(),
+            fc->yco.str_std(), "curve file " + S(k) + " spectrum " + S(j));
+
+    return xf.size();
 }