diff --git a/src/CpintPlot.cpp b/src/CpintPlot.cpp index 56796aa05..a1f3bc793 100644 --- a/src/CpintPlot.cpp +++ b/src/CpintPlot.cpp @@ -16,20 +16,29 @@ * Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ +#include "Zones.h" #include "CpintPlot.h" - #include +#include #include #include #include #include +#include "RideItem.h" #include "LogTimeScaleDraw.h" #include "LogTimeScaleEngine.h" #include "RideFile.h" -CpintPlot::CpintPlot(QString p) : - progress(NULL), needToScanRides(true), path(p), allCurve(NULL), - thisCurve(NULL), grid(NULL) +CpintPlot::CpintPlot( + QString p + ) : + progress(NULL), + needToScanRides(true), + path(p), + thisCurve(NULL), + CPCurve(NULL), + grid(NULL), + zones(NULL) { insertLegend(new QwtLegend(), QwtPlot::BottomLegend); setCanvasBackground(Qt::white); @@ -39,45 +48,55 @@ CpintPlot::CpintPlot(QString p) : setAxisScaleEngine(xBottom, new LogTimeScaleEngine); setAxisScale(xBottom, 1.0 / 60.0, 60); - allCurve = new QwtPlotCurve("All Rides"); - allCurve->setRenderHint(QwtPlotItem::RenderAntialiased); - allCurve->setPen(QPen(Qt::red)); - allCurve->attach(this); - grid = new QwtPlotGrid(); grid->enableX(false); QPen gridPen; gridPen.setStyle(Qt::DotLine); grid->setPen(gridPen); grid->attach(this); -} + + allCurves= QList ::QList(); + allZoneLabels= QList ::QList(); + } struct cpi_file_info { QString file, inname, outname; }; +bool +is_ride_filename(const QString filename) +{ + QRegExp re("^([0-9][0-9][0-9][0-9])_([0-9][0-9])_([0-9][0-9])" + "_([0-9][0-9])_([0-9][0-9])_([0-9][0-9])\\.(raw|srm|csv|tcx)$"); + return (re.exactMatch(filename)); +} + +QString +ride_filename_to_cpi_filename(const QString filename) +{ + return (QFileInfo(filename).completeBaseName() + ".cpi"); +} + static void cpi_files_to_update(const QDir &dir, QList &result) { - QRegExp re("^([0-9][0-9][0-9][0-9])_([0-9][0-9])_([0-9][0-9])" - "_([0-9][0-9])_([0-9][0-9])_([0-9][0-9])\\.(raw|srm|csv|tcx)$"); QStringList filenames = RideFileFactory::instance().listRideFiles(dir); QListIterator i(filenames); while (i.hasNext()) { const QString &filename = i.next(); - if (re.exactMatch(filename)) { - QString inname = dir.absoluteFilePath(filename); - QString outname = dir.absoluteFilePath( - QFileInfo(filename).completeBaseName() + ".cpi"); - QFileInfo ifi(inname), ofi(outname); - if (!ofi.exists() || (ofi.lastModified() < ifi.lastModified())) { - cpi_file_info info; - info.file = filename; - info.inname = inname; - info.outname = outname; - result.append(info); - } - } + if (is_ride_filename(filename)) { + QString inname = dir.absoluteFilePath(filename); + QString outname = + dir.absoluteFilePath(ride_filename_to_cpi_filename(filename)); + QFileInfo ifi(inname), ofi(outname); + if (!ofi.exists() || (ofi.lastModified() < ifi.lastModified())) { + cpi_file_info info; + info.file = filename; + info.inname = inname; + info.outname = outname; + result.append(info); + } + } } } @@ -103,14 +122,15 @@ update_cpi_file(const cpi_file_info *info, QProgressDialog *progress, QStringList errors; RideFile *rideFile = RideFileFactory::instance().openRideFile(file, errors); - assert(rideFile); + if (! rideFile) + return; cpint_data data; data.rec_int_ms = (int) round(rideFile->recIntSecs() * 1000.0); QListIterator i(rideFile->dataPoints()); while (i.hasNext()) { - const RideFilePoint *p = i.next(); - double secs = round(p->secs * 1000.0) / 1000; - data.points.append(cpint_point(secs, (int) round(p->watts))); + const RideFilePoint *p = i.next(); + double secs = round(p->secs * 1000.0) / 1000; + data.points.append(cpint_point(secs, (int) round(p->watts))); } delete rideFile; @@ -118,11 +138,20 @@ update_cpi_file(const cpi_file_info *info, QProgressDialog *progress, assert(out); int total_secs = (int) ceil(data.points.back().secs); - if (total_secs > 7*24*60*60) { + + // don't allow data more than one week + #define SECONDS_PER_WEEK 7 * 24 * 60 * 60 + if (total_secs > SECONDS_PER_WEEK) { fclose(out); return; } - double *bests = (double*) calloc(total_secs + 1, sizeof(double)); + + + QVector ride_bests(total_secs + 1); // was calloc'ed array (unfreed?), changed djconnel + + // initialize ride_bests + for (int i = 0; i < ride_bests.size(); i ++) + ride_bests[i] = 0.0; bool canceled = false; int progress_count = 0; @@ -151,16 +180,27 @@ update_cpi_file(const cpi_file_info *info, QProgressDialog *progress, int dur_secs_bot = qMax((int) floor(dur_secs - data.rec_int_ms / 1000.0), 0); for (int k = dur_secs_top; k > dur_secs_bot; --k) { - if (bests[k] < avg) - bests[k] = avg; + if (ride_bests[k] < avg) + ride_bests[k] = avg; } prev_secs = q->secs; } } - for (int i = 1; i <= total_secs; ++i) { - if (bests[i] != 0) - fprintf(out, "%6.3f %3.0f\n", i / 60.0, round(bests[i])); + // avoid decreasing work with increasing duration + { + double maxwork = 0.0; + for (int i = 1; i <= total_secs; ++i) { + // note index is being used here in lieu of time, as the index + // is assumed to be proportional to time + double work = ride_bests[i] * i; + if (maxwork > work) + ride_bests[i] = round(maxwork / i); + else + maxwork = work; + if (ride_bests[i] != 0) + fprintf(out, "%6.3f %3.0f\n", i / 60.0, round(ride_bests[i])); + } } done: @@ -170,25 +210,37 @@ done: progress_sum += progress_count; } +QDate +cpi_filename_to_date(const QString filename) { + QRegExp rx(".*([0-9][0-9][0-9][0-9])_([0-9][0-9])_([0-9][0-9])" + "_([0-9][0-9])_([0-9][0-9])_([0-9][0-9])\\.cpi$"); + if (rx.exactMatch(filename)) { + assert(rx.numCaptures() == 6); + QDate date = QDate( + rx.cap(1).toInt(), + rx.cap(2).toInt(), + rx.cap(3).toInt() + ); + + if (! date.isValid()) { + return QDate(); + } + else + return date; + } + else + return QDate(); // return value was 1 Jan: changed to null +} + static int -read_one(const char *inname, QVector &bests, QVector &bestDates) +read_one(const char *inname, QVector &bests, QVector &bestDates, QHash *cpiDataInBests) { FILE *in = fopen(inname, "r"); if (!in) return -1; int lineno = 1; char line[40]; - - QRegExp rx(".*([0-9][0-9][0-9][0-9])_([0-9][0-9])_([0-9][0-9])" - "_([0-9][0-9])_([0-9][0-9])_([0-9][0-9])\\.cpi$"); - QDate date; - if (rx.exactMatch(inname)) { - assert(rx.numCaptures() == 6); - date = QDate(rx.cap(1).toInt(), rx.cap(2).toInt(),rx.cap(3).toInt()); - } - else - date = QDate(2009,1,1); - + while (fgets(line, sizeof(line), in) != NULL) { double mins; int watts; @@ -203,28 +255,393 @@ read_one(const char *inname, QVector &bests, QVector &bestDates) } if (bests[secs] < watts){ bests[secs] = watts; - bestDates[secs] = date; + bestDates[secs] = cpi_filename_to_date(QString(inname)); + + // mark the filename as having contributed to the bests + // Note this contribution may subsequently be over-written, so + // for example the first file scanned will always be tagged. + if (cpiDataInBests) + (*cpiDataInBests)[inname] = true; } ++lineno; } fclose(in); + return 0; } static int -read_cpi_file(const QDir &dir, const QFileInfo &raw, QVector &bests, QVector &bestDates) +read_cpi_file(const QDir &dir, const QFileInfo &raw, QVector &bests, QVector &bestDates, QHash *cpiDataInBests) { QString inname = dir.absoluteFilePath(raw.completeBaseName() + ".cpi"); - return read_one(inname.toAscii().constData(), bests, bestDates); + return read_one(inname.toAscii().constData(), bests, bestDates, cpiDataInBests); +} + +// extract critical power parameters which match the given curve +// model: maximal power = cp (1 + tau / [t + t0]), where t is the +// duration of the effort, and t, cp and tau are model parameters +// the basic critical power model is t0 = 0, but non-zero has +// been discussed in the literature +// it is assumed duration = index * seconds +void +CpintPlot::deriveCPParameters() +{ + // bounds on anaerobic interval in minutes + const double t1 = USE_T0_IN_CP_MODEL ? 0.25 : 1; + const double t2 = 6; + + // bounds on aerobic interval in minutes + const double t3 = 10; + const double t4 = 60; + + // bounds of these time valus in the data + int i1, i2, i3, i4; + + // find the indexes associated with the bounds + // the first point must be at least the minimum for the anaerobic interval, or quit + for (i1 = 0; i1 < 60 * t1; i1++) + if (i1 + 1 >= bests.size()) + return; + // the second point is the maximum point suitable for anaerobicly dominated efforts. + for (i2 = i1; i2 + 1 <= 60 * t2; i2++) + if (i2 + 1 >= bests.size()) + return; + // the third point is the beginning of the minimum duration for aerobic efforts + for (i3 = i2; i3 < 60 * t3; i3++) + if (i3 + 1 >= bests.size()) + return; + for (i4 = i3; i4 + 1 <= 60 * t4; i4++) + if (i4 + 1 >= bests.size()) + break; + + // initial estimate of tau + if (tau == 0) + tau = 1; + + // initial estimate of cp (if not already available) + if (cp == 0) + cp = 300; + + // initial estimate of t0: start small to maximize sensitivity to data + t0 = 0; + + // lower bound on tau + const double tau_min = 0.5; + + // convergence delta for tau + const double tau_delta_max = 1e-4; + const double t0_delta_max = 1e-4; + + // previous loop value of tau and t0 + double tau_prev; + double t0_prev; + + // maximum number of loops + const int max_loops = 100; + + // loop to convergence + int iteration = 0; + do { + if (iteration ++ > max_loops) { + fprintf(stderr, "maximum number of loops %d exceeded in cp model extraction\n", max_loops); + break; + } + + // record the previous version of tau, for convergence + tau_prev = tau; + t0_prev = t0; + + // estimate cp, given tau + int i; + cp = 0; + for (i = i3; i <= i4; i++) { + double cpn = bests[i] / (1 + tau / (t0 + i / 60.0)); + if (cp < cpn) + cp = cpn; + } + + // if cp = 0; no valid data; give up + if (cp == 0.0) + return; + + // estimate tau, given cp + tau = tau_min; + for (i = i1; i <= i2; i++) { + double taun = (bests[i] / cp - 1) * (i / 60.0 + t0) - t0; + if (tau < taun) + tau = taun; + } + + // update t0 if we're using that model + #if USE_T0_IN_CP_MODEL + t0 = tau / (bests[1] / cp - 1) - 1 / 60.0; + #endif + + // the following line is debugging code and can be removed + fprintf(stderr, "%d: tau = %.2f; cp = %.2f; t0 = %.2f\n", iteration, tau, cp, t0); + + } while ((fabs(tau - tau_prev) > tau_delta_max) || + (fabs(t0 - t0_prev) > t0_delta_max) + ); + +} + + + +void +CpintPlot::plot_CP_curve( + CpintPlot *thisPlot, // the plot we're currently displaying + double cp, + double tau, + double t0 + ) { + + // detach the CP curve if it exists + if (CPCurve) + CPCurve->detach(); + + // if there's no cp, then there's nothing to do + if (cp <= 0) + return; + + // populate curve data with a CP curve + const int curve_points = 100; + double tmin = USE_T0_IN_CP_MODEL ? 1.0/60 : tau; + double tmax = 180.0; + double cp_curve_power[curve_points]; + double cp_curve_time[curve_points]; + int i; + + for (i = 0; i < curve_points; i ++) { + double x = (double) i / (curve_points - 1); + double t = pow(tmax, x) * pow(tmin, 1-x); + cp_curve_time[i] = t; + cp_curve_power[i] = cp * (1 + tau / (t + t0)); + } + + // generate a plot + QString curve_title; +#if USE_T0_IN_CP_MODEL + curve_title.sprintf("CP=%.1f W; AWC/CP=%.2f m; t0=%.1f s", cp, tau, 60 * t0); +#else + curve_title.sprintf("CP=%.1f W; AWC/CP=%.2f m", cp, tau); +#endif + + CPCurve = new QwtPlotCurve(curve_title); + CPCurve->setRenderHint(QwtPlotItem::RenderAntialiased); + QPen *pen = new QPen(Qt::red); + pen->setWidth(2.0); + pen->setStyle(Qt::DashLine); + CPCurve->setPen(*pen); + CPCurve->setData( + cp_curve_time, + cp_curve_power, + curve_points + ); + CPCurve->attach(thisPlot); + delete pen; + + return; +} + +void CpintPlot::clear_CP_Curves() { + // unattach any existing shading curves and reset the list + if (allCurves.size()) { + QListIterator i(allCurves); + while (i.hasNext()) { + QwtPlotCurve *curve = i.next(); + if (curve) { + curve->detach(); + delete curve; + } + } + allCurves.clear(); + } + + // now delete any labels + if (allZoneLabels.size()) { + QListIterator i(allZoneLabels); + while (i.hasNext()) { + QwtPlotMarker *label = i.next(); + if (label) { + label->detach(); + delete label; + } + } + allZoneLabels.clear(); + } +} + + +void +CpintPlot::plot_allCurve ( + CpintPlot *thisPlot, + int n_values, + const double *power_values + ) { + + clear_CP_Curves(); + + // generate an array of time values + double time_values[n_values]; + for (int t = 1; t <= n_values; t++) + time_values[t - 1] = t / 60.0; + + // generate zones from derived CP value + if (cp > 0) { + QList power_zone; + int n_zones = (*zones)->lowsFromCP(&power_zone, (int) int(cp)); + + QList n_zone; + + // the lowest zone goes to zero power, so mark its start at the last data point + n_zone.append(n_values - 1); + + // start the search at the next-to-lowest zone + int z = 1; + + // search the maximal power curve to extract the zone times + for (int i = n_values; i-- > 0;) { + // if we reach the beginning of the curve OR if we hit a zone boundary, we're done with the present zone + if ((i == 0) || (power_values[i] > power_zone[z])) { + n_zone.append( + (z == n_zones) ? + 0 : + ( + ( + (i == n_values - 1) || + (abs(power_values[i] - power_zone[z]) < abs(power_zone[z] - power_values[i + 1])) + ) ? + i : + i + 1 + ) + ); + + // draw curves for the zone we're leaving, if it spans any segments + if (n_zone[z - 1] > n_zone[z]) { + // define the individual code segments. Note in the old code with a single segment, it was + // part of the class. This curve is not a protected member of the class. djconnel Apr2009 + QwtPlotCurve *curve; + curve = + new QwtPlotCurve((*zones)->getDefaultZoneName(z - 1)); + curve->setRenderHint(QwtPlotItem::RenderAntialiased); + QPen *pen = new QPen(zoneColor(z - 1, n_zones)); + pen->setWidth(2.0); + curve->setPen(*pen); + curve->attach(thisPlot); + QColor brush_color = zoneColor(z - 1, n_zones); + brush_color.setAlpha(64); + curve->setBrush(brush_color); // brush fills below the line + curve->setData( + time_values + n_zone[z], + power_values + n_zone[z], + n_zone[z - 1] - n_zone[z] + 1 + ); + delete pen; + + // add the curve to the list + allCurves.append(curve); + + // render a colored label on the zone + QwtText text((*zones)->getDefaultZoneName(z - 1)); + text.setFont(QFont("Helvetica",24, QFont::Bold)); + QColor text_color = zoneColor(z - 1, n_zones); + text_color.setAlpha(128); + text.setColor(text_color); + QwtPlotMarker *label_mark; + label_mark = new QwtPlotMarker(); + + // place the text in the geometric mean in time, at a decent power + label_mark->setValue( + sqrt(time_values[n_zone[z-1]] * time_values[n_zone[z]]), + (power_values[n_zone[z-1]] + power_values[n_zone[z]]) / 5 + ); + label_mark->setLabel(text); + label_mark->attach(thisPlot); + allZoneLabels.append(label_mark); + } + + if (z < n_zones) + fprintf(stderr, "zone %s: %d watts, index = %d\n", + (*zones)->getDefaultZoneName(z).toAscii().constData(), + power_zone[z], + n_zone[z] + ); + + // if we're to the smallest time, we're done + if (i == 0) + break; + + // increment zone number + if (z < n_zones) + z ++; + + // if we're to the final zone, just jump to the beginning of the plot: we're done + if (z == n_zones) + i = 1; + + // else, we've got to recheck this point for the next zone + else + i ++; + } + } + } + // no zones available: just plot the curve without zones + else { + QwtPlotCurve *curve; + curve = new QwtPlotCurve("maximal power"); + curve->setRenderHint(QwtPlotItem::RenderAntialiased); + QPen *pen = new QPen(Qt::red); + pen->setWidth(2.0); + curve->setPen(*pen); + QColor brush_color = Qt::red; + brush_color.setAlpha(64); + curve->setBrush(brush_color); // brush fills below the line + curve->setData( + time_values, + power_values, + n_values + ); + curve->attach(thisPlot); + delete pen; + allCurves.append(curve); + } + + + // set the x-axis to span the time of the all-time curve, starting at 1 second + thisPlot->setAxisScale( + thisPlot->xBottom, + 1.0 / 60, + time_values[n_values - 1] + ); + + // set the y-axis to go from zero to the maximum power, rounded up to nearest 100 watts + thisPlot->setAxisScale( + thisPlot->yLeft, + 0, + 100 * ceil( power_values[0] / 100 ) + ); + } void -CpintPlot::calculate(QString fileName, QDateTime dateTime) +CpintPlot::calculate(RideItem *rideItem) { + QString fileName = rideItem->fileName; + QDateTime dateTime = rideItem->dateTime; QDir dir(path); QFileInfo file(fileName); + zones = rideItem->zones; if (needToScanRides) { + bests.clear(); + bestDates.clear(); + cpiDataInBests.clear(); + if (CPCurve) { + CPCurve->detach(); + CPCurve = NULL; + } + fflush(stderr); bool aborted = false; QList to_update; @@ -238,10 +655,11 @@ CpintPlot::calculate(QString fileName, QDateTime dateTime) QStringList errors; RideFile *rideFile = RideFileFactory::instance().openRideFile(file, errors); - assert(rideFile); - double x = rideFile->dataPoints().size(); - progress_max += x * (x + 1.0) / 2.0; - delete rideFile; + if (rideFile) { + double x = rideFile->dataPoints().size(); + progress_max += x * (x + 1.0) / 2.0; + delete rideFile; + } } } progress = new QProgressDialog( @@ -268,7 +686,6 @@ CpintPlot::calculate(QString fileName, QDateTime dateTime) } } } - QVector bests; if (!aborted) { QString existing = progress->labelText(); existing.chop(progress->labelText().size() - endingOffset); @@ -284,7 +701,7 @@ CpintPlot::calculate(QString fileName, QDateTime dateTime) while (i.hasNext()) { const QString &filename = i.next(); QString path = dir.absoluteFilePath(filename); - read_one(path.toAscii().constData(), bests, bestDates); + read_one(path.toAscii().constData(), bests, bestDates, &cpiDataInBests); progress->setValue(progress->value() + 1); QCoreApplication::processEvents(); if (progress->wasCanceled()) { @@ -292,32 +709,55 @@ CpintPlot::calculate(QString fileName, QDateTime dateTime) break; } } - } - if (!aborted) { - double *timeArray = new double[bests.size()]; - int maxNonZero = 0; - for (int i = 0; i < bests.size(); ++i) { - timeArray[i] = i / 60.0; - if (bests[i] > 0) maxNonZero = i; - } - if (bests.size() > 1) { - allCurve->setData(timeArray + 1, bests.constData() + 1, - maxNonZero - 1); - setAxisScale(xBottom, 1.0 / 60.0, maxNonZero / 60.0); - } - delete [] timeArray; - needToScanRides = false; - } - delete progress; - progress = NULL; + } + if (!aborted && bests.size()) { + int maxNonZero = 0; + + // check that total work doesn't decrease with time + double maxwork = 0.0; + + for (int i = 0; i < bests.size(); ++i) { + // record the date associated with each point's CPI file, + if (bests[i] > 0) + maxNonZero = i; + + // note index is being used here in lieu of time, as the index + // is assumed to be proportional to time + double work = bests[i] * i; + if ((i > 0) && (maxwork > work)) { + bests[i] = round(maxwork / i); + bestDates[i] = bestDates[i - 1]; + } + else + maxwork = work; + } + + // derive CP model + if (bests.size() > 1) { + // cp model parameters + cp = 0; + tau = 0; + t0 = 0; + + // calculate CP model from all-time best data + deriveCPParameters(); + plot_CP_curve(this, cp, tau, t0); + + plot_allCurve(this, maxNonZero - 1, bests.constData() + 1); + } + needToScanRides = false; + } + delete progress; + progress = NULL; } if (!needToScanRides) { - delete thisCurve; + if (thisCurve) + delete thisCurve; thisCurve = NULL; QVector bests; QVector bestDates; - if (read_cpi_file(dir, file, bests, bestDates) == 0) { + if ((read_cpi_file(dir, file, bests, bestDates, NULL) == 0) && bests.size()) { double *timeArray = new double[bests.size()]; int maxNonZero = 0; for (int i = 0; i < bests.size(); ++i) { @@ -328,18 +768,42 @@ CpintPlot::calculate(QString fileName, QDateTime dateTime) thisCurve = new QwtPlotCurve( dateTime.toString("ddd MMM d, yyyy h:mm AP")); thisCurve->setRenderHint(QwtPlotItem::RenderAntialiased); - thisCurve->setPen(QPen(Qt::green)); + QPen *pen = new QPen(Qt::black); + pen->setWidth(2.0); + thisCurve->setPen(QPen(Qt::black)); thisCurve->attach(this); thisCurve->setData(timeArray + 1, bests.constData() + 1, maxNonZero - 1); } - delete [] timeArray; + delete [] timeArray; } } replot(); } +// delete a CPI file +bool +CpintPlot::deleteCpiFile(QString filename) +{ + // first, get ride of the file + if (! QFile::remove(filename)) + return false; + + // now check to see if this file contributed to the bests + // in the current implementation a false means it does + // not contribute, but a true only means it at one time + // contributed (may not in the end). + if (cpiDataInBests.contains(filename)) { + if (cpiDataInBests[filename]) + needToScanRides = true; + cpiDataInBests.remove(filename); + } + + return true; +} + + void CpintPlot::showGrid(int state) { @@ -347,4 +811,3 @@ CpintPlot::showGrid(int state) grid->setVisible(state == Qt::Checked); replot(); } -