Dan Connelly's MEGA patch.

It includes both powerzones and weekly summary plots.

  Thanks Dan.
This commit is contained in:
Justin F. Knotzke
2009-06-22 02:21:12 +00:00
parent ea58d961a6
commit b2acd45c8d

View File

@@ -16,20 +16,29 @@
* Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include "Zones.h"
#include "CpintPlot.h"
#include <assert.h>
#include <QDebug>
#include <qwt_data.h>
#include <qwt_legend.h>
#include <qwt_plot_curve.h>
#include <qwt_plot_grid.h>
#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 <QwtPlotCurve *>::QList();
allZoneLabels= QList <QwtPlotMarker *>::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<cpi_file_info> &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<QString> 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<RideFilePoint*> 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 <double> 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<double> &bests, QVector<QDate> &bestDates)
read_one(const char *inname, QVector<double> &bests, QVector<QDate> &bestDates, QHash <QString, bool> *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<double> &bests, QVector<QDate> &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<double> &bests, QVector<QDate> &bestDates)
read_cpi_file(const QDir &dir, const QFileInfo &raw, QVector<double> &bests, QVector<QDate> &bestDates, QHash <QString, bool> *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<QwtPlotCurve *> 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<QwtPlotMarker *> 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 <int> power_zone;
int n_zones = (*zones)->lowsFromCP(&power_zone, (int) int(cp));
QList <int> 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<cpi_file_info> 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<double> 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<double> bests;
QVector<QDate> 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();
}