text:
# QCBS R Workshops

This series of 10 workshops walks participants through the steps required to use R for a wide array of statistical analyses relevant to research in biology and ecology. These open-access workshops were created by members of the QCBS both for members of the QCBS and the larger community.

The content of this workshop has been peer-reviewed by several QCBS members. If you would like to suggest modifications, please contact the current series coordinators, listed on the main wiki page

# Workshop 10: Advanced multivariate analyses

Developed by: Monica Granados, Emmanuelle Chrétien, Bérenger Bourgeois, Amanda Winegardner and Xavier Giroux-Bougard. (Material in R script obtained from: Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.)

Summary: In this workshop, you will learn about advanced multivariate analyses for your own data. This workshop concentrates on constrained methods such as redundancy analysis, multivariate regression tree and linear discriminant analysis to explore how environmental variables may be driving patterns in species assemblage across sites.

Make sure to load the following packages (see how in the R script):

| Load the required packages and functions
install.packages("vegan")
install.packages("mvpart")
install.packages("labdsv")
install.packages("plyr")
install.packages("MASS")

# For the two following packages, upload the file provided on the wiki page.
# To do so, go to Packages tab on the bottom right panel of R Studio
# Click on Install Packages
# Choose to install from Package Archive file and upload these two files
install.packages("MVPARTwrap")
install.packages("rdaTest")

library(vegan)
library(mvpart)
library(MVPARTwrap)
library(rdaTest)
library(labdsv)
library(plyr)
library(MASS)

# Introduction to advanced multivariate analyses

The previous workshop presented the basics of multivariate analyses:

• how to choose appropriate distance metrics and transformations
• hierarchical clustering
• unconstrained ordinations
• Principal component analysis
• Principal coordinate Analysis
• Correspondence analysis
• Nonmetric multidimensional scaling

The present workshop builds on this knowledge, and will focus on constrained analyses. All the methods overviewed during the introductory workshop allowed to find patterns in the community composition data or in the descriptors, but not to explore how environmental variables could be driving these patterns. With constrained analyses, such as redundancy analysis (RDA), linear discriminant analysis (LDA) and multivariate regression tree (MRT), one can describe and predict relationships between community composition data and environmental variables.

# Getting started with data

We will continue to use the Doubs river datasets for this workshop. “DoubsSpe.csv” is a data frame of fish community data where the first column contains site names from 1 to 30 and the remaining columns are fish taxa. The taxa columns are populated by fish abundance data (counts). “DoubsEnv.csv” is a data frame of environmental data for the same sites contained in the fish community data frame. Again, the first column contains site names from 1 to 30. The remaining columns contain measurements for 11 abiotic variables. Note that data used in ordination analyses is generally in wide-format.

| Load the Doubs Species and Environmental Data
#Species community data frame (fish abundance): “DoubsSpe.csv”
spe<- spe[-8,] #Site number 8 contains no species and so row 8 (site 8) is removed. Be careful to
#only run this command line once as you are overwriting "spe" each time.

#Environmental data frame: “DoubsEnv.csv”
env<- env[-8,] #Remove corresponding abiotic data for site 8 (since removed from fish data).
#Again, be careful to only run the last line once.

# 1. Explore and prepare the data

We can use summary functions to explore the “spe” data (fish community data) and discover things like the dimensions of the matrix, column headings and summary statistics for the columns. This is a review from Workshop 2.

| Explore DoubsSpe
names(spe) #see names of columns in spe
dim(spe) #dimensions of spe; number of columns and rows
str(spe) #displays internal structure of objects
head(spe) #first few rows of the data frame
summary(spe) #summary statistics for each column; min value, median value, max value, mean value etc.

Look at the species’ distribution frequencies.

| Species distribution of DoubsSpe data
#Species distribution
(ab <- table(unlist(spe))) #note that when you put an entire line of code in brackets like this, the output for that operation is displayed right away in the R console

barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5))

Can see that there is a high frequency of zeros in the abundance data.

See how many absences there are in the fish community data.

| Absences in fish data
sum(spe==0)

Look at the proportion of zeros in the fish community data.

| Proportion of zeros
sum(spe==0)/(nrow(spe)*ncol(spe))

The proportion of zeros in the dataset is ~0.5.

This is high but not uncommon for species abundance data. However, in order to avoid the use of double-zeros as indications of resemblance among sites, we will apply a transformation to the species data. Legendre and Gallagher (2001) proposed five pre-transformations of the species data, four of them being available in vegan in the function decostand().

The Hellinger transformation will be applied to the fish data. It expresses abundances as the square-root of their relative abundance at each site (Borcard et al. 2011).

| Hellinger transformation
spe.hel <- decostand(spe, method="hellinger") # you can also use method="hell"

Explore the environmental data and create a panel of plots to compare collinearity:

| Explore Doubs Env data
names(env)
dim(env)
str(env)
summary(env)
pairs(env, main="Bivariate Plots of the Environmental Data" )

In this case, the environmental data (explanatory variables) are all in different units and need to be standardized prior to computing distance measures to perform most ordination analyses. Standardize the environmental data (11 variables) using the function decostand() in vegan.

| Decostand
env.z <- decostand(env, method="standardize")
apply(env.z, 2, mean) # the data are now centered (means~0)
apply(env.z, 2, sd)   # the data are now scaled (standard deviations=1)

# 2. Canonical analyses

Described first by Rao (1964), canonical analysis is a generic term that for several types of statistical analyses sharing a common goal; to identify the relationship between a multivariate response table (matrix Y, generally describing the species composition of communities) and a multivariate explanatory table (matrix X, generally containing environmental descriptors) by combining ordination and regression concepts. Canonical analyses allow users to test ecological hypothesis concerning the environmental drivers of species composition. Among the diversity of canonical analysis, we will mainly focus here on Redundancy Analysis (RDA).

Redundancy Analysis is a direct extension of multiple regression, as it models the effect of an explanatory matrix X (n x p) on a response matrix Y (n x m). This is done by preforming an ordination of Y to obtain ordination axes that are linear combinations of the variables in X. In RDA, ordination axes are calculating from a PCA of a matrix Yfit, computed by fitting the Y variables to X by multivariate linear regression. Note that the explanatory variables in X can be quantitative, qualitative or binary variables. Prior to RDA, explanatory variables in Y must be centered, standardized (if explanatory variables are not dimensionally homogeneous, i.e. in different units), transformed (to limit the skew of explanatory variables) or normalized (to linearize relationships) following the same principles as in PCA. Collinearity between the X variables should also be reduced before RDA.

In order to obtain the best model of RDA, explanatory variables can be selected by forward, backward or stepwise selection that remove non-significant explanatory variables.

RDA involves two computational steps. In the first step, a matrix of fitted values Yfit is calculated trough the linear equation: Yfit = X[X'X]-1 [X'Y]

In the second step, a PCA of the fitted matrix Yfit is calculated (see equations in the figure below) and produce the canonical eigenvalues and eigenvectors together with the matrix Z containing the canonical axes. These canonical axes correspond to linear combinations of the explanatory variables in X. This linearity of the combinations of the X variables is a fundamental property of RDA. In the analysis of community composition, these canonical axes are interpreted as complex environmental gradients.

Redundancy analysis as a two-step process (from Legendre and Legendre 2012) Several statistics can be computed from RDA:

- The R² measures the strength of the canonical relationship between Y and X by calculating the proportion of the variation of Y explained by the variables in X, - The adjusted R² also measures the strength of the relationship between Y and X, but applies a correction of the R² to take into account the number of explanatory variables, - The F-statistic corresponds to an overall test of significance of an RDA by comparing the computed model to a null model. This test is based on the null hypothesis that the strength of the linear relationship calculated by the R² is not larger than the value that would be obtained for unrelated Y and X matrices of the same size. Note that F-statistics can also be used to sequentially test the significance of each canonical axis.

In R, RDA can be computed using the function rda from package vegan, as follows:

| rda() in vegan for RDA
#Preparing the data prior to RDA
env.z <- subset(env.z, select = -das) # remove the "distance from the source" variable

#Running the RDA
?rda
spe.rda <- rda(spe.hel~., data=env.z)

### Extract the results
summary(spe.rda, display=NULL)

#The results are called using summary:
summary(spe.rda, display=NULL) #display = NULL optional

The summary of the output looks like this:

These results contain the proportion of variance of Y explained by the X variables (constrained proportion, 73.41% here), the unexplained variance of Y (unconstrained proportion, 26.59% here) and then summarize the eigenvalues, the proportions explained and the cumulative proportion of each canonical axis (each canonical axis = each constraining variable, in this case, the environmental variables from env).

To select the significant explanatory variables, you can then perform a forward selection (or backwards or stepwise), using the ordiR2step() function (or using the forward.sel function of package packfor):

| ordiR2step() for forward selection
### Select the significant explanatory variables by forward selection
?ordiR2step
ordiR2step(rda(spe.hel~1, data=env.z), scope= formula(spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
env.signif <- subset(env.z, select = c("alt", "oxy", "dbo"))

The resulting output reads: rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z) with the proportion of variation explained by the three constraining variables being 0.59. So in this case, only three variables are retained by the forward selection, i.e. alt, oxy and dbo. These three variables can be placed in a new data frame env.signif to perform the new RDA including only significant X variables:

| RDA with the two retained variables
spe.rda.signif <- rda(spe.hel~., data=env.signif)
summary(spe.rda.signif, display=NULL)

The explanatory variables (altitude, oxygen and biological oxygen demand) now explain 59% of the variance in Y (species). The adjusted R² of this RDA is calculated using the function RsquareAdj:

(R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared) #Here the strength of the relationship between X and Y corrected for the number of X variables is 0.54. The significance of the model and of each canonical axis can be tested using the function anova (note this is different from retaining significant variables as was done with forward selection, now we're testing the significance of the RDA axes): | anova.cca for testing the significance of axes ?anova.cca anova.cca(spe.rda.signif, step=1000) anova.cca(spe.rda.signif, step=1000, by="axis") #In this case, the RDA model is highly significant (p=0.001) as well as all three canonical axes. To visualize the results of the RDA, triplots can be drawn using the plot(). Note that as in PCA, users can create scaling 1 and scaling 2 triplots. In scaling 1, distance among objects approximate their Euclidean distances while in scaling 2, angles between variables X and Y reflect their correlation. Thus, scaling 1 triplots can be used to interpret distances among objects and scaling 2 triplots to interpret the relationships between X and Y. To plot the scaling 1 triplots of the RDA, the following code can be used: | plot RDAs #Quick plots scaling 1 and 2 windows() plot(spe.rda.signif, scaling=1, main="Triplot RDA (scaling 1)") windows() plot(spe.rda.signif, scaling=2, main="Triplot RDA (scaling 2)") #Advanced plots scaling 1 windows() plot(spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1)) points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1), pch=21, col="black", bg="steelblue", cex=1.2) arrows(0,0, scores(spe.rda.signif, display="species", choices=c(1), scaling=1), scores(spe.rda.signif, display="species", choices=c(2), scaling=1), col="black",length=0) text(scores(spe.rda.signif, display="species", choices=c(1), scaling=1), scores(spe.rda.signif, display="species", choices=c(2), scaling=1), labels=rownames(scores(spe.rda.signif, display="species", scaling=1)), col="black", cex=0.8) arrows(0,0, scores(spe.rda.signif, display="bp", choices=c(1), scaling=1), scores(spe.rda.signif, display="bp", choices=c(2), scaling=1), col="red") text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05, scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05, labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)), col="red", cex=1) #Advanced plots scaling 2 windows() plot(spe.rda.signif, scaling=2, main="Triplot RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1)) points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=2), pch=21, col="black", bg="steelblue", cex=1.2) arrows(0,0, scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2, scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2, col="black",length=0) text(scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2.1, scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2.1, labels=rownames(scores(spe.rda.signif, display="species", scaling=2)), col="black", cex=0.8) arrows(0,0, scores(spe.rda.signif, display="bp", choices=c(1), scaling=2), scores(spe.rda.signif, display="bp", choices=c(2), scaling=2), col="red") text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=2)+0.05, scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)+0.05, labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)), col="red", cex=1)  And the final triplots would look like this: Challenge 1: Run an RDA of the mite environmental variables constraining the mite species abundances. Use: | Load the mite data #Load the mite species and environmental data from vegan package data(mite) mite.spe<-mite mite.spe.hel <- decostand(mite.spe, method="hellinger") data(mite.env) What are the significant explanatory variables? How much of the variation in species data is explain by significant explanatory variables (i.e. those that were selected)? What are the significant axes? What group(s) of sites can you identify? What species are related to each group(s) of sites? Challenge 1: Solution Click to display ⇲ Click to hide ⇱ Click to hide ⇱ Your code probably looks something like this: | RDA of mite data #Initial RDA with ALL of the environmental data mite.spe.rda<-rda(mite.spe.hel~., data=mite.env) #Select significant environmental variables ordiR2step(rda(mite.spe.hel~1, data=mite.env), scope= formula(mite.spe.rda), direction= "forward", R2scope=TRUE, pstep=1000) #Create a new dataframe with only the significant variables that you identified above mite.env.signif <- subset(mite.env, select = c("WatrCont", "Shrub", "Substrate", "Topo", "SubsDens")) #Re-run the RDA with the significant variables and look at the summary mite.spe.rda.signif=rda(mite.spe~., data=mite.env.signif) summary(mite.spe.rda.signif, display=NULL) #Find the R2 adjusted of the model with the retained environmental variables (R2adj <- RsquareAdj(mite.spe.rda.signif)$adj.r.squared)

#Determine the significant canonical (constrained) axes)
anova.cca(mite.spe.rda.signif, step=1000)
anova.cca(mite.spe.rda.signif, step=1000, by="axis")

#Plot the RDA
windows()
plot(mite.spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(mite.spe.rda.signif, display="species", choices=c(1), scaling=1),
scores(mite.spe.rda.signif, display="species", choices=c(2), scaling=1),
labels=rownames(scores(mite.spe.rda.signif, display="species", scaling=1)),
col="grey", cex=0.8)
arrows(0,0,
scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1),
scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1),
col="red")
text(scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
labels=rownames(scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)),
col="red", cex=1)

Five explanatory variables are significant: WatrCont, Shrub, Substrate, Topo, and SubsDens. The proportion of the variance explained by these variables is 32.08% and the adjusted R² of this model is 19.20%. While the first canonical is significant (p<0.001), the overall model is however not significant (p=0.107). Three groups of sites appear on the triplot. One group of sites with high water content (top right) have high abundance of species 9 and 25. Species 6, 12, 19 and 26 are related to another group of sites with hummock topography (bottom left). The last group of sites (top left) have not characteristic species and heterogeneous environmental conditions. See the triplot below:

Partial RDA is a special case of RDA in which the response variables Y are related to explanatory variables X in the presence of additional explanatory variables, W, called covariables. As in partial linear regression, the linear effect of X variables on the Y variables are adjusted for the effects of the covariables W. For this, a RDA of the covariables W on the response variables Y is first performed. The residuals of this RDA are then extracted, i.e. a matrix Yres|W containing the Y response variables in which the effect of W were removed. The partial RDA finally correspond to the RDA of X on Yres|W. All statistics previously presented for RDA also apply for partial RDA.

Partial RDA is thus a powerful tool when users what to assess the effect of environmental variables on species composition while taking into account the species variation due to other environmental variables with no interest. It can also be used to control for well-known linear effects, isolate the effect of a single explanatory variable or analyse related samples. In the example below, we will assess the effect of water chemistry on fish species abundances partialling out the effect of physiography.

In R, partial RDA is performed in the same way as RDA using rda() with the addition of a condition term:

| Partial RDA with rda()
#Divide the env2 dataframe into two dataframes:
envtopo <- env[, c(1:3)] # Physiography : explanatory dataset 1
names(envtopo)
envchem <- env[, c(4:10)] # Water quality : explanatory dataset 2
names(envchem)

#Run the partial RDA
spechem.physio=rda(spe.hel, envchem, envtopo)
summary(spechem.physio, display=NULL)
#or
spechem.physio2=rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo
+ Condition(alt + pen + deb), data=env)

#Extract the results
summary(spechem.physio, display=NULL)

#Calculate the adjusted R2 of the partial RDA
(R2adj <- RsquareAdj(spechem.physio)$adj.r.squared) #Test the significance of the axes in partial RDA anova.cca(spechem.physio, step=1000) anova.cca(spechem.physio2, step=1000, by="axis") #Construct the triplots #Scaling 1 windows(title="Partial RDA scaling 1") plot(spechem.physio, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1)) points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=1), pch=21, col="black", bg="steelblue", cex=1.2) arrows(0,0, scores(spechem.physio, display="species", choices=c(1), scaling=1), scores(spechem.physio, display="species", choices=c(2), scaling=1), col="black",length=0) text(scores(spechem.physio, display="species", choices=c(1), scaling=1), scores(spechem.physio, display="species", choices=c(2), scaling=1), labels=rownames(scores(spechem.physio, display="species", scaling=1)), col="black", cex=0.8) arrows(0,0, scores(spechem.physio, display="bp", choices=c(1), scaling=1), scores(spechem.physio, display="bp", choices=c(2), scaling=1), col="red") text(scores(spechem.physio, display="bp", choices=c(1), scaling=1)+0.05, scores(spechem.physio, display="bp", choices=c(2), scaling=1)+0.05, labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=1)), col="red", cex=1) #Scaling 2 windows(title="Partial RDA scaling 2") plot(spechem.physio, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1)) points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=2), pch=21, col="black", bg="steelblue", cex=1.2) arrows(0,0, scores(spechem.physio, display="species", choices=c(1), scaling=2), scores(spechem.physio, display="species", choices=c(2), scaling=2), col="black",length=0) text(scores(spechem.physio, display="species", choices=c(1), scaling=2), scores(spechem.physio, display="species", choices=c(2), scaling=2), labels=rownames(scores(spechem.physio, display="species", scaling=2)), col="black", cex=0.8) arrows(0,0, scores(spechem.physio, display="bp", choices=c(1), scaling=2), scores(spechem.physio, display="bp", choices=c(2), scaling=2), col="red") text(scores(spechem.physio, display="bp", choices=c(1), scaling=2)+0.05, scores(spechem.physio, display="bp", choices=c(2), scaling=2)+0.05, labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=2)), col="red", cex=1)  This RDA is significant (p<0.001) as well as the two first canonical axis. Water chemistry explained 31.89% of the variance of fish species composition, physiographic covariables explained 41.53% of this variation and the unexplained variation is 26.59%. The adjusted R² of this RDA is 24.13%. The triplot looks like: Challenge 2 Run the partial RDA of the mite environmental variables of the mite species abundances partialling out for the substrate variables (SubsDens, WaterCont and Substrate). Is the model significant? Which are the significant axes? Interpret the obtained triplot. Challenge 2: Solution Click to display ⇲ Click to hide ⇱ Click to hide ⇱ Here is some potential code: | Mite partial RDA #Partial RDA mite.spe.subs=rda(mite.spe.hel ~ Shrub + Topo + Condition(SubsDens + WatrCont + Substrate), data=mite.env) #Summary summary(mite.spe.subs, display=NULL) (R2adj <- RsquareAdj(mite.spe.subs)$adj.r.squared)

#Significant axes
anova.cca(mite.spe.subs, step=1000)
anova.cca(mite.spe.subs, step=1000, by="axis")

#Triplot scaling 1
windows(title="Partial RDA scaling 1")
plot(mite.spe.subs, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=1),
pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
col="black",length=0)
text(scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
labels=rownames(scores(mite.spe.subs, display="species", scaling=1)),
col="black", cex=0.8)
arrows(0,0,
scores(mite.spe.subs, display="bp", choices=c(1), scaling=1),
scores(mite.spe.subs, display="bp", choices=c(2), scaling=1),
col="red")
text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=1)+0.05,
scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)+0.05,
labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)),
col="red", cex=1)

#Triplot scaling 2
windows(title="Partial RDA scaling 2")
plot(mite.spe.subs, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=2),
pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
col="black",length=0)
text(scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
labels=rownames(scores(mite.spe.subs, display="species", scaling=2)),
col="black", cex=0.8)
arrows(0,0,
scores(mite.spe.subs, display="bp", choices=c(1), scaling=2),
scores(mite.spe.subs, display="bp", choices=c(2), scaling=2),
col="red")
text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=2)+0.05,
scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)+0.05,
labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)),
col="red", cex=1)

This RDA is significant (p<0.001) as well as the first canonical axis. Environmental variables explained 9.81% of the variance of mite species composition, substrate covariables explained 42.84% of this variation and the unexplained variation is 47.35%. The adjusted R² of this RDA is 8.33%.

Variation partitioning is a type of analysis that combines RDA and partial RDA to divide the variation of a response variable among two, three or four explanatory data sets. Variation partitioning are generally represented by Venn diagram in which the percentage of explained variance by each explanatory data set (or combination of data stets) is reported.

In the case of two datasets (below): - Fraction a + b +c is the explained variance by the two datasets calculated using a RDA of y by X + W. - Fraction d is the unexplained variance by the two datasets calculated using the same RDA as above. - Fraction a is the explained variance by the X data set only calculated using a partial of y by X with W as covariables. - Fraction c is the explained variance by the W data set only calculated using a partial of y by W with X as covariables. - Fraction b is calculated by subtraction, i.e. b = [a + b] + [b + c] - [a + b + c].

Venn diagram of partition of the variation of a response variable y among two sets of explanatory variables X and W (from Legendre and Legendre 2012).

Variation partitioning is thus an indicated analysis when user what to relate the abundance of species in a community to various type of environmental variables, for example abiotic vs biotic variables, large-scale versus small-scale variables, etc. In the next example, we will partition the variation of fish species composition between chemical and physiographic variables.

In R, variation partitioning is performed using the function varpart(). Venn diagrams can also be drawn using the function plot().

| varpart()
?varpart
vegandocs("partitioning.pdf")

#Variation partitioning with all explanatory variables
spe.part.all <- varpart(spe.hel, envchem, envtopo)
spe.part.all
windows(title="Variation partitioning - all variables")
plot(spe.part.all, digits=2)

The output looks like:

In this case, the chemical variables explain 24.10% of the fish species composition, the physiographic variables explain 11.20% of the fish species composition and the interaction of these two types of variables explained 23.30% of the fish species composition. Note that the varpart() function also identify the fractions that can be tested for significance using the function anova.cca().

Users can also perform variation partitioning between data sets that only contain significant environmental variables:

| varpart with significant variables
#RDA of chemistry variables
spe.chem <- rda(spe.hel~., data=envchem)

#Select significant chemistry variables
R2a.all.chem <- RsquareAdj(spe.chem)$adj.r.squared ordiR2step(rda(spe.hel~1, data=envchem), scope= formula(spe.chem), direction= "forward", R2scope=TRUE, pstep=1000) names(envchem) (envchem.pars <- envchem[, c( 4, 6, 7 )]) #RDA with other environmental variables spe.topo <- rda(spe.hel~., data=envtopo) R2a.all.topo <- RsquareAdj(spe.topo)$adj.r.squared
ordiR2step(rda(spe.hel~1, data=envtopo),
scope= formula(spe.topo), direction= "forward", R2scope=TRUE, pstep=1000)
names(envtopo)
envtopo.pars <- envtopo[, c(1,2)]

#Varpart
spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars)
windows(title="Variation partitioning - parsimonious subsets")
plot(spe.part, digits=2)

#Tests of significance
anova.cca(rda(spe.hel, envchem.pars), step=1000) # Test of fractions [a+b]
anova.cca(rda(spe.hel, envtopo.pars), step=1000) # Test of fractions [b+c]
env.pars <- cbind(envchem.pars, envtopo.pars)
anova.cca(rda(spe.hel, env.pars), step=1000) # Test of fractions [a+b+c]
anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step=1000) # Test of fraction [a]
anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step=1000) # Test of fraction [c]

Now, the chemical variables explain 25.30% of the fish species composition, the physiographic variables explain 14.20% of the fish species composition and the interaction of these two types of variables explained 19.60% of the fish species composition. All these fractions are significant (p<0.001).

Challenge 3: Perform variation partitioning of the mite species abundances with a first dataset for the significant substrate variables (SubsDens, WaterCont and Substrate) and a second dataset for the significant other variables (Shrud and Topo). What proportion of the variation are explained by each dataset? What are the significant fractions ?

Challenge 3: Solution

Click to display ⇲

Click to hide ⇱

Click to hide ⇱

This is what your code may look like:

| varpart with significant variables
str(mite.env)
(mite.subs=mite.env[,c(1,2,3)]) #First set of variables outlined in challenge
(mite.other=mite.env[,c(4,5)]) #Second set of variables outlined in challenge

#RDA for mite.subs
rda.mite.subs <- rda(mite.spe.hel~., data=mite.subs)
R2a.all.subs <- RsquareAdj(rda.mite.subs)$adj.r.squared #Forward selection for mite.subs ordiR2step(rda(mite.spe.hel~1, data=mite.subs), scope= formula(rda.mite.subs), direction= "forward", R2scope=TRUE, pstep=1000) names(mite.subs) (mite.subs.pars <- mite.subs[, c(2, 3)]) #RDA for mite.other rda.mite.other <- rda(mite.spe.hel~., data=mite.other) R2a.all.other <- RsquareAdj(rda.mite.other)$adj.r.squared

#Forward selection for mite.other
ordiR2step(rda(mite.spe.hel~1, data=mite.other),
scope= formula(rda.mite.other), direction= "forward", R2scope=TRUE, pstep=1000)
names(mite.other)
(mite.other.pars <- mite.other[, c(1,2)])

#Variation partitioning
(mite.spe.part <- varpart(mite.spe.hel, ~WatrCont+Substrate, ~Shrub+Topo,
data=mite.env))
windows(title="Variation partitioning - parsimonious subsets")
plot(mite.spe.part, digits=2)

# Tests of all testable fractions
anova.cca(rda(mite.spe.hel~ WatrCont+Substrate, data=mite.env), step=1000) # Test of fractions [a+b]
anova.cca(rda(mite.spe.hel~Shrub+Topo, data=mite.env), step=1000) # Test of fractions [b+c]
(env.pars <- cbind(mite.env[,c(2,3,4,5)]))
anova.cca(rda(mite.spe.hel~ WatrCont+Substrate+Shrub+Topo, data=env.pars), step=1000) # Test of fractions [a+b+c]
anova.cca(rda(mite.spe.hel~WatrCont+Substrate + Condition(Shrub+Topo), data=env.pars), step=1000) # Test of fraction [a]
anova.cca(rda(mite.spe.hel~Shrub+Topo+ Condition(WatrCont+Substrate ), data=env.pars), step=1000) # Test of fraction [c]

In this case, substrate variables explain 14.00% of the mite species composition, the other environmental variables explain 9.1% of the mite species composition and the interaction of these two types of vriables explained 16.90% of the mite species composition. All these fractions are significant (p<0.001).

# 3. Multivariate regression tree

Multivariate regression tree (MRT) is a constrained clustering technique. Introduced by De’ath (2002), MRTs allow the partitioning of a quantitative response matrix by a matrix of explanatory variables constraining (guiding) on where to divide the data of the response matrix. RDA and MRT are both regression techniques, the former explaining the global structure of relationships through a linear model, the latter better highlighting local structures and interactions among variables by producing a tree model.

Advantages of the MRT compared to the RDA:

• does not make assumptions about the shape of the relationships between species and environmental variables (quantitative or categorical),
• is robust in dealing with missing values
• is robust in dealing with collinearity among the explanatory variables
• is insensitive to transformations of the explanatory variables, which allows the use of raw values
• the outcome, the tree, is easy to interpret, especially to a non-scientist audience.

The MRT technique splits the data into clusters of samples similar in their species composition based on environmental value thresholds. It involves two procedures running at the same time: 1) the computation of the constrained partitioning of the data, and 2) the calculation of the relative error of the successive partitioning levels by multiple cross-validations. The function mvpart() from the package mvpart computes both the partition and the cross-validation.

A quick note on MRT terminology:

Leaf: Terminal group of sites

Node: Point where the data splits into two groups. It is characterized by a threshold value of an explanatory variable.

Branch: Each group formed by a split

1- Constrained partitioning of the data

First, the method computes all possible partitions of the sites into two groups. For each quantitative explanatory variable, the sites will be sorted in the ascending values of the variables; for categorical variables, the sites will be aggregated by levels to test all combinations of levels. The method will split the data after the first object, the second object and so on, and compute the sum of within-group sum of squared distances to the group mean (within-group SS) for the response data. The method will retain the partition into two groups minimizing the within-group SS and the threshold value/level of the explanatory variable. These steps will be repeated within the two subgroups formed previously, until all objects form their own group. In other words, when each leaf of the tree contains one object.

2- Cross-validation and pruning the tree

The mvpart function also performs a cross-validation and identifies the best predictive tree. The cross-validation procedure consists in using a subset of the objects to construct the tree, and to allocate the remaining objects to the groups. In a good predictive tree, objects are assigned to the appropriate groups. The cross-validated relative error (CVRE) is the measure of the predictive error. Without cross-validation, one would retain the number of partitions minimizing the variance not explained by the tree (i.e. the relative error: the sum of the within-group SS over all leaves divided by the overall SS of the data). This is the solution maximizing the R2 so to speak. This approach is explanatory rather than predictive.

Let’s create a multivariate regression tree on the Doubs data.

| mvpart()
?mvpart

#Prepare the data: remove “distance from source”
env <- subset(env, select = -das)

# Create the regression tree
doubs.mrt <- mvpart(as.matrix(spe.hel) ~. ,env,
legend=FALSE, margin=0.01, cp=0, xv="pick",
xval=nrow(spe.hel), xvmult=100, which=4)

At this point, you will need to select the tree who's size (number of groups) is appropriate to the aim of your study from the following graph. This step requires the argument xv=“pick”. In other words, you must prune the tree by picking the best-fit tree. Indeed, a fully resolved tree is not the desirable outcome. Instead, one is usually interested in a tree including only informative partitions/groups. It is possible to have an a-priori idea of the number of potential groups to be retained as well.

The graph shows the relative error RE (in green) and the cross-validated relative error CVRE (in blue) of trees of increasing size. The red dot indicates the solution with the smallest CVRE, and the orange dot shows the smallest tree within one standard error of CVRE. It has been suggested that instead of choosing the solution minimizing CVRE, it would be more parsimonious to opt for the smallest tree for which the CVRE is within one standard error of the tree with the lowest CVRE (Breiman et al. 1984). The green bars at the top indicate the number of times each size was chosen during the cross-validation process.

This graph is interactive, which means you will have to click on the blue point corresponding your choice of tree size. Once you do so, the corresponding multivariate regression tree will appear. If you click on the orange dot, the following tree appears.

The statistics at the bottom of the figure are: the residual error (the reciprocal of the R2 of the model, in this case 43.7%), the cross-validated error, and the standard error. This tree has only two leaves separated by one node. This node splits the data into two groups at the threshold altitude value of 361.5m.

Each leaf is characterized by a small barplot showing the abundances of the species, its number of sites and its relative error.

We can compare this tree with the 10-group solution, as suggested by the CVRE criterion, or choose a solution in between, e.g. with 4 leaves to compare.

| Compare trees
# Using the CVRE criterion
doubs.mrt.cvre <- mvpart(as.matrix(spe.hel)~., env,
legend=FALSE, margin=0.01, cp=0,xv="pick",
xval=nrow(spe.hel), xvmult=100,which=4)

# Choosing ourself the best number of partitions
doubs.mrt.4 <- mvpart(as.matrix(spe.hel)~., env,
legend=FALSE, margin=0.01, cp=0, xv="pick",
xval=nrow(spe.hel), xvmult=100,which=4)

The 10-group solution has a high EXPLANATORY power but its predictive power (indicated by the cross-validated error) is just slightly better than that of the 2-group solution. The 4-group solution seems to be a good compromise.

| MRT summary
summary(doubs.mrt)

CP stands for “complexity parameter”, which is the equivalent of the variance explained by each node. The CP at nsplit 0 is the R2 of the whole tree. The summary then outlines, for each node, the best threshold values to split the data. While informative, this output is very dense. A more detailed and yet more manageable output can be generated by using the wrapper from the function MRT() of the MVPARTwrap package. Plus, this other function allows identification of discriminant species.

| Find discriminant and indicator species
# Find discriminant species with MRT results
doubs.mrt.wrap<-MRT(doubs.mrt,percent=10,species=colnames(spe.hel))
summary(doubs.mrt.wrap)

# Extract indval p-values
doubs.mrt.indval<-indval(spe.hel,doubs.mrt$where) doubs.mrt.indval$pval

# Extract indicator species of each node, with its indval
doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval<=0.05)]
doubs.mrt.indval$indcls[which(doubs.mrt.indval$pval<=0.05)]

The main discriminant species of the first split are TRU, VAI and ABL. TRU and VAI contribute highly to the left leaf, and ABL is the most indicative species of the sites at lower altitude (<361.5m). This output also indicates which sites are included in each leaf.

The second part of the code allows us to test the significance of the indicator value of each species through a permutation test. For each significant indicator species, we extracted the leaf number and the indicator value. In this particular case, TRU, VAI and LOC are all significant species of the left leaf, TRU having the highest indicator value (0.867).

Challenge 4: Run the multivariate regression tree for the mite data. Select the minimum size of tree within one SE of the CVRE. What is the proportion of variance explained by this tree? How many leaves contain this tree? What are the discriminant species?

Challenge 4 - Solution

Click to display ⇲

Click to hide ⇱

Click to hide ⇱

| Create a MRT with the mite data
mite.mrt<-mvpart(data.matrix(mite.spe.hel)~.,mite.env,
legend=FALSE,margin=0.01,cp=0,xv="pick",
xval=nrow(mite.spe.hel),xvmult=100,which=4)
summary(mite.mrt)

mite.mrt.wrap<-MRT(mite.mrt,percent=10,species=colnames(mite.spe.hel))
summary(mite.mrt.wrap)

mite.mrt.indval<-indval(mite.spe.hel,mite.mrt$where) mite.mrt.indval$pval

mite.mrt.indval$maxcls[which(mite.mrt.indval$pval<=0.05)]
mite.mrt.indval$indcls[which(mite.mrt.indval$pval<=0.05)]

25.6% of the variation in the mite species assemblage across sites is explained by the partition of the sites based on water content of the substrate (at 385.1 mg/l). LCIL is a discriminant species of sites with higher water content, and has an indicator value of 0.715.

# 4. Linear discriminant analysis

Linear discriminant analysis (LDA) is a constrained (canonical) technique that allows you to determine how well your independent set of variables explains an a priori grouping. This grouping may have been obtained from a previous clustering analysis (see Workshop 8) or from a hypothesis (e.g. grouping is based on sites at different latitudes or different treatments). An LDA can also be used to classify new data into these pre-determined groups. You can imagine some useful applications of this technique including assessing which population a fish should be classified in based on morphology or classifying whether a new paper is a freshwater, marine or terrestrial study based on the abstract of papers in those pre-determined biomes.

LDA computes discriminant functions from standardized descriptors. These coefficients quantify the relative contributions of the (standardized) explanatory variables to the discrimination of objects. Identification functions can be computed from the original (not standardized) descriptors to classify new data into pre-determined groups. Let’s continue to work with the Doubs fish data. First we must ensure that the within-group covariance matrices of the explanatory variables are homogeneous – a condition necessary for the application of LDA.

First we want to make an a priori classification that is independent from the environmental data set. We know that there is a general relationship that indicates environmental variables change with latitude (Budyko 1969). Here we will classify our Doubs fish sites based on latitude to determine how well the environmental factors explain our latitude grouping. Our groups are determined by simply dividing the range of latitudes equally into three groups and then assigning each site to a group depending on where they fall along the divided range.

| Load spatial data and classify Doubs fish sites based on latitude
#load spatial data to determine groups
spa <- spa[,-8]

#View spatial data
View (spa)

numbers<-(1:30)
numbers<-numbers[!numbers%in%8]
spa$site<-numbers #make groups based on lattitude y<82=group1, 82<y<156=group2, y>156=group3 spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3) #order by site spa.group<-spa.group[with(spa.group, order(site)), ] Generally, we would first want to check that the within-group covariance matrices of the explanatory variables are homogeneous by verifying multivariate homogeneity of within-group covariance (MHV). For the purposes of this workshop we will by pass it but more information can be found in Borcard et al. (2011). Once we run the LDA we can use the result object to determine 1. What groups the sites are classified in based on the environmental data. 2. What are the posterior probabilities of that the sites to belong to the groups. 3. The percentage of correct classification based on our latitudinal grouping. | Run the linear discriminant analysis (LDA) #run LDA LDA<-lda(env,spa.group[,4]) #classification of the objects based on LDA spe.class <- predict(LDA)$class

#posterior probabilities of the objects to belong to the groups
spe.post <- predict(LDA)$posterior #table of prior versus predicted classifications spe.table <- table(spa.group[,4], spe.class) #proportion of correct classification diag(prop.table(spe.table, 1)) The results suggest that the environmental factors explain the first, lower latitude, group and the group 3 perfectly but only 83% of the group 2 sites were predicted correctly. What does that tell us about our classification? Perhaps there are stronger delineations in the lower and higher latitude and the group 2 is a mix of both? Now what we have some new sites and we want to classify them based on the relationship we have established between our latitudinal grouping and environmental factors using the LDA. Using the predict() function we can load in a new matrix with sites and classify them using the LDA object. Load in the classifyme.csv file, which contains dummy data from 5 new sites. | Load the classify me data and predict grouping of new data #predicting classification of new data #read in new sites classify.me<-read.csv("classifyme.csv", header = T) #predict grouping of new data predict.group<-predict(LDA, newdata=classify.me) #give classification for each new site group.new<-predict.group$class

Our new sites, in order, have been classified in groups 1,1, 1, 3 and 3 respectively.

Challenge 5: Run an LDA for the mite env data (only first two vars) based on four latitudinal groups you create from the mite.xy data set. What group was group 2 most incorrectly grouped into? What proportion of sites was correctly classified in group 1? group 2?

Challenge 5: Solution

Click to display ⇲

Click to hide ⇱

Click to hide ⇱

| LDA on mite data
mite.xy$site<-seq(1:70) (max(mite.xy[,2])-min(mite.xy[,2]))/4 mite.xy.group<-ddply(.data=mite.xy, .variables=.(x, y, site), .fun= summarise, group = if(y <= 2.5) 1 else if (y <= 4.9) 2 else if (y <= 7.3) 3 else 4) mite.xy.group<-mite.xy.group[with(mite.xy.group, order(site)), ] LDA.mite<-lda(mite.env[,1:2],mite.xy.group[,4]) mite.class <- predict(LDA.mite)$class
mite.post <- predict(LDA.mite)\$posterior
mite.table <- table(mite.xy.group[,4], mite.class)
diag(prop.table(mite.table, 1))

# 5. Some other useful ordination methods

| other methods
?cca #(constrained correspondence analysis)
# Constrained Correspondence Analysis (CCA) is a canonical ordination method similar to RDA that preserve
# Chi-square distances among object (instead of Euclidean distances in RDA). This method is well suited for the
# analysis of large ecological gradients.

?CCorA # Canonical Correlation Analysis

# Canonical Correlation Analysis (CCorA) differs from RDA given that the two matrices are considered symmetric
# while in RDA the Y matrix is dependent on the X matrix. The main use of this technique is to test the
# significance of the correlation between two multidimensional data sets, then explore the structure of the data by
# computing the correlations (which are the square roots of the CCorA eigenvalues) that can be found between
# linear functions of two groups of descriptors.

#Coinertia Analysis (CoIA) is a symmetric canonical ordination method that is appropriate to compare pairs
# of data sets that play equivalent roles in the analysis. The method finds a common space onto which the objects
# and variables of these data sets can be projected and compared. Compared to CCorA, co-inertia analysis
# imposes no constraint regarding the number of variables in the two sets, so that it can be used to compare
# ecological communities even when they are species-rich. Co-inertia analysis is not well-suited, however, to
# analyse pairs of data sets that contain the same variables, because the analysis does not establish one-to-one
# correspondences between variables in the two data sets; the method does not ‘know’ that the first variable is the
# same in the first and the second data sets, and likewise for the other variables.

help(mfa, package=ade4) # Multiple Factorial Analysis

# Multiple factor analysis (MFA) can be used to compare several data sets describing the same objects. MFA
# consists in projecting objects and variables of two or more data sets on a global PCA, computed from all data
# sets, in which the sets receive equal weights.

# Spatial analysis can be performed using packages AEM and PCNM : http://r-forge.r-project.org/R/?group_id=195

# References

Alday & Marrs (2014). A simple test for alternative states in ecological restoration: the use of principal response curves. Journal of Vegetation Science, 17, 302-311.

Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.

Breiman, L., J. H. Friedman, et al. (1984). Classification and Regression Trees. Belmont, California, USA, Wadsworth International Group.

Budyko, M.I. (1969) The effect of solar radiation variations on the climate of the Earth. Tellus, 21(5), 611-619.

Clarke & Warwick (2001). Change in Marine Communities: An Approach to Statistical Analysis and Interpretation 2nd edition. Primer-E Ltd.

De'ath, G. (2002). Multivariate regression trees : a new technique for modeling species-environment relationships. Ecology, 83(4), 1105–1117.

Gotelli & Ellison (2004). A Primer of Ecological Statistics. Sinaeuer Associates Inc., Sunderland MA.

Legendre & Legendre (2012). Numerical Ecology 3rd edition. Elsevier Science BV, Amsterdam.

Poulin, Andersen & Rochefort (2013) A new approach for tracking vegetation change after restoration: a case study with peatlands. Restoration Ecology, 21, 363-371.


text:
# American Institute of Mathematical Sciences

January  2011, 15(1): 231-254. doi: 10.3934/dcdsb.2011.15.231

## Study on the stability and bifurcations of limit cycles in higher-dimensional nonlinear autonomous systems

 1 School of Mathematics and Computer Science, Fujian Normal University, Fuzhou, 350007, China 2 Department of Applied Mechanics and Engineering, Sun Yat-sen University, Guangzhou, 510275, China, China

Received  May 2009 Revised  March 2010 Published  October 2010

A semi-analytical procedure for studying stability and bifurcations of limit cycles in higher-dimensional nonlinear autonomous dynamical systems is developed. This procedure is based mainly on the incremental harmonic balance (IHB) method. It is composed of three key steps, namely, the determination of limit cycles by IHB method, the calculation of transition matrix by precise integration (PI) algorithm and the discrimination of limit cycle stability by Floquet theory. As an application, the procedure is used to investigate the dynamics of the limit cycle of a three-dimensional nonlinear autonomous system. The symmetry-breaking bifurcation, the first and the second period-doubling bifurcations of the limit cycle are identified. The critical parameter values corresponding to these bifurcations are calculated. The phase portraits and bifurcation points agree well with those of direct numerical integrations by using Runge-Kutta method.
Citation: Jianhe Shen, Shuhui Chen, Kechang Lin. Study on the stability and bifurcations of limit cycles in higher-dimensional nonlinear autonomous systems. Discrete & Continuous Dynamical Systems - B, 2011, 15 (1) : 231-254. doi: 10.3934/dcdsb.2011.15.231
##### References:

show all references

##### References:
 [1] Julien Barral, Yan-Hui Qu. On the higher-dimensional multifractal analysis. Discrete & Continuous Dynamical Systems - A, 2012, 32 (6) : 1977-1995. doi: 10.3934/dcds.2012.32.1977 [2] Aymen Balti, Valentina Lanza, Moulay Aziz-Alaoui. A multi-base harmonic balance method applied to Hodgkin-Huxley model. Mathematical Biosciences & Engineering, 2018, 15 (3) : 807-825. doi: 10.3934/mbe.2018036 [3] Alexandru Kristály. Asymptotically critical problems on higher-dimensional spheres. Discrete & Continuous Dynamical Systems - A, 2009, 23 (3) : 919-935. doi: 10.3934/dcds.2009.23.919 [4] N. Kamran, K. Tenenblat. Periodic systems for the higher-dimensional Laplace transformation. Discrete & Continuous Dynamical Systems - A, 1998, 4 (2) : 359-378. doi: 10.3934/dcds.1998.4.359 [5] Misha Bialy. Maximizing orbits for higher-dimensional convex billiards. Journal of Modern Dynamics, 2009, 3 (1) : 51-59. doi: 10.3934/jmd.2009.3.51 [6] Xinru Cao. Global bounded solutions of the higher-dimensional Keller-Segel system under smallness conditions in optimal spaces. Discrete & Continuous Dynamical Systems - A, 2015, 35 (5) : 1891-1904. doi: 10.3934/dcds.2015.35.1891 [7] Jiashan Zheng, Yifu Wang. A note on global existence to a higher-dimensional quasilinear chemotaxis system with consumption of chemoattractant. Discrete & Continuous Dynamical Systems - B, 2017, 22 (2) : 669-686. doi: 10.3934/dcdsb.2017032 [8] Randa Ben Mahmoud, Hichem Chtioui. Prescribing the scalar curvature problem on higher-dimensional manifolds. Discrete & Continuous Dynamical Systems - A, 2012, 32 (5) : 1857-1879. doi: 10.3934/dcds.2012.32.1857 [9] Youshan Tao. Global dynamics in a higher-dimensional repulsion chemotaxis model with nonlinear sensitivity. Discrete & Continuous Dynamical Systems - B, 2013, 18 (10) : 2705-2722. doi: 10.3934/dcdsb.2013.18.2705 [10] Andrey Gogolev. Smooth conjugacy of Anosov diffeomorphisms on higher-dimensional tori. Journal of Modern Dynamics, 2008, 2 (4) : 645-700. doi: 10.3934/jmd.2008.2.645 [11] Maoan Han, Tonghua Zhang. Some bifurcation methods of finding limit cycles. Mathematical Biosciences & Engineering, 2006, 3 (1) : 67-77. doi: 10.3934/mbe.2006.3.67 [12] Elena Bonetti, Pierluigi Colli, Gianni Gilardi. Singular limit of an integrodifferential system related to the entropy balance. Discrete & Continuous Dynamical Systems - B, 2014, 19 (7) : 1935-1953. doi: 10.3934/dcdsb.2014.19.1935 [13] Jueyou Li, Guoquan Li, Zhiyou Wu, Changzhi Wu, Xiangyu Wang, Jae-Myung Lee, Kwang-Hyo Jung. Incremental gradient-free method for nonsmooth distributed optimization. Journal of Industrial & Management Optimization, 2017, 13 (4) : 1841-1857. doi: 10.3934/jimo.2017021 [14] Inmaculada Baldomá, Ernest Fontich, Rafael de la Llave, Pau Martín. The parameterization method for one- dimensional invariant manifolds of higher dimensional parabolic fixed points. Discrete & Continuous Dynamical Systems - A, 2007, 17 (4) : 835-865. doi: 10.3934/dcds.2007.17.835 [15] Yunming Zhou, Desheng Shang, Tonghua Zhang. Seventeen limit cycles bifurcations of a fifth system. Conference Publications, 2007, 2007 (Special) : 1070-1081. doi: 10.3934/proc.2007.2007.1070 [16] Maoan Han, Yuhai Wu, Ping Bi. A new cubic system having eleven limit cycles. Discrete & Continuous Dynamical Systems - A, 2005, 12 (4) : 675-686. doi: 10.3934/dcds.2005.12.675 [17] Salomón Rebollo-Perdomo, Claudio Vidal. Bifurcation of limit cycles for a family of perturbed Kukles differential systems. Discrete & Continuous Dynamical Systems - A, 2018, 38 (8) : 4189-4202. doi: 10.3934/dcds.2018182 [18] Nobu Kishimoto. Resonant decomposition and the $I$-method for the two-dimensional Zakharov system. Discrete & Continuous Dynamical Systems - A, 2013, 33 (9) : 4095-4122. doi: 10.3934/dcds.2013.33.4095 [19] Jiashan Zheng. Boundedness of solutions to a quasilinear higher-dimensional chemotaxis-haptotaxis model with nonlinear diffusion. Discrete & Continuous Dynamical Systems - A, 2017, 37 (1) : 627-643. doi: 10.3934/dcds.2017026 [20] Ruijun Zhao, Yong-Tao Zhang, Shanqin Chen. Krylov implicit integration factor WENO method for SIR model with directed diffusion. Discrete & Continuous Dynamical Systems - B, 2019, 24 (9) : 4983-5001. doi: 10.3934/dcdsb.2019041

2018 Impact Factor: 1.008

format: plain


text:
# Deduction of centripetal force

1. Feb 9, 2010

### Misr

I have questions about the deduction of law of centripetal force

http://img42.imageshack.us/img42/5841/centripetalforce.gif [Broken]
Why do we assume that the change in direction of velocity of the object is V(t+dt)
And similarly the change in direction of radius?

Also what is the triangle of velocity ? and how to draw it?
http://img203.imageshack.us/img203/4669/triangleofvelocity2.gif [Broken]
Do we use only the magnitude of these directions to express the sides triangle of velocity bec. The triangle sides can't be negative while the directions could be +ve or negative?

Thanks

Last edited by a moderator: May 4, 2017
2. Feb 10, 2010

### Staff: Mentor

The velocity (for uniform circular motion) is always tangent to the circle. V(t+dt) is the velocity at time t+dt, not the change in velocity. The change in velocity is V(t+dt) - V(t).

You already have a drawing of the triangle of velocity. It represents the vector addition: V(t+dt) = V(t) + ΔV. The "trick" is realizing that the angle between the two velocity vectors is θ.

When drawing a vector, the length of the arrow represents the magnitude, which is always positive. The direction of the vector is represented by the direction of the arrow.

Last edited by a moderator: Apr 24, 2017
3. Feb 10, 2010

### Misr

Thanks this really helped me but actually i can't understand this point

You already have a drawing of the triangle of velocity. It represents the vector addition: V(t+dt) = V(t) + ΔV. The "trick" is realizing that the angle between the two velocity vectors is θ.

----
I have another question
can we prove that the triangle of velocity in the second circle is similar to the triangular cab (ABC) in the first circle?

Thanks again

4. Feb 11, 2010

### Staff: Mentor

Two points:
(1) The angle between the two velocity vectors is θ. To understand that, realize that the velocity vector is always perpendicular to the radius. Thus when the radius moves by an angle θ, so must the velocity vector.
(2) The triangles are similar. In the first circle you have an isosceles triangle with angle θ formed by the radii; In the second circle you have an isosceles triangle with angle θ formed by the velocity vectors. Two isosceles triangles with the same angle are similar.

Another point: For a small change in angle dθ, the length AB (in the first diagram) will equal rdθ.

5. Feb 11, 2010

### Misr

mmm actually i can't imagine can u draw an illustration?
i'll put a proof later to make sure

Thanks so much

Last edited: Feb 11, 2010
6. Feb 11, 2010

### Misr

bec ΔV = V(t+dt) - V(t)
therefore V(t+dt) = V(t) + ΔV ?Right?

7. Feb 12, 2010

### Staff: Mentor

The velocity vector is always a fixed angle (90°) with respect to the radius. If the radius points in the +x direction, the velocity vector points in the +y direction. When the radius has moved to be 30° (for example) above the x-axis, where must the velocity vector point?
Right.

8. Feb 12, 2010

### Misr

Thanks so much
I want to ask if this works ?

http://img215.imageshack.us/img215/3401/proof3f.jpg [Broken]
http://img693.imageshack.us/img693/8039/proofe.jpg [Broken]
http://img705.imageshack.us/img705/5425/proof2t.jpg [Broken]

I'm very sorry 4 the handwriting but i'm in a hurry

Thanks again

Last edited by a moderator: May 4, 2017
9. Feb 12, 2010

### Misr

Also i made an animated picture to demonstrate ur words
http://img30.imageshack.us/img30/1725/centripetalanimated.gif [Broken]
Is my understanding right?

I want to make sure if my two posts are right

Last edited by a moderator: May 4, 2017
10. Feb 12, 2010

### Staff: Mentor

Your animation looks fine to me.

11. Feb 13, 2010

### Misr

Good
I think i got it now
but i have another way and also i want to make sure if its right ?

[url]http://img24.imageshack.us/i/prooffz.jpg/[/url]

IS this right

Thanks so much

12. Feb 14, 2010

### Staff: Mentor

I'm afraid that I don't understand what you are doing in that image. Could you define your triangles and angles, and then describe in words what you are trying to prove?

13. Feb 14, 2010

### Misr

I'm very sorry its my fault I should have explained more
Actually I don't really need to put this proof but I'm just eager to know if its right or wrong
http://img246.imageshack.us/img246/745/centripetalforce2.jpg [Broken]
This pink line is an extension to the velocity vector
And I drew a green radius to help me in the proof
I don't want to prove any thing new I just want to prove that the angle between the velocity vectors is equal to the angle between the two radii .
I understand what u said but I want to make sure if this method also works .

Since the velocity vector is at point b is parallel to ax
Therefore m(<MXY)=m(MXA)=m(<b) = 90 degrees

Therefore XA = XY (According to a theorem in geometry )

And MX is a common side in the two trianles MXY and MXA
Therefore triangle MXY is congruent to triangle MXA

Therefore m(<2)=m(<3)

Since m(<1) = 1/2 m(<AMY) (angle of tangency and centeral angle subtended by the same arc)

Therefore m(<1) = m(<2)
I hope its obvious now
Thanks

Last edited by a moderator: May 4, 2017
14. Feb 14, 2010

### Staff: Mentor

Looks good to me. You can also say that m(<1) = m(<2) since they are both complementary to angle MAX.

15. Feb 15, 2010

### Misr

Yeah this way is better and much more easier
this is a very interesting discussion
nothing is well-illustrated in my txt book

thanks so so much this really helped

16. Feb 15, 2010

### ideasrule

As an aside, there's a much easier derivation of this.

$$\vec v \cdot \vec r =C$$

Derive both sides:

$$\vec a \cdot \vec r + \vec v \cdot \vec v =0$$

Get rid of the vectors:

$$-ar + v^2 = 0$$

17. Feb 15, 2010

### Staff: Mentor

Nice. Combine that with a similar treatment of $\vec v \times \vec r = C$ to show that the acceleration is radial.


text:
# Nano's question at Yahoo! Answers regarding Newton's Law of Cooling

#### MarkFL

Staff member
Here is the question:

I need help with this problem !?

Sandy manages a ceramics shop and uses a 800°F kiln to fire ceramic greenware. After turning off her kiln, she must wait until its temperature gauge reaches 155°F before opening it and removing the ceramic pieces. If room temperature is 70°F and the gauge reads 700°F in 10 minutes, how long must she wait before opening the kiln? Assume the kiln cools according to Newton's Law of Cooling.
(Round you anwer to the nearest whole minute)
a) 512 mins B)94 mins c)369 mins D)146 mins
I have posted a link there to this topic so the OP can see my work.

#### MarkFL

Staff member
Hello Nano,

Newton's law of Cooling states that the time rate of change of the temperature $T$ of an object is proportional to the difference between the ambient temperature $M$ and the temperature of the object. Stated mathematically, this is:

$$\displaystyle \frac{dT}{dt}=-k(T-M)$$ where $$\displaystyle T(0)=T_0,\,0<k\in\mathbb{R}$$ and $$\displaystyle T>M$$.

The ODE is separable and may be written:

$$\displaystyle \frac{1}{T-M}\,dT=-k\,dt$$

Integrating, using the boundaries, and dummy variables of integration, we find:

$$\displaystyle \int_{T_0}^{T(t)}\frac{1}{u-M}\,du=-k\int_0^t v\,dv$$

$$\displaystyle \ln\left(\frac{T(t)-M}{T_0-M} \right)=-kt$$

If we know another point $$\displaystyle \left(t_1,T_1 \right)$$ we may now find $k$:

$$\displaystyle -k=\frac{1}{t_1}\ln\left(\frac{T_1-M}{T_0-M} \right)$$

Hence, we find:

$$\displaystyle t=\frac{t_1\ln\left(\frac{T(t)-M}{T_0-M} \right)}{\ln\left(\frac{T_1-M}{T_0-M} \right)}$$

Using the following given data:

$$\displaystyle t_1=10,\,T(t)=155,\,M=70,\,T_0=800,\,T_1=700$$

we find:

$$\displaystyle t=\frac{10\ln\left(\frac{155-70}{800-70} \right)}{\ln\left(\frac{700-70}{800-70} \right)}=\frac{10\ln\left(\frac{17}{146} \right)}{\ln\left(\frac{63}{73} \right)}\approx145.9628332694268$$

Therefore d) is the correct answer.


text:
37 views
sum of how many term of series

6+12+18+24+… is 1800?
| 37 views

sun of n terms of AP = N/2{2a+(n-1)*d} where a= first term .d=common difference ,n= number of terms

now put values in the formula a=6, d=6 , n=?,sum=1800

N/2{2*6+(n-1)*6}=1800

$n^{2}+n-600=0$

solving the quadratic equation you will get 24 terms so answer is 24
by (42 points) 1


text:
# Definition:Total Ordering

(Redirected from Definition:Linear Ordering)

## Definition

Let $\mathcal R \subseteq S \times S$ be a relation on a set $S$.

### Definition 1

$\RR$ is a total ordering on $S$ if and only if:

$(1): \quad \RR$ is an ordering on $S$
$(2): \quad \RR$ is connected

That is, $\RR$ is an ordering with no non-comparable pairs:

$\forall x, y \in S: x \mathop \RR y \lor y \mathop \RR x$

### Definition 2

$\RR$ is a total ordering on $S$ if and only if:

$(1): \quad \RR \circ \RR = \RR$
$(2): \quad \RR \cap \RR^{-1} = \Delta_S$
$(3): \quad \RR \cup \RR^{-1} = S \times S$

## Also known as

Some sources call this a linear ordering, or a simple ordering.

If it is necessary to emphasise that a total ordering $\preceq$ is not strict, then the term weak total ordering may be used.

## Also see

• Results about total orderings can be found here.


text:
## Double Doors Regional P5 - Tudor Eats Ice Cream Cones

View as PDF

Points: 15
Time limit: 1.4s
Java 2.75s
Memory limit: 256M

Problem type
Allowed languages
Ada, Assembly, Awk, Brain****, C, C#, C++, COBOL, CommonLisp, D, Dart, F#, Forth, Fortran, Go, Groovy, Haskell, Intercal, Java, JS, Kotlin, Lisp, Lua, Nim, ObjC, OCaml, Octave, Pascal, Perl, PHP, Pike, Prolog, Python, Racket, Ruby, Rust, Scala, Scheme, Sed, Swift, TCL, Text, Turing, VB, Zig

Tudor likes ice cream cones!

One day, Tudor travels to an ice cream convention to eat ice cream. There are stands, numbered through . Stand is serving flavor .

Each of the stands charges a certain price for ice cream, which Tudor interprets as a signal of the quality of the ice cream. Tudor has strong flavor preferences and will decide to order ice cream from two stands within certain flavor constraints, choosing to either order from the two stands with the highest price, in the event he wishes to try some gourmet ice cream, or the two stands with the lowest price, in the event he wishes to save on money.

As the convention goes on, stands may change their prices. Tudor's flavor preferences may also change. Help Tudor eat ice cream!

#### Input Specification

The first line contains two positive integers, and . You may assume and are no larger than . Furthermore, will be a perfect square and will not be equal to 1.

The next line contains positive integers, the costs of each of the stands in order. No stand will exceed in cost.

Each of the next lines will contain three integers. The first integer in each line will be either , , or . The remaining integers on the line will be encoded. All encoded integers should be decoded by XOR'ing them with the answer of the most recent query. In the event there were no prior answers, XOR the encoded integers with .

If the first integer in the line is 1, the next two integers are and , indicating that stand should have value . The decoded value will satisfy . The decoded value will satisfy .

If the first integer in the line is 2, the next two integers are and , indicating that Tudor wishes to know the sum of the prices of the two cheapest stands that offer a flavor between and . The decoded values and will satisfy .

If the first integer in the line is 3, the next two integers are and , indicating that Tudor wishes to know the sum of the prices of the two most expensive stands that offer a flavor between and . The decoded values and will satisfy .

#### Output Specification

For each query, output the desired cost on a separate line.

#### Sample Input

4 4
1 4 2 5
2 1 1
2 7 7
1 3 6
3 6 5

#### Sample Output

5
7
6


text:
# Flushing all current figures in matplotlib

To flush all current figures in matplotlib, use close('all') method.

## Steps

• Set the figure size and adjust the padding between and around the subplots.

• Create a figure with the title "First Figure".

• Create another figure with the title "Second Figure".

• To close all figures, use close('all').

• To display the figure, use show() method.

## Example

from matplotlib import pyplot as plt

plt.rcParams["figure.figsize"] = [7.50, 3.50]
plt.rcParams["figure.autolayout"] = True

plt.figure("First Figure")
plt.figure("Second Figure")

# plt.close('all')

plt.show()

## Output

Notice that we have commented the line −

plt.close('all')

Hence, it will display two figures −

Uncomment the line plt.close('all') and run the code again. It will flush all the current figures.


text:
$$\renewcommand\AA{\unicode{x212B}}$$

# Mantid 3.10.1 Release Notes¶

This is a patch release that corrects some significant issues since version 3.10.0.

There is no common theme to the fixes contained in this patch, rather it is collection of small but significant fixes. Please see below for the full list of changes.

## Citation¶

Please cite any usage of Mantid as follows:

• Mantid 3.10.1: Manipulation and Analysis Toolkit for Instrument Data.; Mantid Project. doi: 10.5286/Software/Mantid3.10.1

• Arnold, O. et al. Mantid-Data Analysis and Visualization Package for Neutron Scattering and mu-SR Experiments. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 764 (2014): 156-166 doi: 10.1016/j.nima.2014.07.029 (download bibtex)

## Changes in this version¶

• 19988 updated the amplitude in GausOsc

• 19995 Add support for 2d character arrays for sample name

• 20001 Fix type conversion warning on MSVC

• 20011 Copy units over logs

## Summary of impact¶

Issue

Impact

Solution

Side Effect Probability

19988

The initial amplitude in the GausOsc functionis unrealistic

Change initial amplitude to a more realistic value is 0.2

low

19995

SNS DAS now writes the sample name as a 2d character

Add support for 2d character arrays for sample name

low

20001

Type conversion warning on MSVC

Fix warning

low

20011

The logs do not get their units copied over.

Copy units over logs

low


text:
nLab Tate spectrum

Contents

Context

Representation theory

representation theory

geometric representation theory

Contents

Idea

For ∞-actions of finite groups $G$ on objects $E$ in stable (∞,1)-categories, then the homotopy cofiber $X^{t G}$ of the norm map is called the Tate construction, sitting in a homotopy fiber sequence

$X_G \stackrel{norm}{\longrightarrow} X^G \longrightarrow X^{t G} \,.$

(e.g Lurie, def. 6.1.6.24)

Depending on which edition you have, chapter 6 may be chapter 7.

For the stable (∞,1)-category of spectra this is accordingly called the Tate spectrum.

References

The general abstract discussion is due to

Review with an eye towards topological cyclic homology is in

Last revised on July 25, 2017 at 11:57:16. See the history of this page for a list of all contributions to it.


text:
# Reinventing the basic calculator functions

I decided to reinvent the basic calculator functions using only:

• the set of positive integers

• greater than, less than, and equality properties

• increment and decrement

(in other words, the counting on fingers method)

def is_positive_int(x):
if type(x) == int and x >= 0:
return True
else:
return False

if is_positive_int(x) and is_positive_int(y):
for i in range(y):
x += 1
return x
else:
return False

def multiply(x,y):
if is_positive_int(y) and is_positive_int(y):
m = 0
for i in range(y):
return m
else:
return False

def power(x,y):
""" x to the y power """
if is_positive_int(x) and is_positive_int(y):
p = 1
for i in range(y):
p = multiply(p,x)
return p
else:
return False

def subtract(x,y):
""" x - y """
if is_positive_int(x) and is_positive_int(y) and x > y:
for i in range(y):
x -= 1
return x
else:
return False

def divide(x,y):
""" x / y """
if is_positive_int(x) and is_positive_int(y):
floor = 0
remainder = 0
for i in range(x):
remainder += 1
if remainder == y:
floor += 1
remainder = 0
return {'floor':floor, 'modulus':remainder}
else:
return False

def floor(x,y):
""" a // b"""
return divide(x,y)['floor']

def modulus(x,y):
""" x % y"""
return divide(x,y)['modulus']

• Instead of if type(x) == int you could've use isinstance(x, int). It's not that isinstance is good, mind you — it's just less bad than checking equality of types. – Grajdeanu Alex. Aug 18 '16 at 10:13

Let's say you have this situation:

if ...:
...
return ...
else:
return ...

What if the if condition is True? Some statements will be executed, and then the function will return. In other words, if the condition is True, execution never leaves the if block. Therefore, you really don't need the else. Just put the return directly under the if:

if ...:
...
return ...
return ...

It is generally a bad idea to use the pattern type(...) == .... Usually, we use isinstance(..., ...). That way, subclasses are allowed. Of course, you might keep it if you don't want subclasses to be allowed, but I don't see a reason to restrict the user. I'm assuming you want only integers because only integers can be passed to range, but subclasses of integers can indeed also be passed. For example, booleans are integers. range(False, True) is perfectly valid. Also, True + False + False is valid, so it won't mess up your arithmetic. That makes your is_positive_int function look like this:

def is_positive_int(x):
if isinstance(x, int) and x >= 0:
return True
return False

But what do you see? If a condition is True, we return True. If it is False, we return False. In other words, we are really just returning the condition, so why not do it like this:

def is_positive_int(x):
return isinstance(x, int) and x >= 0

I see a very similar pattern in all of your functions:

if is_positive_int(x) and is_positive_int(y):
...
else:
return False

First of all, we can remember what I mentioned earlier about not leaving the if block unless the condition is False. Therefore, I would actually switch the if and else:

if not is_positive_int(x) or not is_positive_int(y):
return False
...

We remove a level of nesting that way. The problem is that we use the same code in a bunch of different functions, so we can make a decorator. A decorator is a function which, when given another function (a) to wrap, returns a wrapper function (b) for a. b overrides a and does its special thing, but it usually calls a later. This is my idea:

def guarantee_positive(function):
# Each function takes two arguments
def wrapper(x, y):
if is_positive_int(x) and is_positive_int(y):
# Deal with it normally
return function(x, y)
return False
# Now wrapper(x, y) will be called instead of function(x, y)
return wrapper

To use it, add @guarantee_positive before each function definition. For example:

@guarantee_positive
def divide(x, y):
floor = 0
remainder = 0
for i in range(x):
remainder += 1
if remainder == y:
floor += 1
remainder = 0
return {'floor': floor, 'modulus': remainder}

That makes for a much cleaner function. Of course, subtract has an extra condition, but I'll get to that later.

With all that about function decorators, why do you need them both to be positive? I can see why you might want y to be positive: because range doesn't work the same with negative numbers, but you aren't using i most of the time, so just use the absolute value of y. I'll give an example:

def multiply(x, y):
m = 0
function = add if y > 0 else subtract
for _ in range(abs(y)):
m = function(m, x)
return m

It doesn't require a decorator, and it is capable of more. You can now use it with any combination of negative and positive:

-3,  3 -> -9
-3, -3 ->  9
3, -3 -> -9
3,  3 ->  9

Your add function would be like this:

def add(x, y):
# Add 1 if y is positive; subtract 1 if it is negative
y_add = (1 if y >= 0 else -1)
for _ in range(abs(y)):
return x

Again, it's shorter than the original, and it does more.

There are some functions where it isn't quite that easy. For example, power() would be difficult to implement a negative y (x doesn't need to be positive even with the current implementation), but most of the others don't need the decorator. I won't tell you all of them, but I have given examples of the more common hangups.

Back to the decorator, I think it is better to raise an exception when something is wrong than to return False. After all, they are called exceptions for a reason. Normally, the function does one thing, but there are exceptions (or there should be).

I see a lot of for i in range(y) when i isn't used. The common way of signifying that i is just a place holder is using _ as the variable name.

It looks like you cheated with divide() :) You don't actually have division, you have divmod. I think a name similar to that would be more appropriate. another thing is that you return a dictionary. I prefer something closer to how the built-in function does it:

return floor, remainder

That way, you can do:

floor, remainder = divide(...)

Instead of the equivalent:

result = divide(...)
floor, remainder = result['float'], result['modulus']

Either one should have a doc string that explains what the return value is. If I called divide(...) when it had a doc string of just x / y, I would assume that I would get x / y, not a dictionary with the floored result and the remainder. Your floor and modulus functions would then be:

def floor(x, y):
"""a // b"""
return divide(x, y)[0]

def modulus(x, y):
"""x % y"""
return divide(x, y)[1]

Except that divide will hopefully have a more descriptive name.

I did notice that floor is inside of divide(), but I am assuming that is just a typo.

You seem to be following PEP 8 fairly well, but I noticed that you don't have any whitespace between the key-value pairs in your dictionary. PEP 8 recommends that you put a space between the colon and the value. I don't see it in the PEP, but I think you should put a space between a comma and the following object. For example:

def divide(x, y):

not

def divide(x,y):

Without the space, it feels smooshed, like you're golfing (code golfing, not with a golf club).

Why not go a step further and implement this as a class, this way you learn something about magic methods as well :)

class Fingers(int):
"""A class implementing counting with fingers instead of the more effective built-in int methods"""

def __init__(self, x):
if not isinstance(x, int):
raise TypeError("Input must be 'int' (or subclass)")
self.x = x

def __pos__(self):
return self

def __neg__(self):
return Fingers(-self.x)

def __abs__(self):
return Fingers(abs(self.x))

"""a + b, a and b may be negative"""
y_add = (1 if other >= 0 else -1)
x = self.x
for _ in range(abs(other)):
return Fingers(x)

def __sub__(self, other):
return self + (-other)

def __mul__(self, other):
"""a * b, a and b may be negative"""
m = Fingers(0)
for _ in range(abs(other)):
m = m + self
return m

def __pow__(self, exponent):
"""a**b, b must be positive"""
if exponent < 0:
raise NotImplementedError
p = Fingers(1)
for _ in range(y):
p = p * self
return p

def __divmod__(self, other):
floor = Fingers(0)
remainder = Fingers(0)
for _ in range(self):
remainder = remainder + Fingers(1)
if remainder == other:
floor = floor + Fingers(1)
remainder = Fingers(0)
return floor, remainder

def __floordiv__(self, other):
"""a//b"""
return self.__divmod__(other)[0]

def __mod__(self, other):
"""a%b"""
return self.__divmod__(other)[1]

if __name__ == "__main__":
x = Fingers(3)
y = Fingers(2)
assert x + y == Fingers(5)
assert x * y == Fingers(6)
assert x**y == Fingers(9)
assert divmod(x, y) == (Fingers(1), Fingers(1))
assert x // y == Fingers(1)
assert x % y == Fingers(1)

(I did not implement __iadd__ etc, so x += 3 will not work. They are trivial to implement, though.)

• There was a decorator or maybe a base class that allowed you to implement only some operators; reflected and in-place would be filled in automagically, but I can't find it now. – Daerdemandt Aug 18 '16 at 14:04
• Found it! functools.total_ordering is the name. – Daerdemandt Oct 10 '16 at 16:47

It has already been said that you should've probably implemented this as a class. There's even something in the docs about that.

However, there's also one concept you could've made use of: Church numerals.

You could make your numbers have something like .to_church() method, so that CustomNumeral(3).to_church(f) would return something similar to lambda x: f(f(f(x))), it allows for very expressive illustration for how addition naturally arises from incrementing and multiplication - from addition.

• Hozier.to_church() print Hozier.sins sharpenKnife() return deathlessDeath – DasBeasto Aug 18 '16 at 19:26


text:
precursor-0.1.0.0: Prelude replacement

Precursor.Control.Alternative

Synopsis

# Alternative class

class Applicative f => Alternative f where #

A monoid on applicative functors.

If defined, some and many should be the least solutions of the equations:

• some v = (:) <\$> v <*> many v
• many v = some v <|> pure []

Minimal complete definition

Methods

empty :: f a #

The identity of <|>

(<|>) :: f a -> f a -> f a infixl 3 #

An associative binary operation

some :: f a -> f [a] #

One or more.

many :: f a -> f [a] #

Zero or more.

Instances

(<|>) :: Alternative f => forall a. f a -> f a -> f a #

An associative binary operation

empty :: Alternative f => forall a. f a #

The identity of <|>

some :: Alternative f => forall a. f a -> f [a] #

One or more.

many :: Alternative f => forall a. f a -> f [a] #

Zero or more.

# Utility functions

optional :: Alternative f => f a -> f (Maybe a) #

One or none.

afilter :: (Monad m, Alternative m) => (a -> Bool) -> m a -> m a Source #

A generalized version of filter, which works on anything which is both a Monad and Alternative.

\(Blind p) xs -> filter p xs === afilter p xs

asum :: (Foldable t, Alternative f) => t (f a) -> f a #

The sum of a collection of actions, generalizing concat.

guard :: Alternative f => Bool -> f () #

guard b is pure () if b is True, and empty if b is False.

ensure :: Alternative f => (a -> Bool) -> a -> f a Source #

ensure allows you to attach a condition to something

eitherA :: Alternative f => f a -> f b -> f (Either a b) Source #

eitherA is especially useful for parsers.

toAlt :: (Alternative f, Foldable t) => t a -> f a Source #

Convert any Foldable to an Alternative

mapAlt :: (Monad m, Alternative m, Foldable f) => (a -> f b) -> m a -> m b Source #

Map a function over a monad, and concat the results. This is a generalized form of the function mapMaybe.


text:
Welcome :Guest
Congratulations!!!

Top 5 Contributors of the Month
Easy Web
Imran Ghani

### is not declared. It may be inaccessible due to its protection level.

Hi

I try to create a new webb site using VS 2010 and .Net Framework 4.0, I wanted it to be a copy of an exsisting .Net Framework 3.5 based site, so I copied the aspx page code and the codebehind code from the old version when I created the new one, but now I get these errors...

is not declared. It may be inaccessible due to its protection level.

It seem I get this error on every control I have on the page, why is that and how do I fix it?

### '' is not declared. It may be inaccessible due to its protection level

Hello

I hope I am posting to the correct forum. I see lots of reports of similar errors, but no fix for me. I have a new VS2010 website project. I am unable to call public routines in a public class in AppCode.

To debug, I built a new blank website project with a public class Class1 in AppCode and a public function Foo returning a string "Hello World".

A sub CallFoo calls Foo without error if CallFoo is in Class1. If  CallFoo is moved to a new class Class2 I get the above error.

I have rebuilt this test site from scratch targeting Framework 3.5 and 4 with the same error. Oddly, I have found I am able to reference public constants without error.

Any suggestions appreciated.

Jen

### '__o' is not declared. It may be inaccessible due to its protection level.

Hi

I have 10-15 errors of this kind

'__o' is not declared. It may be inaccessible due to its protection level.

at my <%: start tags

Im new to this, so why do i get this error and how do i fix it.
I have looked at the numbers 0 that its not the letter oO, and the other way bag
And its ok...

Code:

Netto:

<%: (Model.Cart.ComputeTotalValue() * 0.8D).ToString("c")%>

Moms 25%:

<%: (Model.Cart.ComputeTotalValue() * 0.2D).ToString("c")%>

### 'Microsoft.SharePoint.MailMessage' is inaccessible due to its protection level

down vote favorite 'Microsoft.SharePoint.MailMessage' is inaccessible due to its protection level On this code: MailMessage mail = new MailMessage(); mail.From = "myemail@xxxxx.edu"; mail.To = "myemail@xxxx.edu"; mail.Subject = "Testing Code"; mail.BodyText = what; mail.Priority = MailPriority.High; Smtp.Send(mail, "smtp.xxxxx.edu"); How can remidy this? changes to web.config? Any way to circumvent in code? Thanks.

I am sure there is an easy answer to this but i am very new,

Child Page: code behind file

protected void Page_Load(object sender, EventArgs e)
{

}
protected void Button1_Click(object sender, EventArgs e)
{
SiteMaster master = (SiteMaster)this.Master;
}

Master: Code Behind

protected void Page_Load(object sender, EventArgs e)
{

}

public string NavigateUrl
{

set
{
}
get
{

### SSIS Package Protection Level Errors - decrypt protected XML node "DTS:Password" with error 0x800900

ErrorÃÂ 7ÃÂ Error loading GetData.dtsx: Failed to decrypt protected XML node "DTS:Password" with error 0x8009000B "Key not valid for use in specified state.". You may not be authorized to access this information. This error occurs when there is a cryptographic error. Verify that the correct key is available.ÃÂ ÃÂ ÃÂ c:\Packages\GetData.dtsxÃÂ 1ÃÂ 1ÃÂ  ErrorÃÂ 8ÃÂ Validation error. Data Flow Task:ÃÂ OLEDB SourceÃÂ [5310]: The AcquireConnection method call to the connection manager "SourceServer" failed with error code 0xC0202009.ÃÂ GetData.dtsxÃÂ 0ÃÂ 0ÃÂ Environment I need to execute - BIDSHave Solution, with files linked with VSS. When I click on package, aforesaid errors show up and package opens in designed.Protection Level : Encrypt Sesitive with User KeyThis means it'd run only on machine and by user where and who created respectively.So I changed it to Encrypt Sesitive with Passoword and provided Passoword. No luckSo I changed it to Do not save sensitive, as anyways in the connection manager, it has Expression saying @Variablename (with appropriate correct syntax), and variable is mapped to config file. I tired with Config file and also by hardcoding the value in the connection manager's connection string after removing the variable name, also tried by hardcoding value in variable value and disabling config file.Nothing works. I read 3 threads discussing this issue by othe

### how to store the values in a object variable (package level scope) in a script component

I have declared a variable of object type with package level scope.  I want to assign values to this variable inside a script component.  My object should hold 5 columns in a row.  How to implement this ?

### Cannot create MDX calculation with one variable in the calculation from a leaf level

I am new to MDX/SSAS and have a problem with setting up two calculations. I have spent several days researching, endless rewriting of MDX queries, and did not find an answer (or didn’t recognize the post as an answer to my problem :)

I need help with two calculations.

One, a calculation, let’s call it  [Dollar Variance], that computes variance at all levels by this formula:

([Actual Work Billings] - [Budgeted Work Billings]) / [Work Days in Month] = [Dollar Variance]

The problem is that while the [Actual Work Billings] and the [Budgeted Work Billings] are summed at the particular dimension intersects chosen, for the variance to be calculated properly, the

### How to get the Package level Variable collections (Global Variable) in Script Task

Hi All,

I have the requirement like read all the Variables in Script Component Task, which includes variables available at package level and their values, global variables declared at package.

Currently are two ways to read the variables but which are not suitable for my requirement

1) Mentioning the variables at Script designer window as ReadOnlyVariables or ReadWriteVariables

<

### Msg 137, Level 15, State 2, Procedure AddCMEXperience, Line 49 Must declare the scalar variable "@Re

Msg 102, Level 15, State 1, Procedure AddCMEXperience, Line 15
Incorrect syntax near ','.

Msg 137, Level 15, State 2, Procedure AddCMEXperience, Line 49
Must declare the scalar variable "@ReasonLeft".

Here is the Stored Procedure that is in question here

USE [ResponsedDatabase]
GO
/****** Object:  StoredProcedure [dbo].[AddNewEmployer]    Script Date: 10/03/2010 14:24:51 ******/
SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO

(
@MemberID bigint,
@EmployerID bigint,
@DateStarted smalldatetime,
@IsContinuing bit,
@DateLeft smalldatetime,
@PositionStarted nvarchar(50),
@PositionLeft nvarchar(50),
@StartingSalary nvarchar(20),
@EndingSalary nvarchar(20),
@TimeWorked float,
@JobResponsibilities text,
@ReasonLeft, nvarchar(200)
)

AS

INSERT CMExperience
(
MemberID,
EmployerID,
&nbs

### Cant see package level variable in script component

I am currently having a issue accessing a package level variable.

My package level variable is called RecordSet1 and type object.

I was able to use a Data Flow and populate the variable with data with a record set destination.

So, now I am writing a script component that makes a recordset as a source(I am new to ssis and I know this has performance issues, I am only doing this to get familiar with scripting cause my manager told me too) So I am not able to access any of the variables. My code is below. I am getting a error on this.Variables.RecordSet1. I looked all over online and unable to really find a solution. I know in script component, you can use the this.Variables.VariableName to access the variables but even my intellisense isn't showing anything.

Error 1 'Variables' does not contain a definition for 'RecordSet1' and no extension method 'RecordSet1' accepting a first argument of type 'Variables' could be found (are you missing a using directive or an assembly reference?) C:\Users\Dhinesh\AppData\Local\Temp\SSIS\ee6fc0aaa4f447e0b024451e97f51dec\main.cs 42 39 sc_95443d5e406b44d

### A Low Level Look at ASP.NET Architecture

ASP.NET is a powerful platform for building Web applications, that provides a tremendous amount of flexibility and power for building just about any kind of Web application. Most people are familiar only with the high level frameworks like WebForms and WebServices which sit at the very top level of the ASP.NET hierarchy. In this article I'll describe the lower level aspects of ASP.NET and explain how requests move from Web Server to the ASP.NET runtime and then through the ASP.NET HTTP pipeline to process requests.

### Local variable scoping in C#

In my previous post, Compiler-generated scopes for local variable declarations, I briefly touched on the issue of multiple meanings applied to the same name. In this post, I'll aim to flush out the compiler's rules with regards to binding names in their local scopes.

### Compiler-generated scopes for local variable declarations

I was tasked with understanding and fixing a bug on error reporting with foreach iteration variables the other day, and it got me thinking about local variable scoping rules in C# in general. First, the bug.

### Retrieving a variable from a previous page

Hi,

I'm looking to have a search page where a user will select an item from a listview, then then taken to another page that wil display data based on their original selection.

A user will type in a string and be presented with a list of usernames that are like that string, when they select a user name who's profile they would like to view I would like the page that they're taken to, to be that users profile.

How would I transfer the variable from the search page to the profile page, and when I work on the profile page, how can I access that variable?

Thanks.

### Embeb a variable number of videos

I need to provide some alternate content for my non-Flash users.

I have a directory of video files that is scanned using the DirectoryInfo() function, and I want to use the results to embed all files ending in ".mov" in QuickTime players on a web form.

How would you approach this? Right now I am using the count on the array of videos to determine how many times to repeat the embed code.  But this seems inelegant.

Is there a way to have a single QuickTime player.ascx control that gets duplicated according to the number of videos available for display?

### How to take value of bound field of grid in a variable or in a textbox.

How to take value of bound field of grid in a variable or in a textbox.

I m using this code in .aspx page for grid.

<asp:GridView ID="Grd_Mprocess" runat="server" AutoGenerateColumns="False"

Width="100%" AllowPaging="True"

onpageindexchanging="Grd_Mprocess_PageIndexChanging"

AllowSorting="True" PageSize="5"  >

<RowStyle CssClass="odd" />

<Columns>


text:
## 3.3 Lorenz Curve (svylorenz)

Though not an inequality measure in itself, the Lorenz curve is a classic instrument of distribution analysis. Basically, it is a function that associates a cumulative share of the population to the share of the total income it owns. In mathematical terms,

$L(p) = \frac{\int_{-\infty}^{Q_p}yf(y)dy}{\int_{-\infty}^{+\infty}yf(y)dy}$

where $$Q_p$$ is the quantile $$p$$ of the population.

The two extreme distributive cases are

• Perfect equality:
• Every individual has the same income;
• Every share of the population has the same share of the income;
• Therefore, the reference curve is $L(p) = p \text{ } \forall p \in [0,1] \text{.}$
• Perfect inequality:
• One individual concentrates all of society’s income, while the other individuals have zero income;
• Therefore, the reference curve is

$L(p)= \begin{cases} 0, &\forall p < 1 \\ 1, &\text{if } p = 1 \text{.} \end{cases}$

In order to evaluate the degree of inequality in a society, the analyst looks at the distance between the real curve and those two reference curves.

The estimator of this function was derived by Kovacevic and Binder (1997)Kovacevic, Milorad, and David Binder. 1997. “Variance Estimation for Measures of Income Inequality and Polarization - the Estimating Equations Approach.” Journal of Official Statistics 13 (1): 41–58. http://www.jos.nu/Articles/abstract.asp?article=13141.:

$L(p) = \frac{ \sum_{i \in S} w_i \cdot y_i \cdot \delta \{ y_i \le \widehat{Q}_p \}}{\widehat{Y}}, \text{ } 0 \le p \le 1.$

Yet, this formula is used to calculate specific points of the curve and their respective SEs. The formula to plot an approximation of the continuous empirical curve comes from Lerman and Yitzhaki (1989)Lerman, Robert, and Shlomo Yitzhaki. 1989. “Improving the Accuracy of Estimates of Gini Coefficients.” Journal of Econometrics 42 (1): 43–47. http://EconPapers.repec.org/RePEc:eee:econom:v:42:y:1989:i:1:p:43-47..

A replication example

In October 2016, (Jann 2016Jann, Ben. 2016. Estimating Lorenz and concentration curves in Stata.” University of Bern Social Sciences Working Papers 15. University of Bern, Department of Social Sciences.) released a pre-publication working paper to estimate lorenz and concentration curves using stata. The example below reproduces the statistics presented in his section 4.1.

# load the convey package
library(convey)

library(survey)

# load the stata-style webuse library
library(webuse)

# load the NLSW 1988 data
webuse("nlsw88")

# coerce that tbl_df to a standard R data.frame
nlsw88 <- data.frame( nlsw88 )

# initiate a linearized survey design object
des_nlsw88 <- svydesign( ids = ~1 , data = nlsw88 )
## Warning in svydesign.default(ids = ~1, data = nlsw88): No weights or
## probabilities supplied, assuming equal probability
# immediately run the convey_prep function on the survey design
des_nlsw88 <- convey_prep(des_nlsw88)

# estimates lorenz curve
result.lin <- svylorenz( ~wage, des_nlsw88, quantiles = seq( 0, 1, .05 ), na.rm = TRUE )

# note: most survey commands in R use Inf degrees of freedom by default
# stata generally uses the degrees of freedom of the survey design.
# therefore, while this extended syntax serves to prove a precise replication of stata
# it is generally not necessary.
section_four_one <-
data.frame(
estimate = coef( result.lin ) ,
standard_error = SE( result.lin ) ,
ci_lower_bound =
coef( result.lin ) +
SE( result.lin ) *
qt( 0.025 , degf( subset( des_nlsw88 , !is.na( wage ) ) ) ) ,
ci_upper_bound =
coef( result.lin ) +
SE( result.lin ) *
qt( 0.975 , degf( subset( des_nlsw88 , !is.na( wage ) ) ) )
)
estimate standard_error ci_lower_bound ci_upper_bound
0 0.0000000 0.0000000 0.0000000 0.0000000
0.05 0.0151060 0.0004159 0.0142904 0.0159216
0.1 0.0342651 0.0007021 0.0328882 0.0356420
0.15 0.0558635 0.0010096 0.0538836 0.0578434
0.2 0.0801846 0.0014032 0.0774329 0.0829363
0.25 0.1067687 0.0017315 0.1033732 0.1101642
0.3 0.1356307 0.0021301 0.1314535 0.1398078
0.35 0.1670287 0.0025182 0.1620903 0.1719670
0.4 0.2005501 0.0029161 0.1948315 0.2062687
0.45 0.2369209 0.0033267 0.2303971 0.2434447
0.5 0.2759734 0.0037423 0.2686347 0.2833121
0.55 0.3180215 0.0041626 0.3098585 0.3261844
0.6 0.3633071 0.0045833 0.3543192 0.3722950
0.65 0.4125183 0.0050056 0.4027021 0.4223345
0.7 0.4657641 0.0054137 0.4551478 0.4763804
0.75 0.5241784 0.0058003 0.5128039 0.5355529
0.8 0.5880894 0.0062464 0.5758401 0.6003388
0.85 0.6577051 0.0066148 0.6447333 0.6706769
0.9 0.7346412 0.0068289 0.7212497 0.7480328
0.95 0.8265786 0.0062686 0.8142857 0.8388715
1 1.0000000 0.0000000 1.0000000 1.0000000

For additional usage examples of svylorenz, type ?convey::svylorenz in the R console.