From 6fa809139478e74851f94a5ff251fc4acc505993 Mon Sep 17 00:00:00 2001 From: Jiachen Lu Date: Thu, 15 Sep 2022 12:29:17 -0500 Subject: [PATCH 1/2] create folders --- fpca_NHANES.Rproj => code/fpca_NHANES.Rproj | 0 fpca_analysis.R => code/fpca_analysis.R | 0 plot_fd_new1.R => code/plot_fd_new1.R | 410 +++++++++--------- plot_fpca_new1.R => code/plot_fpca_new1.R | 180 ++++---- .../NHANES_sample_data.rds | Bin .../NHANES_sample_data_day_average.rds | Bin 6 files changed, 295 insertions(+), 295 deletions(-) rename fpca_NHANES.Rproj => code/fpca_NHANES.Rproj (100%) rename fpca_analysis.R => code/fpca_analysis.R (100%) rename plot_fd_new1.R => code/plot_fd_new1.R (96%) mode change 100755 => 100644 rename plot_fpca_new1.R => code/plot_fpca_new1.R (96%) mode change 100755 => 100644 rename NHANES_sample_data.rds => data/NHANES_sample_data.rds (100%) rename NHANES_sample_data_day_average.rds => data/NHANES_sample_data_day_average.rds (100%) diff --git a/fpca_NHANES.Rproj b/code/fpca_NHANES.Rproj similarity index 100% rename from fpca_NHANES.Rproj rename to code/fpca_NHANES.Rproj diff --git a/fpca_analysis.R b/code/fpca_analysis.R similarity index 100% rename from fpca_analysis.R rename to code/fpca_analysis.R diff --git a/plot_fd_new1.R b/code/plot_fd_new1.R old mode 100755 new mode 100644 similarity index 96% rename from plot_fd_new1.R rename to code/plot_fd_new1.R index 6f2f6c8..039f56b --- a/plot_fd_new1.R +++ b/code/plot_fd_new1.R @@ -1,205 +1,205 @@ -# 7/11/21 change color -plot.fd.new <- function (x, y, Lfdobj = 0, href = TRUE, titles = NULL, xlim = NULL, - ylim = NULL, xlab = NULL, ylab = NULL, ask = FALSE, nx = NULL, - axes = NULL, family = NULL, ...) { - fdobj <- x - if (!(is.fd(fdobj) || is.fdPar(fdobj))) - stop(paste("First argument is neither a functional data or a ", - "functional parameter object.")) - if (is.fdPar(fdobj)) - fdobj <- fdobj$fd - if (is.null(axes)) { - if (is.null(fdobj$basis$axes)) { - Axes <- TRUE - axFun <- FALSE - } - else { - if (!inherits(fdobj$basis$axes, "list")) - stop("fdobj$basis$axes must be a list; ", - "class(fdobj$basis$axes) = ", class(fdobj$basis$axes)) - if (!(inherits(fdobj$basis$axes[[1]], "character") || - inherits(fdobj$basis$axes[[1]], "function"))) - stop("fdobj$basis$axes[[1]] must be either a function or the ", - "name of a function; class(fdobj$basis$axes[[1]]) = ", - class(fdobj$basis$axes[[1]])) - Axes <- FALSE - axFun <- TRUE - axList <- c(fdobj$basis$axes, ...) - } - } - else { - if (is.logical(axes)) { - Axes <- axes - axFun <- FALSE - } - else { - if (!inherits(axes, "list")) - stop("axes must be a logical or a list; class(axes) = ", - class(axes)) - if (!(inherits(axes[[1]], "character") || inherits(axes[[1]], - "function"))) - stop("axes[[1]] must be either a function or the ", - "name of a function; class(axes[[1]]) = ", - class(axes[[1]])) - Axes <- FALSE - axFun <- TRUE - axList <- c(axes, ...) - } - } - Lfdobj <- int2Lfd(Lfdobj) - if (!inherits(Lfdobj, "Lfd")) - stop("Second argument is not a linear differential operator.") - coef <- fdobj$coefs - coefd <- dim(coef) - ndim <- length(coefd) - nbasis <- coefd[1] - if (is.null(nx)) - nx <- max(c(501, 10 * nbasis + 1)) - nrep <- coefd[2] - if (ndim > 2) - nvar <- coefd[3] - else nvar <- 1 - basisobj <- fdobj$basis - rangex <- basisobj$rangeval - if (missing(y)) { - y <- nx - } - else { - if (is.numeric(y)) - y <- as.vector(y) - } - Y <- y - if (length(y) == 1) { - if (y >= 1) { - y <- seq(rangex[1], rangex[2], len = round(y)) - } - else { - stop("'y' a single number less than one.") - } - } - if (min(y) < rangex[1] || max(y) > rangex[2]) - stop("Values in Y are outside the basis range.") - if (is.null(xlim)) { - xlim <- rangex - } - else { - rangex[1] <- max(rangex[1], xlim[1]) - rangex[2] <- min(rangex[2], xlim[2]) - if (length(Y) == 1) - y <- seq(rangex[1], rangex[2], len = round(Y)) - } - fdmat <- eval.fd(y, fdobj, Lfdobj) - rangey <- range(fdmat) - if (is.null(ylim)) - ylim <- rangey - fdnames = fdobj$fdnames - fdlabelslist = fdlabels(fdnames, nrep, nvar) - xlabel = fdlabelslist$xlabel - ylabel = fdlabelslist$ylabel - casenames = fdlabelslist$casenames - varnames = fdlabelslist$varnames - if (is.null(xlab)) - xlab <- xlabel - if (is.null(ylab)) - ylab <- ylabel - if (ndim < 2) { - plot(y, fdmat, type = "l", xlim = xlim, ylim = ylim, - xlab = xlab, ylab = ylab, axes = Axes - , xaxt="n", ...) - # xaxis labels and locations - x.axis <- seq(from=0, to=24, by=1) - x.axis.lab <- paste0(x.axis,":00") - x.axis.loc <- x.axis/24*288 - axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) - - if (axFun) - do.call(axList[[1]], axList[-1]) - if (zerofind(fdmat) && href) - abline(h = 0, lty = 2) - } - if (ndim == 2) { - if (!ask) { - matplot(y, fdmat, type = "l", xlim = xlim, - ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes - , xaxt="n" - #, col="blue1" - , family=family, ps=11 - , ... - ) - # xaxis labels and locations - # x.axis <- seq(from=0, to=24, by=1) - x.axis <- seq(from=0, to=24, by=4) - x.axis.lab <- paste0(x.axis,":00") - x.axis.loc <- x.axis/24*288 - axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) - - if (axFun) - do.call(axList[[1]], axList[-1]) - if (zerofind(fdmat) && href) - abline(h = 0, lty = 2) - } - else { - op <- par(ask = FALSE) - on.exit(par(op)) - cat("Multiple plots: Click in the plot to advance to the next") - for (irep in 1:nrep) { - plot(y, fdmat[, irep], type = "l", xlim = xlim, - ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes, - ...) - if (axFun) - do.call(axList[[1]], axList[-1]) - if (irep < 2) - par(ask = ask) - if (!is.null(casenames)) - title(casenames[irep]) - else title(paste("Case", irep)) - if (zerofind(ylim) && href) - abline(h = 0, lty = 2) - } - } - } - if (ndim == 3) { - if (!ask) { - for (ivar in 1:nvar) { - matplot(y, fdmat[, , ivar], type = "l", - xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, - ask = FALSE, axes = Axes, ...) - if (axFun) - do.call(axList[[1]], axList[-1]) - if (!is.null(varnames)) - title(varnames[ivar]) - else title(paste("Variable", ivar)) - if (zerofind(ylim) && href) - abline(h = 0, lty = 2) - } - } - else { - op <- par(ask = FALSE) - on.exit(par(op)) - cat("Multiple plots: Click in the plot to advance to the next") - for (irep in 1:nrep) { - for (ivar in 1:nvar) { - plot(y, fdmat[, irep, ivar], type = "l", - xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, - axes = Axes, ...) - if (axFun) - do.call(axList[[1]], axList[-1]) - if (!is.null(casenames)) - titlestr = casenames[irep] - else titlestr = paste("Case", irep) - if (!is.null(varnames)) { - titlestr = paste(titlestr, " ", varnames[ivar]) - } - else { - titlestr = paste(titlestr, " ", "Variable", - ivar) - } - title(titlestr) - if (zerofind(ylim) && href) - abline(h = 0, lty = 2) - } - } - } - } - "done" -} +# 7/11/21 change color +plot.fd.new <- function (x, y, Lfdobj = 0, href = TRUE, titles = NULL, xlim = NULL, + ylim = NULL, xlab = NULL, ylab = NULL, ask = FALSE, nx = NULL, + axes = NULL, family = NULL, ...) { + fdobj <- x + if (!(is.fd(fdobj) || is.fdPar(fdobj))) + stop(paste("First argument is neither a functional data or a ", + "functional parameter object.")) + if (is.fdPar(fdobj)) + fdobj <- fdobj$fd + if (is.null(axes)) { + if (is.null(fdobj$basis$axes)) { + Axes <- TRUE + axFun <- FALSE + } + else { + if (!inherits(fdobj$basis$axes, "list")) + stop("fdobj$basis$axes must be a list; ", + "class(fdobj$basis$axes) = ", class(fdobj$basis$axes)) + if (!(inherits(fdobj$basis$axes[[1]], "character") || + inherits(fdobj$basis$axes[[1]], "function"))) + stop("fdobj$basis$axes[[1]] must be either a function or the ", + "name of a function; class(fdobj$basis$axes[[1]]) = ", + class(fdobj$basis$axes[[1]])) + Axes <- FALSE + axFun <- TRUE + axList <- c(fdobj$basis$axes, ...) + } + } + else { + if (is.logical(axes)) { + Axes <- axes + axFun <- FALSE + } + else { + if (!inherits(axes, "list")) + stop("axes must be a logical or a list; class(axes) = ", + class(axes)) + if (!(inherits(axes[[1]], "character") || inherits(axes[[1]], + "function"))) + stop("axes[[1]] must be either a function or the ", + "name of a function; class(axes[[1]]) = ", + class(axes[[1]])) + Axes <- FALSE + axFun <- TRUE + axList <- c(axes, ...) + } + } + Lfdobj <- int2Lfd(Lfdobj) + if (!inherits(Lfdobj, "Lfd")) + stop("Second argument is not a linear differential operator.") + coef <- fdobj$coefs + coefd <- dim(coef) + ndim <- length(coefd) + nbasis <- coefd[1] + if (is.null(nx)) + nx <- max(c(501, 10 * nbasis + 1)) + nrep <- coefd[2] + if (ndim > 2) + nvar <- coefd[3] + else nvar <- 1 + basisobj <- fdobj$basis + rangex <- basisobj$rangeval + if (missing(y)) { + y <- nx + } + else { + if (is.numeric(y)) + y <- as.vector(y) + } + Y <- y + if (length(y) == 1) { + if (y >= 1) { + y <- seq(rangex[1], rangex[2], len = round(y)) + } + else { + stop("'y' a single number less than one.") + } + } + if (min(y) < rangex[1] || max(y) > rangex[2]) + stop("Values in Y are outside the basis range.") + if (is.null(xlim)) { + xlim <- rangex + } + else { + rangex[1] <- max(rangex[1], xlim[1]) + rangex[2] <- min(rangex[2], xlim[2]) + if (length(Y) == 1) + y <- seq(rangex[1], rangex[2], len = round(Y)) + } + fdmat <- eval.fd(y, fdobj, Lfdobj) + rangey <- range(fdmat) + if (is.null(ylim)) + ylim <- rangey + fdnames = fdobj$fdnames + fdlabelslist = fdlabels(fdnames, nrep, nvar) + xlabel = fdlabelslist$xlabel + ylabel = fdlabelslist$ylabel + casenames = fdlabelslist$casenames + varnames = fdlabelslist$varnames + if (is.null(xlab)) + xlab <- xlabel + if (is.null(ylab)) + ylab <- ylabel + if (ndim < 2) { + plot(y, fdmat, type = "l", xlim = xlim, ylim = ylim, + xlab = xlab, ylab = ylab, axes = Axes + , xaxt="n", ...) + # xaxis labels and locations + x.axis <- seq(from=0, to=24, by=1) + x.axis.lab <- paste0(x.axis,":00") + x.axis.loc <- x.axis/24*288 + axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) + + if (axFun) + do.call(axList[[1]], axList[-1]) + if (zerofind(fdmat) && href) + abline(h = 0, lty = 2) + } + if (ndim == 2) { + if (!ask) { + matplot(y, fdmat, type = "l", xlim = xlim, + ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes + , xaxt="n" + #, col="blue1" + , family=family, ps=11 + , ... + ) + # xaxis labels and locations + # x.axis <- seq(from=0, to=24, by=1) + x.axis <- seq(from=0, to=24, by=4) + x.axis.lab <- paste0(x.axis,":00") + x.axis.loc <- x.axis/24*288 + axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) + + if (axFun) + do.call(axList[[1]], axList[-1]) + if (zerofind(fdmat) && href) + abline(h = 0, lty = 2) + } + else { + op <- par(ask = FALSE) + on.exit(par(op)) + cat("Multiple plots: Click in the plot to advance to the next") + for (irep in 1:nrep) { + plot(y, fdmat[, irep], type = "l", xlim = xlim, + ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes, + ...) + if (axFun) + do.call(axList[[1]], axList[-1]) + if (irep < 2) + par(ask = ask) + if (!is.null(casenames)) + title(casenames[irep]) + else title(paste("Case", irep)) + if (zerofind(ylim) && href) + abline(h = 0, lty = 2) + } + } + } + if (ndim == 3) { + if (!ask) { + for (ivar in 1:nvar) { + matplot(y, fdmat[, , ivar], type = "l", + xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, + ask = FALSE, axes = Axes, ...) + if (axFun) + do.call(axList[[1]], axList[-1]) + if (!is.null(varnames)) + title(varnames[ivar]) + else title(paste("Variable", ivar)) + if (zerofind(ylim) && href) + abline(h = 0, lty = 2) + } + } + else { + op <- par(ask = FALSE) + on.exit(par(op)) + cat("Multiple plots: Click in the plot to advance to the next") + for (irep in 1:nrep) { + for (ivar in 1:nvar) { + plot(y, fdmat[, irep, ivar], type = "l", + xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, + axes = Axes, ...) + if (axFun) + do.call(axList[[1]], axList[-1]) + if (!is.null(casenames)) + titlestr = casenames[irep] + else titlestr = paste("Case", irep) + if (!is.null(varnames)) { + titlestr = paste(titlestr, " ", varnames[ivar]) + } + else { + titlestr = paste(titlestr, " ", "Variable", + ivar) + } + title(titlestr) + if (zerofind(ylim) && href) + abline(h = 0, lty = 2) + } + } + } + } + "done" +} diff --git a/plot_fpca_new1.R b/code/plot_fpca_new1.R old mode 100755 new mode 100644 similarity index 96% rename from plot_fpca_new1.R rename to code/plot_fpca_new1.R index 6c1016c..e9e1638 --- a/plot_fpca_new1.R +++ b/code/plot_fpca_new1.R @@ -1,90 +1,90 @@ -### based on source code of plot.pca.fd -### revisions: -### 1) colored add/subtract pca components -### 2) change xaxis labels to clock time of 0:00 - 23:00 -### 3) add the argument of 'flip' to determine which components to flip the sign -### warning: depends on the results and interpretation, may not suitable for other data or analysis. - -plot.fpca.new <- function (x, nx = 128, pointplot = TRUE, harm = 0, expand = 0, - cycle = FALSE, flip=c(1,2,4) - , ...) { - pcafd <- x - - # xaxis labels and locations - x.axis <- seq(from=0, to=24, by=1) - x.axis.lab <- paste0(x.axis,":00") - x.axis.loc <- x.axis/24*288 - - if (!(inherits(pcafd, "pca.fd"))) - stop("Argument 'x' is not a pca.fd object.") - harmfd <- pcafd[[1]] - basisfd <- harmfd$basis - rangex <- basisfd$rangeval - { - if (length(nx) > 1) { - x <- nx - nx <- length(x) - } - else x <- seq(rangex[1], rangex[2], length = nx) - } - fdmat <- eval.fd(x, harmfd) - meanmat <- eval.fd(x, pcafd$meanfd) - dimfd <- dim(fdmat) - nharm <- dimfd[2] - plotsPerPg <- sum(par("mfrow")) - harm <- as.vector(harm) - if (harm[1] == 0) - harm <- (1:nharm) - if (length(dimfd) == 2) { - for (jharm in 1:length(harm)) { - #jharm=1 - if (jharm == 2) { - op <- par(ask = TRUE) - on.exit(par(op)) - } - iharm <- harm[jharm] - if (expand == 0) { - fac <- sqrt(pcafd$values[iharm]) - } - else { - fac <- expand - } - vecharm <- fdmat[, iharm] - #if (jharm == 3){ # fPCA3 - if (jharm %in% flip){ - pcmat <- cbind(meanmat + fac * (-vecharm) - , meanmat - fac * (-vecharm) - ) - } - #else{ # fPCA1,2,4 flip the sign - else{ - pcmat <- cbind(meanmat + fac * vecharm - , meanmat - fac * vecharm - ) - } - - if (pointplot) - plottype <- "p" - else plottype <- "l" - percentvar <- round(100 * pcafd$varprop[iharm], 1) - plot(x, meanmat, type = "l", ylim = c(min(pcmat), - max(pcmat)), ylab = paste("Harmonic", iharm), - main = paste("PCA function", iharm, "(Percentage of variability", - percentvar, ")") - , xaxt="n", xlab=NA , ... - ) - axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) - - if (pointplot) { - points(x, pcmat[, 1], pch = "+", col="darkgoldenrod1") - points(x, pcmat[, 2], pch = "-", col="blue1") - } - else { - lines(x, pcmat[, 1], lty = 2, col="darkgoldenrod1") - lines(x, pcmat[, 2], lty = 3, col="blue1") - } - } - } -# this is not a complete version as original plot.pca.fd function, the rest starting from here is deleted. - invisible(NULL) -} +### based on source code of plot.pca.fd +### revisions: +### 1) colored add/subtract pca components +### 2) change xaxis labels to clock time of 0:00 - 23:00 +### 3) add the argument of 'flip' to determine which components to flip the sign +### warning: depends on the results and interpretation, may not suitable for other data or analysis. + +plot.fpca.new <- function (x, nx = 128, pointplot = TRUE, harm = 0, expand = 0, + cycle = FALSE, flip=c(1,2,4) + , ...) { + pcafd <- x + + # xaxis labels and locations + x.axis <- seq(from=0, to=24, by=1) + x.axis.lab <- paste0(x.axis,":00") + x.axis.loc <- x.axis/24*288 + + if (!(inherits(pcafd, "pca.fd"))) + stop("Argument 'x' is not a pca.fd object.") + harmfd <- pcafd[[1]] + basisfd <- harmfd$basis + rangex <- basisfd$rangeval + { + if (length(nx) > 1) { + x <- nx + nx <- length(x) + } + else x <- seq(rangex[1], rangex[2], length = nx) + } + fdmat <- eval.fd(x, harmfd) + meanmat <- eval.fd(x, pcafd$meanfd) + dimfd <- dim(fdmat) + nharm <- dimfd[2] + plotsPerPg <- sum(par("mfrow")) + harm <- as.vector(harm) + if (harm[1] == 0) + harm <- (1:nharm) + if (length(dimfd) == 2) { + for (jharm in 1:length(harm)) { + #jharm=1 + if (jharm == 2) { + op <- par(ask = TRUE) + on.exit(par(op)) + } + iharm <- harm[jharm] + if (expand == 0) { + fac <- sqrt(pcafd$values[iharm]) + } + else { + fac <- expand + } + vecharm <- fdmat[, iharm] + #if (jharm == 3){ # fPCA3 + if (jharm %in% flip){ + pcmat <- cbind(meanmat + fac * (-vecharm) + , meanmat - fac * (-vecharm) + ) + } + #else{ # fPCA1,2,4 flip the sign + else{ + pcmat <- cbind(meanmat + fac * vecharm + , meanmat - fac * vecharm + ) + } + + if (pointplot) + plottype <- "p" + else plottype <- "l" + percentvar <- round(100 * pcafd$varprop[iharm], 1) + plot(x, meanmat, type = "l", ylim = c(min(pcmat), + max(pcmat)), ylab = paste("Harmonic", iharm), + main = paste("PCA function", iharm, "(Percentage of variability", + percentvar, ")") + , xaxt="n", xlab=NA , ... + ) + axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) + + if (pointplot) { + points(x, pcmat[, 1], pch = "+", col="darkgoldenrod1") + points(x, pcmat[, 2], pch = "-", col="blue1") + } + else { + lines(x, pcmat[, 1], lty = 2, col="darkgoldenrod1") + lines(x, pcmat[, 2], lty = 3, col="blue1") + } + } + } +# this is not a complete version as original plot.pca.fd function, the rest starting from here is deleted. + invisible(NULL) +} diff --git a/NHANES_sample_data.rds b/data/NHANES_sample_data.rds similarity index 100% rename from NHANES_sample_data.rds rename to data/NHANES_sample_data.rds diff --git a/NHANES_sample_data_day_average.rds b/data/NHANES_sample_data_day_average.rds similarity index 100% rename from NHANES_sample_data_day_average.rds rename to data/NHANES_sample_data_day_average.rds From b66c4313fd7e759dc2b20ff6e5fd28d7f1948c33 Mon Sep 17 00:00:00 2001 From: Jiachen Lu Date: Thu, 15 Sep 2022 14:28:49 -0500 Subject: [PATCH 2/2] add fPCs figure to Readme --- README.md | 11 ++++++----- image/fPCs.png | Bin 0 -> 11921 bytes 2 files changed, 6 insertions(+), 5 deletions(-) create mode 100644 image/fPCs.png diff --git a/README.md b/README.md index c07ca8a..89a7421 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,8 @@ This study aimed at characterizing the rest-activity patterns in the general pop ## Findings -* Four distinct rest-activity profiles were identified to describe aspects of overall amplitude, early rising time, prolonged daytime activity and biphasic pattern. [insert the profiles figure here. need access to upload figure] +* Four distinct rest-activity profiles were identified to describe aspects of overall amplitude, early rising time, prolonged daytime activity and biphasic pattern. +![fPCs](image/fPCs.png) * Rest-activity profiles are highly associated with age, race, education levels and household income levels. * Rest-activity profiles differ on weekdays and weekends by demographics and socioeconomics. * Lower overall activity is associated with higher odds of self-reported poor or fair health. @@ -25,7 +26,7 @@ This study aimed at characterizing the rest-activity patterns in the general pop In a nationally representative sample of US adults, we identified four distinct profiles for the 24-h rest-activity cycle. We found considerable variation in these profiles across different subgroups by age, gender, race/ethnicity, SES and work status. We also observed associations between rest-activity profiles and self-rated health status. ## Data -A sample data of 50 randomly selected subjects from NHANES 2011-2014 with accelerometer data can be found [here](NHANES_sample_data.rds). +A sample data of 50 randomly selected subjects from NHANES 2011-2014 with accelerometer data can be found [here](data/NHANES_sample_data.rds). This dataset includes the following variables: - SEQN @@ -36,13 +37,13 @@ This dataset includes the following variables: ## Codes -The main code to run fpca analysis on the sample data can be found [here](fPCA_analysis.R) +The main code to run fpca analysis on the sample data can be found [here](code/fPCA_analysis.R) -R function to [plot a functional object](plot_fd_new1.R) +R function to [plot a functional object](code/plot_fd_new1.R) - code adapted from plot.fd{fda} - changes are made to present labels of xasis as clock time -R function to generate Figure 1 in the published paper of [functional principal components](plot_fpca_new1.R) +R function to generate Figure 1 in the published paper of [functional principal components](code/plot_fpca_new1.R) - code adapted from plot.pca.fd{fda} - add argument 'flip' to decide which component(s) to have sign(s) flipped ('+' to '-' or vice versa) - use different colors to represent component's sign diff --git a/image/fPCs.png b/image/fPCs.png new file mode 100644 index 0000000000000000000000000000000000000000..d9e846375baddb6949b5d395bc2862d289f3eca4 GIT binary patch literal 11921 zcmeI&XH-+cz9{gW5PI(&LN8K6uaZ!t3J547qJSVxdXa8I5m7)uRJw_Rf}jE-2uKNn zfKn9@qzNd!3xpDq{o=XzeR${GSI(!m)>|)YvUVoxKQlW!lgaGg%-wS~X3Pu%3;+Ny zTbP^J0{|7}C1jTF??Ed0n3s})pF4a03?&Nyx&Yt>0FwX!0`fNiU0p!e3((C1bV0z9 z5wP?J98&yu1LfEWoQC@B%4iC~>-g*FVC&biJ~4z3O#+n01p`bhGMpvnDAu zQfi?T@cJ8+MLBr=ebjqpQPQMWJtcuz@<5g@CFxQMXHkOcC$lEOdM}{fOSe8tmy)Xo z!AZ->$@5c_K9n>&8S;BFVs|oec`{3PGRtc+3!JR?nyk;7tgoN^TXb>~1i>U4FoOns z#7wCOtj_`)nZZUDu#pvPWCI)d!A2plQ59@70(4 zT}mx~Tc#YAZacD4t`L2&`K2%bVCne#gmed$djbF)urM(^A62wboaQlQTHJm87l-kR zVX#kXyV7&&w|RH8)Lj~Rc0kQ&Yz&Z7Q4p2)1 z2ul|(koGECiN+VEKbVNF+uW~OzF4?AYk#x1*R=$7>ND@jFo`x--&t5ziCLo2Cls$V zgkA}Iau=Jj7znEhtGRTB$Wyadyq!1BJa=^-oFmh!9;K8})W>1LU(lhZx@lhP)* z4mamp_wLsWuI0*bYkWCAymzM3DqdY_mo4Bj?qZ9e=IG97z4O(vb^LZE$(sXDJ|REujYE?=!mhplkXE5f7W_gw-d@@KJ$J|HGq(7Aroj>x z{6$)w4e4+DxKLFN_&k)Hh>1<4+87jlB4qF>(zV8RIY2U_Ql_-(c>VXRL^zTpU_4(Y zX5O7pB57#hkHraPwcV5v#h&-JX5HhP4&E{rIyUI#+6XoGKt1=Qo)$9L53=Gt`O&uc z{Oo74x@|Mh1fv;oPz&lGMGfzNwbey2MRnxc!r~dghav zv52CoIcXkX?kOMN3nb8gVD7Fvb4C5$hnTqcU0<`SZun_)A=gg%3_quwi!A#w@BNiC zb$%&saUC@1vmZ*Z!zz*1wJ8Ab^jidVM1!Gc^>G2*bkM%|E0!vy+`4ql=?rPBj{ zv6w--LmKz?I`$VbEC0vE-&b!+kaNr~yE~TZ;Ua5o#{)W>m61lRoMfKq!PYmv!Cc5n zApCZyY0`u7)c1G6YopUIF{C3uW2P1za*@*s{YF%Ul((0iF{1@LtSNB_32c6Dw2o6Aw ztQ;Zx$7a0V8W9rPY^^Uhvu(J!Z(mmsgGALLm9Og(%A)7%hR&N6gX_<>#Wg=4~{+r#iAq+oTF z!R~iOkm%9^%^Of%7E%}^=hL}?2}xphOL5zRw*0d*+p;VE$^Gr28`!%5;3o3F{ZYSQ z#6h{D%EhZe*xq8{RY;D96G9AzKyw_;TLU*LWPU_TmIKsqOv51uzCf79#mZ)C!2$9< zs0oBh(nC|zWFYF^kVquI$miZ?Q1jR=MmUBI1%yB#Y)*QjocQI7!wxi%&QFDmoi5nA z-;JrAr!ZQJh03lGWU;HHUBNowMH&E`#qUZaY<&e!#uep*a%Qbx23eDr>#JAcWHk)Lw+ZaEsBumb8zx1>y!E)gU$-=xjxI7LM5@ z?(`+F!U8eCZ}q^|3m@J`CO5x<2C!niBrvyj$6$_3EI~3B01VODpYOX_*EGDPx=c6> zil)O?998}@oN%*?U*n*n!}UJw+DcryK4dKweRM9AdnkZx@GQ)gof`7)(D&q~2d2#e z;eh%BWoztOS4R2T@(^(M$k@&8IM<+VARPK(|-(R_0-I zgb)=9D0hA$EZh)5}>Q7&!a3BO=j}?IYV&hWZ zh$Jn@)RipP8|iJFRXX}-+CD0hhD%CB&0B*-DvM1MZZo1HCj#!P_wrx}d+?b%(W3!m zQep9doUmn+3j`J(@4Q|Nw}omWib z^1}xvaXO$M;QM;Nq)?yFyD(>RbBEiQ-Rb#TZ3rEY8bg)uLW@}YZ|06tskkF$W0{`Z z5YMgBy`p4}BzoSOoP`>6$bz|x&FZGe9w!GlB)Pns-NRKsuP8cHw=&Odzk+sGq2cpF z?SE(HaodahW<2hS+z*M6pKM9WK6q^LD!`#q#Ark=;M1hSNtBD|is^zYQ1L!6%JQ0k zzy~#=_Gy9?)^*}AucS{`Rw>MStr-9Bj}T0UDKuA4v7Orj9lN7P9JyU3oYw) z>M1ru&j^ro=jo|Nv7N#}D5WEMm`&O86cynX=G4j zMNb#acL8@BxXG%-oyL|neJ`^8BctlzqY3iW2f;?t7`lGf?=Ysf?n*$@*ShT2VKz1$ zRMNCHi3yx)MYMuTr>u!qRzI|E0<@{~v#@+dOrYtDMsC*|Ft#|<;xItsb1SMf*-)1A zNZ-UFgs9@84z20t7px+iq{~#e>!`RoM6VA6nNaf)@^Ow&`!Txg66Q9reRt|hJGpoI z-o*rm{ZknA4Wz8cJ2KX71E%OH3O~MDDtW#366$*;>H$Yr07$SI;(+1C@7;A0j#C_w zbix$%-F^eIE(3O5#c;V4C`t`cKlW(WU7_T4-gcgjp;z?T5pq(?{RSSHbUj0RMJ0Eb zx%a6Txm0LqUO{I5G9&yJh}`?4{@k?!ST8^m!+u{qm3|*`j4*H-&XIM;ezx%=Pr5C@ zF44Ia&akCt2Yt%>`I8thu*A}4Yw+b@f7oreNPc55*X+8Cz1LfCX!|oGN)>*|Ke@O{ zNR2D*)Jy|Wl;rMA#a|8gpW>sx>{|@q%|JXN{gB9U7qu8w(n*QL&^cN?t=DT;O6I|w z*}XWHd7j3)BR5Vzb)XD5Pp~xv?(Tt)t)B8IWvz#}ND~&__m3OFr(5+bm ze#Edbvw8<1Q{d$u$j9hA(SVX0fVc>M&G++Wdwo-_?<2MS1e!C7V@laco`h__CZ|rk z052n>y9Qp0Mm)bnTkuso44G40N|ju5iQ!5CaQ0h4wEyOrs*uFD)3oMekdJ&?Y3r;V zzN=yk#a|AjF`er_bBbS4nIN7} zk1HkJKb^(yE|GF=47}qkCVL%QB5jduA=; z*+SJhX(Q38l1 z-=QO|Yw49@_Pe}Ym>T#Ko1}{EK=k#ooD|B*UlJcYQujIT2p$%q^etqW;iEo1XsO@P zgufbmX&7g8^k5yNUS|FMWjjpZ{M~uUbcwFVnQ5*NSzC#epQ$TM9bdS_wa^7|8kYdM z;;w+pj-h%0)U*}|3d7JTcWbQPeo)T-NnN{J2PMM7W{PCo)}&3XX7_+&BnyzAJIB?v z+jKDHi=C5g37ikNgxuozY#!8HLJjGlHeOowv_r&~S$iiNP%DSHJeOlM)0dvupx+*t z+l&wIk#<*Nn7kmGO6PP9%D*Ct&hC1jo^~l0bQrWMi(;;qhn>dw&adQAf7eS{*gh@@lq9-hg1H^&468rLTaZXYD+v>U`g$a zdX-JBEHEn%_&yo7SJTRlG%Vb_#&cplx({G=H zIuG&);*nK^NX-j3)qr$S%D?%m&0CTA7{F8rH&Tdn*QVhU?7l-P8|%#5%~!8DxO2G= zLKzd^|G9I^+G880XzQ4SpuOJmp35&#NR)te;eIPMG!lT^$zR`dOz?#Ez#i_;@Uh+G zE^i62$`F-}akVB5(0Pwb&8(K#?w0%j-jS(H*KP6U4#$NyF>+Gk9$(EEpZ}ftN&C8#OP9FGS znkD>=A7#UIEdGy~h~HfbF=p3Kg)Iy zCxM_@;JW9DolH*-(B~+&J&XM6Q*7B&4uvhUV-gKH?oXe4qXwkJF>SwN@T_KNS3kVG zE7ZIEVd&`TQ2@QplZ$gxOyjljfzyMMkY&W>+c=(ZEukBFXN?;#psZ%!H(5QXcFkLx zFfcYP!s$>WQr-XX6&%6)?kWR!JY(wzYlztykXD%6TLweD?5ftAqEyL-+T9GZl5l{i z0!J@}*+BB%7!6@Jf;A@+IFf;ly0|a&!PDFGGF)%&O9;vz%+?Mo2yZUOl*)vMU#o zI8TRvE#qn`wU;Uj6pCD_j-Q#%B*%Y^*Js4BKV9HMwG4qI%>#yIPJ19pBkRjq=8{kZz*WJ=XkB~m3Lln_%r)3l3OJZIc zrF4oyk_URCnBYm<;xRZ82sJ+Q8PnsBqa3KygQ3T0qj1Q`V-JA4&y~jkBF=6Kj3(q6 zo{~Au+T@|rkEDlGtIQT^d3oTCWbk>_L_jNJ+AJx&GU@zCnVA`*ANk7seZZ4x$O5zH zT$Up$WAm^Y+7Q3%s07)WW*@vwAjh=_VA-sIk9@csJo82p*(228f#jhLVk@RA@)M$T z&SOS6y12(aRT^_lsDmf$6xJL#WVDR>K$3Y;9+7pho}c#574Ahrm-o(;TyCDCy|OMV zG^dbd4}DE%=|j$@hMzZ)I!Mae1-vYVVtJsBE3*aOG zb|Wb=DKs_Pm0n_d8;k9Gel`v=q_p%xPbBx-NjyEKxqzbteY%NGZND;*RLViKIFjkvMGFwJ6Cg%AT|ySjQo+7|4vD^ zkc!JYHeCjy&C)JxxrZWICu%vE%_fjeUsniOefrvnugvi4Czj;bKH0C^Z0F5NjGrrIS(7YIadrq zWL2)i0E;b5^yBELm((!C)^&3{xx-2kj)_(_h7v@t`avLlvHTp>;I=&@o!iWRRwX{l z4U|IuXYlqfGUq@Y9BC(Pu;IdhW%7MeF-9{}e4aMf%iW?>Hm>g59U1qpS|=+5JJkwC zB~Z<~g}H;AZfLPrX{0g+xA?aBYfr8+V9gVE*O#GNgBSNxsJ|yqV@rayPacoviR;s+ zVK;uB9!g$dFW{=So@$6|1&CwA>qtA`@V4C9h;{oz&mO{ub-@6)6S-KyH0kvhF{uXA zBCcYMRqAuTa&7%#Uz*9$41*b8{r z?$pU)W$1}ky(K3e)44hFUu((P&+0;@Lt z3^wO|?~1e~m)M%Adm-73Qv{s|4!V6b6<0q6r@wLGZ|G!oAc3`+Y2vpjFGVl=!ujtD z>>E~wU5I)BQMZ~{o!@*t%zwL17ORXW@2eASxq%1MVu~DXo?v!IZxw5w=C!jQ=gr5? zlA5+9siGvjp5HAr^RdTA4DaXt{rbR|0@9o;>|as#U;4&BtL^{80gLQ-a&tum_UdRP zy1XJbcR6YL^A0!3&q3JbOJMZT4K5#zlM{=GZ9##he>G`cyP4AzkS^}q5HauNV$LK~ zZ7D6kUKvLxh1wNId~BRVD?y5$HW#x`1_av`LVsSF%AHXm93gES1!LC!Q=}D?(*7aR z#e#n&x{gDw%vME#yq+2K{O>GI+u@zIvaR{)vo&BGbR7HhXK4P`AUeus{br=yg>JO* zb!Sn7>0~mqV3yDNRcw< z0eT8Y%bP@(K=a(dOJ*46DqzI@jW17)<9|~0W%2)@=#hn2;}`d?`sv1Trt$EX-U`9c zczT*p95s15<=2FW>)cJ})2)|(E%%IR)Ct|g?(@Bq8pcU#KBOZtcJDoN)-Xo2x+fu} z3Aa9)xjku&x%|2Lj?ZB&IQ?nxCF8>fq#KTT(^XdyA3jW?yL9j+FPMMLg-x0*X)Y?WrI;;f~hZe=asH_GC`Pa*WJcE;j+q`Pc#&2y9{TiJ1ewPs8g zvsLxzFs7!b{tb(kjDL9Wf?@njD0JHsz5a?VQc8-*NQ$3_&gLMlVd~s`JQH3(>r1!i z+0L3?n4!uBBJ}W$0M&MKpvHFM;=PFMFq94a3ge5s5y;*U72iUIHm^!%tN86OmuLi|1kqxej-YL z5)~Q_$0!N!9wVx0*f7k4u@QgbVwe*E=i|jk@$-Ug7@fC^2A804psQS}PP`%B7}Sws zgd#lZ-_XI%6*@(p7b%|(gJ@ge&OzR3s?`)V6=F5%FaQ&?Ec2Pwfo5?jE;$PQUYs0A zDyhOnLGlVC**L}~`g5y#ke!w2TOP^@-(ZMk)e#QFZrls=(qPDV{u|z7f=1wh>#co_ zismsF=wp6aE8tym+x6nEO2ehE!}#WS6fVT2DF)Dn1%%@)$6TIA-f8Sbv{U){$ypXzc`pm@gW!|MW3 zuL3c8^tdKy&`@Yb8e#4H;xGPInCc!I#;?li6?EOj((z#g>w~w*=^{89IzG>fq1>sZ z27q=`$1INdq`S_eM-%$26r!{^Kz&ip0uzDH zG%k4rQ&nbO^p?^6brQveDmXoK`NLvn$mE4^h3l7~=QKdCw5uWLvdC4|Mg>bfnI0_g z?G<3WREEw?gUQSMy%QT?#Ri(wN*;+lOok2N|BUF}dnOzoJNz>N7KPYVahrqkZGn3I z=<(1iniqrT+X~C_VpX3> zAMHot1l3DIP81xXOSZ_VWxP&>TjU=v-)HXzepA7O$GVbxdnjbU@bDU|`X;l0=jX4? zZ=W~BV@$I)_+JIVczw3#@Ot)sw7S^O4O|s)Z#7mqOxZC;W%1%(BULKJ6k$_c(Zc8!?-&%kma7=X@n9`Max*hbzoV?D z08Og$AG_bD`drAJB&nhtu-jIH$oiyeQl`AC48)Bj6dhy-gj?Ho=xG_#`Ye<40;g`q zC_dfE)OYieff)ZSupDAGy@yD ziiQPC<4BXwA-d=0HB?AqTCS;&{B#+engD%l2g2_}-Tkw^i@?C;nQPL9zEdU!!vW3c z-BF<+K3&6JPLwyTzvN23E7ABFQseT6s}F&us+rWGlTY0l7WWyrP-P%%Vd$NTHyNh) z;5zh%pwDg{0>UzJpw{OFvEsq!5jmIZf*yUXjZR}n1(ivD=bs(mb3-szv*9<}r7oq{ zOP9(`nDT-;Szl~fJ1r*3SCF)V$f-4{ZUjG7A2mueI{eaImR6>VrwY(`zn5p93~FlH zOc1u7Y~7tXKrl7I!ys81vwM+=z6l;T0ngO1$5~OGVYlwm{5kd76QKJPjz=x*m z$!GJ!1f~hjc(F(=@zJ|7kkSCqq+W!&scKf<3q~0|ennU`ZB(zD0XZcyDo3p`JoV}j zmtKQR9q~0U$PY%GL(_;5&RLJhPPnV zNLQUjso6_h>Lb}a=FOn3L^?Aldt&&7bLOWAvREM!TIN@O?m02a`@g!2lJLcXudR49a zU&tA{iEU2|`k4{X+-uLBY2Stny2$JhQ#pg?wQML%@On z(%59Z?+CV&$0RWv+oYFmp0(39R1--l1P}P_c>yyO=K5%O7ePC4NVG@$N)xcY5mcNk zEZnBFLZ}K&e7ue73p?K3{M-+~y1XI$NH(=8h*c;JW=9P3I|iQZNkq@(iv7MY>kUi& zLuB#|ghcv!L53+yT80v26iqC>?K9|Xx~^gCBEWtND>kQyfs{}T`Oe%kGODY*HK#oZ z7)d6_(|j$nnudH#*D;)dv8~XmR(}n+#iJ?UQH$}nOjpLbJQZLta?+pw(l1MOvchB@6S?PAE&g%)D`RVVuFFRVt@`UipPw&9x? zXK8vjqqjL8St$P}X;h_o)Siu;Q@p6-^a1FfCse=1F)4GYj~0Y5Aj92*)Fm%%#%eyu z_neDjaayB0#$7K}`)Z;k3IfqcPdzF;QxOgq2B57U>7Xd*ibzs9HLQc-3ldQS>irDd zVZmN+*m0}PkKh3Xrr$)9jRzq5Kjho^QN3=$N<~2P!S@Tr*?_=tEhdgS%YzuXf|;mA ziMYVA{TjBAy+|d60qFDj#KmQx?))3dZC6!&_rpUs}IBL+}S=hS1ll42NKdo*?<~qpZZ{!wF z^Pr5MA6m7C6O#bhN56#J4rW77q4x`YnGA{vanwU~jnT%GPnrqs38TS>#;!^O-ujI2 zUVgyTi6MZ#fx*O@FgnKh>`N34+|3*O4a3j()qzzd^cWVh!|TT%D5F$qZUz1euWPVejLyzleRjs#MiNi88~|D*f~vu&f%Y5s7U(x-VQZvSGA7 zsC0IG2#-`B@YvW1@~*qxMR4b%5Y8MLJkIR)Ycr^2L7x3tI#9Q3qA+R9Rzg`H+;)0c z%+0$3oZ-Y+sdag|7Ep2BfCckW@lf7kLk-7jTPGU=!#Tl1RMpeDW8{)XR3PWvCmqy| zDFtIEt}(HB16HA!XLK0;EJsX;|<|IUBF!AA9kB@-lk zJ;-aY3*zhuihiXIrR-N&20jTv#+B~~jOYA_53$1KnM13Esd^ZOphVCKS}jDyvj&_I zhA4Dn&Je&5GvJqsTMf%O69}8bOsfi(c0B48Kf4vM=z;G~6vmZH0B@k*E72u^-?d9$ zwiO6c3goC@S}m!~-ZqnwZyC>ftUVD>yrc9E%#;HRXvz<$dM>6-K}$@yF-LsWBV;6* zZfKyhRKP7Lhlq{IC!zP)0sj4#u?hybAq(#$*v*ZKw6gcQ>PPUaA)&(|Y~_OLO`l_9 zTKe`yEl=2S@xIG=^2oNUTZeKjXo~{+`B3t1=co_YwDzwN!ExxU2fh&dS>Mj@=84k6 zb4?Hi!%rMiD5B=M!|qRo1~5|*XlK4gOi$c*-6f+N|H?XozL{GJPd=@u8x!dnewBeD z((64v>$dU3XL6DOAE{~)v+~5@w#ub$D|jn#C_7Y(bS-oWJ#dwgz!cchT-~umWq=M- zAZHsF3=zgbMZaC?-5W@jqsK%Wy?25giUyKfblN#AjDd)s#|=l|!m`gwKT;8rWeFeM zPA5g~^drR}h_&|wnYExa8AS-ZDFaOjkb|(95CR;+nG