From fec4b2983ca3cc178532302b708f9b1c2f557a1d Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Wed, 8 Mar 2017 12:38:15 +0100
Subject: [PATCH] new function plotter respects prange ... but still does not
bisect at limits.
---
pub/itest/curveplot.f2t | 15 ++++
pub/lib/commands.cpp | 8 +-
pub/lib/curve.cpp | 29 +++++--
pub/lib/curve.hpp | 2 +-
pub/lib/fit.cpp | 4 +-
pub/lib/obj.hpp | 2 +-
pub/lib/olf.cpp | 35 ++++++++-
pub/lib/olf.hpp | 10 ++-
pub/lib/plot.cpp | 163 +++++++++-------------------------------
9 files changed, 118 insertions(+), 150 deletions(-)
diff --git a/pub/itest/curveplot.f2t b/pub/itest/curveplot.f2t
index 59b3acd9..0b723b0a 100755
--- a/pub/itest/curveplot.f2t
+++ b/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
diff --git a/pub/lib/commands.cpp b/pub/lib/commands.cpp
index 22bb2dc5..0e0132a6 100644
--- a/pub/lib/commands.cpp
+++ b/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();
diff --git a/pub/lib/curve.cpp b/pub/lib/curve.cpp
index ece8756d..b3f2265e 100644
--- a/pub/lib/curve.cpp
+++ b/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());
}
}
diff --git a/pub/lib/curve.hpp b/pub/lib/curve.hpp
index 9b0c781b..433de2e8 100644
--- a/pub/lib/curve.hpp
+++ b/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);
diff --git a/pub/lib/fit.cpp b/pub/lib/fit.cpp
index 49247480..5daae25d 100644
--- a/pub/lib/fit.cpp
+++ b/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] ) {
diff --git a/pub/lib/obj.hpp b/pub/lib/obj.hpp
index 7a049a56..9a8e6f1a 100644
--- a/pub/lib/obj.hpp
+++ b/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 ) ); }
diff --git a/pub/lib/olf.cpp b/pub/lib/olf.cpp
index 23008a1d..e403d608 100644
--- a/pub/lib/olf.cpp
+++ b/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;
diff --git a/pub/lib/olf.hpp b/pub/lib/olf.hpp
index e343cb38..327b1f6d 100644
--- a/pub/lib/olf.hpp
+++ b/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"; }
};
diff --git a/pub/lib/plot.cpp b/pub/lib/plot.cpp
index 6db000e2..68478e16 100644
--- a/pub/lib/plot.cpp
+++ b/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();
}
--
GitLab