Files
GoldenCheetah/test/charts/eCP (LM).gchart
Mark Liversedge 972919fb87 Some .gcharts to test
.. added to test/charts

[skip ci]
2016-05-23 15:34:31 +01:00

19 lines
5.3 KiB
Plaintext

{
"CHART":{
"VERSION":"1",
"VIEW":"home",
"TYPE":"39",
"PROPERTIES":{
"title":"eCP (LM) ",
"subtitle":" ",
"widthFactor":"2",
"heightFactor":"2",
"style":"0",
"resizable":"0",
"script":"##\n## LM fit of eCP model\n##\n## And visualise all 3 energy systems\n##\nlibrary(minpack.lm)\n##library(corrplot)\n\nGC.page(height=600, width=900)\npar(bty=\"n\")\n## two plots are shown\n## 1) PD model fit and 3 energy systems\n## 2) Residuals from the fit\n##par(mfrow=c(1,2))\n\n##-------------------------------------------------\n## THE eCP MODEL\n##-------------------------------------------------\n\n# eCP function used to compute the residuals during \n# the fit -- note inf and nan are checked for as they\n# can occur with zero values for parms and will break\n# the regression fit, so we return zero instead\necp <- function(x,paa,paadec,cp,tau,taudel,cpdel,cpdec,cpdecdel) {\n\n x = x\/60.00\n\n rt = (\n\n paa*(1.20-0.20*exp(-1*x))*exp(paadec*(x))\n + ( cp * (1-exp(taudel*x)) * \n (1-exp(cpdel*x)) * \n (1+cpdec*exp(cpdecdel\/x)) * \n (tau\/(x)))\n + ( cp * (1-exp(taudel*x)) * \n (1-exp(cpdel*x)) * \n (1+cpdec*exp(cpdecdel\/x)) * \n (1))\n )\n\n if (is.infinite(rt) || is.nan(rt)) rt <- 0;\n return(rt)\n\n}\n\n## set an empty plot\nplot(x=c(1,1,1), y=c(1,1,1), pch=20,\n type=\"n\",\n xlim=c(1,5400),\n ylim=c(1,900),\n log=\"x\",\n xlab=\"seconds\", \n ylab=\"watts\")\n\n##-------------------------------------------------\n## PERFORM THE FIT\n##-------------------------------------------------\n\ncompares <- GC.season.meanmax(compare=TRUE)\n\nfor (compare in compares) {\n\n# get meanmax and make a dataframe\nmm <- head(compare$meanmax$power, n=5400)\ny <- tail(mm, length(mm)-1)\nx <- seq(1,length(y),1)\ndata <- data.frame(x,y)\n\n## fit using nls\nfit <- nlsLM(y~ecp(x,paa,paadec,cp,tau,taudel,\n cpdel,cpdec,cpdecdel), \n data = data,\n start = list(paa = 811.241,\n paadec = -2.08695, \n cp = 280.131,\n tau = 1.20846,\n taudel = -4.8,\n cpdel = -0.9,\n cpdec = -0.583937,\n cpdecdel = -180),\n trace = FALSE)\n\n# capture parameter estimates\nPAA <- coef(fit)[[1]]\nPAADEC <- coef(fit)[[2]]\nCP <- coef(fit)[[3]]\nTAU <- coef(fit)[[4]]\nTAUDEL <- coef(fit)[[5]]\nCPDEL <- coef(fit)[[6]]\nCPDEC <- coef(fit)[[7]]\nCPDECDEL <- coef(fit)[[8]]\n\n\n##-------------------------------------------------\n## FORMULAS TO PLOT CURVES\n##-------------------------------------------------\n\n## function to plot full PD model\ncp <- function(x) {\n return (ecp(x,PAA,PAADEC,CP,TAU,\n TAUDEL,CPDEL,CPDEC,CPDECDEL))\n}\n\n## aerobic\nacp <- function(x,paa,paadec,cp,tau,taudel,cpdel,cpdec,cpdecdel) {\n\n x = x\/60.00\n rt = (CP * (1-exp(TAUDEL*x)) * \n (1-exp(CPDEL*x)) * \n (1+CPDEC*exp(CPDECDEL\/x)) * \n (1))\n\n if (is.infinite(rt) || is.nan(rt)) rt <- 0;\n return(rt)\n\n}\n\n## immediate\nicp <- function(x) {\n x = x\/60.00\n rt = PAA *(1.20-0.20*exp(-1*x))*exp(PAADEC*(x))\n if (is.infinite(rt) || is.nan(rt)) rt <- 0;\n return(rt)\n}\n\n## intermediate\nncp <- function(x) {\n x = x\/60.00\n rt = (CP * (1-exp(TAUDEL*x)) * \n (1-exp(CPDEL*x)) * \n (1+CPDEC*exp(CPDECDEL\/x)) * \n (TAU\/(x)))\n if (is.infinite(rt) || is.nan(rt)) rt <- 0;\n return(rt)\n}\n\n\n\n##-------------------------------------------------\n## PLOT THE INPUTS AND RESULTS\n##-------------------------------------------------\n\n## plot the mean max first, so we can overlay\n## the model curves\nlines(x=seq(1,length(y),1), y=y, pch=20, \n add=TRUE, log=\"x\", type=\"p\",\n col=compare$color)\n\n# grid lines\ngrid(col=\"#404040\", lty=\"solid\")\n\n## plot the eCP curve using the fit object\nlines(x,fitted(fit), col=compare$color, lty=\"dotted\")\n\n## 3 energy systems\nt <- length(mm)\n#curve(acp, from=1, to=t, add=TRUE, \n# lty=\"dashed\", col=\"gold\")\n\n##curve(icp, from=1, to=t, add=TRUE, \n## lty=\"dashed\", col=\"cyan\")\n\n#curve(ncp, from=1, to=t, add=TRUE, \n# lty=\"dashed\", col=\"green\")\n\n# plot both glycolysis systems together\n##fn <- function(x) { return(ncp(x)+acp(x)) }\n##curve(fn,\n## from=1, to=t, add=TRUE,\n## lty=\"solid\", col=\"white\")\n\n\n## display results on the plot\nif (length(compares)==1) {\n cp <- coef(fit)[[3]]\n wprime <- cp * coef(fit)[[4]] * 60\n text(x=1000,y=700, paste(\"CP=\", round(cp,0)))\n text(x=1000,y=600, paste(\"W'=\", round(wprime,0)))\n}\n\n##-------------------------------------------------\n## PLOT THE RESIDUALS\n##-------------------------------------------------\n#res <- resid(fit)\n#plot(y,res,col=\"green\",pch=20, \n# ylab=\"Watts\", xlab=\"Residuals\")\n#grid(lty=\"solid\", col=\"#404040\")\n\n##-------------------------------------------------\n## CORRELATION MATRIX\n##-------------------------------------------------\n\n##cm <- summary(fit, correlation=TRUE)$correlation\n##corrplot(cm, bg=\"black\", order=\"hclust\")\n\n}\n\n ",
"state":" ",
"showConsole":"0",
"__LAST__":"1",
}
}
}