text:
# Solving the chinese postman problem

Some pre-Halloween post today. It started actually while I was in Barcelona : kids wanted to go back to some store we’ve seen the first day, in the gothic part, and I could not remember where it was. And I said to myself that would be quite long to do all the street of the neighborhood. And I discovered that it was actually an old problem. In 1962, Meigu Guan was interested in a postman delivering mail to a number of streets such that the total distance walked by the postman was as short as possible. How could the postman ensure that the distance walked was a minimum?

A very close notion is the concept of traversable graph, which is one that can be drawn without taking a pen from the paper and without retracing the same edge. In such a case the graph is said to have an Eulerian trail (yes, from Euler’s bridges problem). An Eulerian trail uses all the edges of a graph. For a graph to be Eulerian all the vertices must be of even order.

An algorithm for finding an optimal Chinese postman route is:

1. List all odd vertices.
2. List all possible pairings of odd vertices.
3. For each pairing find the edges that connect the vertices with the minimum weight.
4. Find the pairings such that the sum of the weights is minimised.
5. On the original graph add the edges that have been found in Step 4.
6. The length of an optimal Chinese postman route is the sum of all the edges added to the total found in Step 4.
7. A route corresponding to this minimum weight can then be easily found.

For the first steps, we can use the codes from Hurley & Oldford’s Eulerian tour algorithms for data visualization and the PairViz package. First, we have to load some R packages

require(igraph) require(graph) require(eulerian) require(GA)

Then use the following function from stackoverflow,

make_eulerian = function(graph){ info = c("broken" = FALSE, "Added" = 0, "Successfull" = TRUE) is.even = function(x){ x %% 2 == 0 } search.for.even.neighbor = !is.even(sum(!is.even(degree(graph)))) for(i in V(graph)){ set.j = NULL uneven.neighbors = !is.even(degree(graph, neighbors(graph,i))) if(!is.even(degree(graph,i))){ if(sum(uneven.neighbors) == 0){ if(sum(!is.even(degree(graph))) > 0){ info["Broken"] = TRUE uneven.candidates <- !is.even(degree(graph, V(graph))) if(sum(uneven.candidates) != 0){ set.j <- V(graph)[uneven.candidates][[1]] }else{ info["Successfull"] <- FALSE } } }else{ set.j <- neighbors(graph, i)[uneven.neighbors][[1]] } }else if(search.for.even.neighbor == TRUE & is.null(set.j)){ info["Added"] <- info["Added"] + 1 set.j <- neighbors(graph, i)[ !uneven.neighbors ][[1]] if(!is.null(set.j)){search.for.even.neighbor <- FALSE} } if(!is.null(set.j)){ if(i != set.j){ graph <- add_edges(graph, edges=c(i, set.j)) info["Added"] <- info["Added"] + 1 } } } (list("graph" = graph, "info" = info))}

Then, consider some network, with 12 nodes

g1 = graph(c(1,2, 1,3, 2,4, 2,5, 1,5, 3,5, 4,7, 5,7, 5,8, 3,6, 6,8, 6,9, 9,11, 8,11, 8,10, 8,12, 7,10, 10,12, 11,12), directed = FALSE)

To plot that network, use

V(g1)$name=LETTERS[1:12] V(g1)$color=rgb(0,0,1,.4) ly=layout.kamada.kawai(g1) plot(g1,vertex.color=V(newg)$color,layout=ly) Then we convert it to some traversable graph by adding 5 vertices eulerian = make_eulerian(g1) eulerian$info broken Added Successfull 0 5 1 g = eulerian$graph as shown below ly=layout.kamada.kawai(g) plot(g,vertex.color=V(newg)$color,layout=ly)

We cut those 5 vertices in two part, and therefore, we add 5 artificial nodes

A=as.matrix(as_adj(g)) A1=as.matrix(as_adj(g1)) newA=lower.tri(A, diag = FALSE)*A1+upper.tri(A, diag = FALSE)*A for(i in 1:sum(newA==2)) newA = cbind(newA,0) for(i in 1:sum(newA==2)) newA = rbind(newA,0) s=nrow(A) for(i in 1:nrow(A)){ Aj=which(newA[i,]==2) if(!is.null(Aj)){ for(j in Aj){ newA[i,s+1]=newA[s+1,i]=1 newA[j,s+1]=newA[s+1,j]=1 newA[i,j]=1 s=s+1 }}}

We get the following graph, where all nodes have an even number of vertices !

newg=graph_from_adjacency_matrix(newA) newg=as.undirected(newg) V(newg)$name=LETTERS[1:17] V(newg)$color=c(rep(rgb(0,0,1,.4),12),rep(rgb(1,0,0,.4),5)) ly2=ly transl=cbind(c(0,0,0,.2,0),c(.2,-.2,-.2,0,-.2)) for(i in 13:17){ j=which(newA[i,]>0) lc=ly[j,] ly2=rbind(ly2,apply(lc,2,mean)+transl[i-12,]) } plot(newg,layout=ly2)

Our network is now the following (new nodes are small because actually, they don’t really matter, it’s just for computational reasons)

plot(newg,vertex.color=V(newg)$color,layout=ly2, vertex.size=c(rep(20,12),rep(0,5)), vertex.label.cex=c(rep(1,12),rep(.1,5))) Now we can get the optimal path n <- LETTERS[1:nrow(newA)] g_2 <- new("graphNEL",nodes=n) for(i in 1:nrow(newA)){ for(j in which(newA[i,]>0)){ g_2 <- addEdge(n[i],n[j],g_2,1) }} etour(g_2,weighted=FALSE) [1] "A" "B" "D" "G" "E" "A" "C" "E" "H" "F" "I" "K" "H" "J" "G" "P" "J" "L" "K" "Q" "L" "H" "O" "F" "C" [26] "N" "E" "B" "M" "A" or edg=attr(E(newg), "vnames") ET=etour(g_2,weighted=FALSE) parcours=trajet=rep(NA,length(ET)-1) for(i in 1:length(parcours)){ u=c(ET[i],ET[i+1]) ou=order(u) parcours[i]=paste(u[ou[1]],u[ou[2]],sep="|") trajet[i]=which(edg==parcours[i]) } parcours [1] "A|B" "B|D" "D|G" "E|G" "A|E" "A|C" "C|E" "E|H" "F|H" "F|I" "I|K" "H|K" "H|J" "G|J" "G|P" "J|P" [17] "J|L" "K|L" "K|Q" "L|Q" "H|L" "H|O" "F|O" "C|F" "C|N" "E|N" "B|E" "B|M" "A|M" trajet [1] 1 3 8 9 4 2 6 10 11 12 16 15 14 13 26 27 18 19 28 29 17 25 24 7 22 23 5 21 20 Let us try now on a real network of streets. Like Missoula, Montana. I will not try to get the shapefile of the city, I will just try to replicate the photography above. If you look carefully, you will see some problem : 10 and 93 have an odd number of vertices (3 here), so one strategy is to connect them (which explains the grey line). But actually, to be more realistic, we start in 93, and we end in 10. Here is the optimal (shortest) path which goes through all vertices. Now, we are ready for Halloween, to go through all streets in the neighborhood ! # Game of Friendship Paradox In the introduction of my course next week, I will (briefly) mention networks, and I wanted to provide some illustration of the Friendship Paradox. On network of thrones (discussed in Beveridge and Shan (2016)), there is a dataset with the network of characters in Game of Thrones. The word “friend” might be abusive here, but let’s continue to call connected nodes “friends”. The friendship paradox states that People on average have fewer friends than their friends This was discussed in Feld (1991) for instance, or Zuckerman & Jost (2001). Let’s try to see what it means here. First, let us get a copy of the dataset download.file("https://www.macalester.edu/~abeverid/data/stormofswords.csv","got.csv") GoT=read.csv("got.csv") library(networkD3) simpleNetwork(GoT[,1:2]) Because it is difficult for me to incorporate some d3js script in the blog, I will illustrate with a more basic graph, Consider a vertex $v\in V$ in the undirected graph $G=(V,E)$ (with classical graph notations), and let $d(v)$ denote the number of edges touching it (i.e. $v$ has $d(v)$ friends). The average number of friends of a random person in the graph is $$\mu = \frac{1}{n_V}\sum_{v\in V} d(v)=\frac{2 n_E}{n_V}$$ The average number of friends that a typical friend has is $$\frac{1}{n_V}\sum_{v\in V} \left(\frac{1}{d(v)}\sum_{v'\in E_v} d(v')\right)$$But $$\sum_{v\in V} \left(\frac{1}{d(v)}\sum_{v'\in E_v} d(v')\right)=\sum_{v,v' \in G} \left( \frac{d(v')}{d(v)}+\frac{d(v)}{d(v')}\right)$$ $$=\sum_{v,v' \in G}\left(\frac{d(v')^2+d(v)^2}{d(v)d(v')}\right)=\sum_{v,v' \in G} \left(\frac{(d(v')-d(v))^2}{d(v)d(v')}+2\right){\color{red}{\succ}}\sum_{v,v' \in G} \left(2\right)=\sum_{v\in V} d(v)$$ Thus,$$\frac{1}{n_V}\sum_{v\in V} \left(\frac{1}{d(v)}\sum_{v'\in E_v} d(v')\right)\succ \frac{1}{n_V}\sum_{v\in V} d(v)$$ Note that this can be related to the variance decomposition $$\text{Var}[X]=\mathbb{E}[X^2]-\mathbb{E}[X]^2$$i.e.$$\frac{\mathbb{E}[X^2]}{\mathbb{E}[X]} =\mathbb{E}[X]+\frac{\text{Var}[X]}{\mathbb{E}[X]}\succ\mathbb{E}[X]$$(Jensen inequality). But let us get back to our network. The list of nodes is M=(rbind(as.matrix(GoT[,1:2]),as.matrix(GoT[,2:1]))) nodes=unique(M[,1]) and we each of them, we can get the list of friends, and the number of friends friends = function(x) as.character(M[which(M[,1]==x),2]) nb_friends = Vectorize(function(x) length(friends(x))) as well as the number of friends friends have, and the average number of friends friends_of_friends = function(y) (Vectorize(function(x) length(friends(x)))(friends(y))) nb_friends_of_friends = Vectorize(function(x) mean(friends_of_friends(x))) We can look at the density of the number of friends, for a random node, Nb = nb_friends(nodes) Nb2 = nb_friends_of_friends(nodes) hist(Nb,breaks=0:40,col=rgb(1,0,0,.2),border="white",probability = TRUE) hist(Nb2,breaks=0:40,col=rgb(0,0,1,.2),border="white",probability = TRUE,add=TRUE) lines(density(Nb),col="red",lwd=2) lines(density(Nb2),col="blue",lwd=2) and we can also compute the averages, just to check mean(Nb) [1] 6.579439 mean(Nb2) [1] 13.94243 So, indeed, people on average have fewer friends than their friends. # Classification from scratch, neural nets 6/8 Sixth post of our series on classification from scratch. The latest one was on the lasso regression, which was still based on a logistic regression model, assuming that the variable of interest $Y$ has a Bernoulli distribution. From now on, we will discuss technique that did not originate from those probabilistic models, even if they might still have a probabilistic interpretation. Somehow. Today, we will start with neural nets. Maybe I should start with a disclaimer. The goal is not to replicate well designed R functions, used for predictive modeling. It is simply to get a basic understanding of what’s going on. ## Networs, nodes and edges First of all, neurals nets are nets, or networks. I will skip the parallel with “neural” stuff because it does not help me understanding what is happening (all apologies for my poor knowledge on biology, and cells) So, it’s about some network. Networks have nodes, and edges (possibly connected) that connect nodes, or maybe, to more specific (at least it helped me understanding what’s going on), some sort of flow network, In such a network, we usually have sources (here multiple) sources (here $\color{red}\{s_1,s_2,s_3\}$), on the left, on a sink (here $\{\color{blue}t\}$), on the right. To continue with this metaphorical introduction, information from the sources should reach the sink. An usually, sources are explanatory variables, $\{\mathbf{x}_1,\cdots,\mathbf{x}_p\}$, and the sink is our variable of interest $\mathbf{y}$. And we want to create a graph, from the sources to the sink. We will have directed edges, with only one (unique) direction, where we will put weights. It is not a flow, the parallel with flow will stop here. For instance, the most simple network will be the following one, with no layer (i.e no node between the source and the sink) The output here is a binary variable $y\in\{0,1\}$ (it can also be $y\in\{-1,+1\}$ but here, it’s not a big deal). In our network, our output will be $y\in(0,1)$, because it is more easy to handly. For instance, consider $y=f($something$)$, for some function $f$ taking values in $(0,1)$. One can consider the sigmoid function$$f(x)=\frac{1}{1+e^{-x}}=\frac{e^{x}}{e^{x}+1}$$which is actually the logistic function (so we should not be surprised to have results somehow close the logistic regression…). This function $f$ is called the activation function, and there are thousands of such functions. If $y\in\{-1,+1\}$, people consider the hyperbolic tangent$$f(x)=\tanh(x)={\frac {(e^{x}-e^{-x})}{(e^{x}+e^{-x})}}$$or the inverse tangent function $$f(x)=\tan ^{-1}(x)$$And as input for such function, we consider a weighted sum of incoming nodes. So here$$y_i=f\left(\sum_{j=1}^p\omega_j x_{j,i}\right)$$We can also add a constant actually$$y_i=f\left(\omega_0+\sum_{j=1}^p\omega_j x_{j,i}\right)$$So far, we are not far away from the logistic regression. Except that our starting point was a probabilistic model, in the sense that the later was interpreted as a probability (the probability that $Y=1$) and we wanted the model with the highest likelihood. But we’ll talk about selection of weights later one. First, let us construct our first (very simple) neural network. First, we have the sigmoid function sigmoid = function(x) 1 / (1 + exp(-x)) The consider some weights. In our model with seven explanatory variables, with need 7 weights. Or 8 if we include the constant term. Let us consider $\mathbf{\omega}=\mathbf{1}$, weights_0 = rep(1,8) X = as.matrix(cbind(1,myocarde[,1:7])) y_5_1 = sigmoid(X %*% weights_0) that’s kind of stupid because all our predictions are 1, here. Let us try something else. Like $\mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}$. It is optimized, somehow, but we needed something to visualize what’s going on weights_0 = lm(PRONO~.,data=myocarde)$coefficients

then use

y_5_1 = sigmoid(X %*% weights_0)

In order to see if we get a “good” prediction, let use plot the ROC curve, and compare it with the one we got with a (simple) logistic regression

library(ROCR) pred = ROCR::prediction(y_5_1,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) reg = glm(PRONO~.,data=myocarde,family=binomial(link = "logit")) y_0 = predict(reg,type="response") pred0 = ROCR::prediction(y_0,myocarde$PRONO) perf0 = ROCR::performance(pred0,"tpr", "fpr") plot(perf0,add=TRUE,col="red")

That’s not bad for a very first attempt. Except that we’ve been cheating here, since we did use $\mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}$. How, for real, should we choose those weights?

## Using a loss function

Well, if we want an “optimal” set of weights, we need to “optimize” an objective function. So we need to quantify the loss of a mistake, between the prediction, and the observation. Consider here a quadratic loss function

loss = function(weights){ mean( (myocarde$PRONO-sigmoid(X %*% weights))^2) } It might be stupid to use a quadratic loss function for a classification, but here, it’s not the point. We just want to understand what is the algorithm we use, and the loss function $\ell$ is just one parameter. Then we want to solve$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,f(\omega_0+\mathbf{x}_i^T\mathbf{\omega})\right)\right\rbrace$$Thus, consider weights_1 = optim(weights_0,loss)$par

(where the starting point is the OLS estimate). Again, to see what’s going on, let us visualize the ROC curve

y_5_2 = sigmoid(X %*% weights_1) pred = ROCR::prediction(y_5_2,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red") That’s not amazing, but again, that’s only a first step. ## A single layer Let us add a single layer in our network. Those nodes are connected to the sources (incoming from sources) from the left, and then connected to the sink, on the right. Those nodes are not inter-connected. And again, for that network, we need edges (i.e series of weights). For instance, on the network above, we did add one single layer, with (only) three nodes. For such a network, the prediction formula is $$\mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3\omega_h f_h\left(\omega_{h,0}+ \sum_{j=1}^p \omega_{h,j} x_j\right)\right)$$or more synthetically$$\mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)$$Usually, we consider the same activation function everywhere. Don’t ask me why, I find that weird. Now, we have a lot of weights to choose. Let us use again OLS estimates weights_1 <- lm(PRONO~1+FRCAR+INCAR+INSYS+PAPUL+PVENT,data=myocarde)$coefficients X1 = as.matrix(cbind(1,myocarde[,c("FRCAR","INCAR","INSYS","PAPUL","PVENT")])) weights_2 <- lm(PRONO~1+INSYS+PRDIA,data=myocarde)$coefficients X2=as.matrix(cbind(1,myocarde[,c("INSYS","PRDIA")])) weights_3 <- lm(PRONO~1+PAPUL+PVENT+REPUL,data=myocarde)$coefficients X3=as.matrix(cbind(1,myocarde[,c("PAPUL","PVENT","REPUL")]))

In that case, we did specify edges, and which sources (explanatory variables) should be used for each additional node. Actually, here, other techniques could be have been used, like using a PCA. Each node will then be one of the components. But we’ll use that idea later one…

X = cbind(sigmoid(X1 %*% weights_1), sigmoid(X2 %*% weights_2), sigmoid(X3 %*% weights_3))

But we’re not done here. Those were weights from the source to the know nodes, in the layer. We still need the weights from the nodes to the sink. Here, let use use a simple average

weights = c(1/3,1/3,1/3) y_5_3 <- sigmoid(X %*% weights)

Again, we can plot the ROC curve to see what we’ve done…

pred = ROCR::prediction(y_5_3,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red") ## On back propagation Now, we need some optimal selection of those weights. Observe that with only 3 nodes, there are already $(7+1)\times3+3=27$ parameters in that model! Clearly, parcimony is not the major issue when you start using neural nets! If $$p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)$$we want to solve$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,p(\mathbf{x}_i)\right)\right\rbrace$$for some loss function, which is$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i-p(\mathbf{x}_i))^2 \right\rbrace$$for the quadratic norm, or$$\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i\log p(\mathbf{x}_i)+[1-y_i]\log [1-p(\mathbf{x}_i)]) \right\rbrace$$if we want to use cross-entropy. For convenience, let us center all the variable we create, otherwise, we get numerical problems. center = function(z) (z-mean(z))/sd(z) loss = function(weights){ weights_1 = weights[0+(1:7)] weights_2 = weights[7+(1:7)] weights_3 = weights[14+(1:7)] weights_ = weights[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) mean( (myocarde$PRONO-sigmoid(X %*% weights_))^2)}

Now that we have our objective function, consider some starting points. We can consider weights from a PCA, and then use a gradient descent algorithm,

pca = princomp(myocarde[,1:7]) W = get_pca_var(pca)$contrib weights_0 = c(W[,1],W[,2],W[,3],c(-1,rep(1,3)/3)) weights_opt = optim(weights_0,loss)$par

The prediction is then obtained using

weights_1 = weights_opt[0+(1:7)] weights_2 = weights_opt[7+(1:7)] weights_3 = weights_opt[14+(1:7)] weights_ = weights_opt[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) y_5_4 = sigmoid(X %*% weights_)

And as previously, why not plot the ROC curve of that model

pred = ROCR::prediction(y_5_4,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf,add=TRUE,col="red") That’s not too bad. But with 27 coefficients, that’s what we would expect, no? ## Using nnet() function That’s more or less what is done in neural nets functions. Let us now have a look at some dedicated R functions. library(nnet) myocarde_minmax = myocarde minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:7) myocarde_minmax[,j] = minmax(myocarde_minmax[,j]) Here, variables are linearly transformed, to take values in $(0,1)$. Then we can construct a neural network with one single layer, and three nodes, model_nnet = nnet(PRONO~.,data=myocarde_minmax,size=3) summary(model_nnet) a 7-3-1 network with 28 weights options were - b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1 -9.60 -1.79 21.00 14.72 -20.45 -5.05 14.37 -17.37 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 i5->h2 i6->h2 i7->h2 4.72 2.83 -3.37 -1.64 1.49 2.12 2.31 4.00 b->h3 i1->h3 i2->h3 i3->h3 i4->h3 i5->h3 i6->h3 i7->h3 -0.58 -6.03 25.14 18.03 -1.19 7.52 -19.47 -12.95 b->o h1->o h2->o h3->o -1.32 29.00 -10.32 26.27 Here, it is the complete full network. And actually, there are (online) some functions that can he used to visualize that network library(devtools) source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r') plot.nnet(model_nnet) Nice, isn’t it? We clearly see the intermediary layer, with three nodes, and on top the constants. Edges are the plain lines, the darker, the heavier (in terms of weights). ## Using neuralnet() Other R functions can actually be considered. library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet) Again, for the same network structure, with one (hidden) layer, and three nodes in it. ## Network with multiple layers The good thing is that it’s not possible to add more layers. Like two layers. Nodes from the first layer are no longuer connected with the sink, but with nodes in the second layer. And those nodes will then be connected to the sink. We now have something like $$p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_h f_h\left(\omega_{h,0}+ \mathbf{z}_h^T\mathbf{\omega}_h\right)\right)$$where$$\mathbf{z}_h=f\left( \omega_{h,0}+ \sum_{j=1}^{k_h} \omega_{h,j} f_{h,j}\left(\omega_{h,j,0}+ \mathbf{x}^T\mathbf{\omega}_{h,j}\right)\right)$$I may be rambling here (a little bit) but that’s a lot of parameters. Here is the visualization of such a network, library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet) ## Application Let us get back on our simple dataset, with only two covariates. library(neuralnet) df_minmax =df df_minmax$y=(df_minmax$y=="1")*1 minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:2) df_minmax[,j] = minmax(df[,j]) X = as.matrix(cbind(1,df_minmax[,1:2])) Consider only one layer, with two nodes model_nnet = neuralnet(formula(lm(y~.,data=df_minmax)), df_minmax,hidden=c(2)) plot(model_nnet) Here, we did not specify it, but the activation function is the sigmoid (actually, it is called logistic here) model_nnet$act.fct function (x) { 1/(1 + exp(-x)) }   attr(,"type") [1] "logistic" f=model_nnet$act.fct The weights (on the figure) can be obtained using w0 = model_nnet$weights[[1]][[2]][,1] w1 = model_nnet$weights[[1]][[1]][,1] w2 = model_nnet$weights[[1]][[1]][,2]

Now, to get our prediction,
we should use$$p(\mathbf{x})=f\left( \omega_0+ \omega_1 f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_1 f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})\right)$$which can be obtained using

f(cbind(1,f(X%*%w1),f(X%*%w2))%*%w0) [,1] [1,] 0.7336477343 [2,] 0.7317999050 [3,] 0.7185803540 [4,] 0.7404005280 [5,] 0.7518482779 [6,] 0.4939774149 [7,] 0.4965876378 [8,] 0.7101714888 [9,] 0.5050760026 [10,] 0.5049877644

Unfortunately, it is not the output of the model here,

neuralnet::prediction(model_nnet) Data Error: 0; $rep1 x1 x2 y 1 0.1250 0.0000000000 0.02030470787 2 0.0625 0.1176470588 0.89621706711 3 0.9375 0.2352941176 0.01995171956 4 0.0000 0.4705882353 1.10849420363 5 0.5000 0.4705882353 -0.01364966058 6 0.3125 0.5294117647 -0.02409150561 7 0.6875 0.8235294118 0.93743057765 8 0.3750 0.8823529412 1.01320924782 9 1.0000 0.9058823529 1.04805134309 10 0.5625 1.0000000000 1.00377379767 If anyone has a clue, I’d be glad to know what went wrong here… I find that odd to have outputs outside the $(0,1)$ interval, but the output is neither$$p(\mathbf{x})=\omega_{0,0}+ \omega_{0,1} f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_{0,2} f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})$$ cbind(1,f(X%*%w1),f(X%*%w2))%*%w0 [,1] [1,] 1.01320924782 [2,] 1.00377379767 [3,] 0.93743057765 [4,] 1.04805134309 [5,] 1.10849420363 [6,] -0.02409150561 [7,] -0.01364966058 [8,] 0.89621706711 [9,] 0.02030470787 [10,] 0.01995171956 # Traffic Flow of Kota Kinabalu (with R) This morning, we had our first practicals on network flows, using an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map, m=matrix(c(0,5.995910, 116.105520, 1,5.992737, 116.093718, 2,5.992066, 116.109883, 3,5.976947, 116.095760, 4,5.985766, 116.091580, 5,5.988940, 116.080112, 6,5.968318, 116.080764, 7,5.977454, 116.075460, 8,5.974226, 116.073604, 9,5.969651, 116.073753, 10,5.972341, 116.069270, 11,5.978818, 116.072880),3,12) we can be visualized below library(OpenStreetMap) map = openmap(c(lat= 6.000, lon= 116.06), c(lat= 5.960, lon= 116.12)) map=openproj(map) plot(map) points(t(m[3:2,]),col="black", pch=19, cex=3 ) text(t(m[3:2,]),c("s",1:10,"t"),col="white") If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least). To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper library(devtools) install_github("ropensci/tabulizer") library(tabulizer) location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf' out <- extract_tables(location) with Windows, it seems to be necessary to download another package first library(devtools) install_github("ropensci/tabulizerjars") install_github("ropensci/tabulizer") library(tabulizer) location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf' out <- extract_tables(location) Now we can get out data frame with capacities B1=as.data.frame(out[[2]]) B2=as.data.frame(out[[3]]) E=data.frame(from=B1[3:20,"V3"], to=B1[3:20,"V4"]) E=E[-c(6,8),] capacity=as.character(B2$V3[-1]) capacity[6]="843" capacity[4]="2913" E$capacity=as.numeric(capacity) We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read) plot(map) points(t(m[3:2,]),col="black", pch=19, cex=3 ) B=data.frame(i=as.character(c("s",paste("V",1:10,sep=""),"t")), x=m[3,],y=m[2,]) for(i in 1:nrow(E)){ i1=which(B$i==as.character(E$from[i])) i2=which(B$i==as.character(E$to[i])) segments(B[i1,"x"],B[i1,"y"],B[i2,"x"],B[i2,"y"],lwd=3) } text(t(m[3:2,]),c("s",1:10,"t"),col="white") To get the graph with capacities, an alternative is to use library(igraph) g=graph_from_data_frame(E) E(g)$label=E$capacity plot(g) but it does not respect geographical locations of nodes. It can actually be done using plot(g, layout=as.matrix(B[,c("x","y")])) To get a better understanding of the capacities of the road, use plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$capacity/200)

From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using

> (m=max_flow(graph=g, source="s", target="t")) $value [1] 2571 $flow [1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422 [16] 1149

Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers  max flow min cut theorem to… and application of the Shortest Path… (“the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables.

E$flux1=m$flow E(g)$label=E$flux1 plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$flux1/200) That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value E$flux2=c(1422,1149,1422,1149,0,0,0,0, 1149,1422,1149,0,0,1149,1422,1149) E(g)$label=E$flux2 plot(g, layout=as.matrix(B[,c("x","y")]), edge.width=E$flux2/200) Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu. location <- 'http://www.worldresearchlibrary.org/up_proc/pdf/999-150486366625-30.pdf' out <- extract_tables(location) dim(out[[3]]) B1=as.data.frame(out[[3]]) E=data.frame(from=B1[2:61,"V2"], to=B1[2:61,"V3"], capacity=B1[2:61,"V4"]) E$capacity=as.numeric( as.character(E$capacity)) library(igraph) g=graph_from_data_frame(E) m=max_flow(graph=g, source="S", target="T") E$flux1=m$flow E(g)$label=E$flux1 plot(g, edge.width=E$flux1/200, edge.arrow.size=0.15)

Here the value of the maximal flow is 4017, just as they found in the original paper

# Traffic Flow of Kota Kinabalu (Malaysia)

For the second practicals of our course on networks and flows, we will study traffic flow of Kota Kinabalu (Malaysia), following several papers published by Noraini Abdullah and Ting Kien Hua, such as max flow min cut theorem to minimize traffic congestion in Kota Kinabalu, traffic congestion problem of road networks in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu.

# Networks and Flows #2

After the Traveling Salesman part, we will see tools used to study flows in networks this week. Slides are now online (from slide 42).

Slides are probably full of typos, not to say mistakes. All apologies.

I will also include some slides on the simplex method (since I am not sure that all students have seen it)

# Métro: centralité et robustesse

Demain matin, nous aurons un TP pour le cours sur les réseaux, les flux et les transports. En particulier, en nous inspirant des travaux de Sybil Derrible, nous allons commencer par étudier la centralité dans les différents systèmes de métro, mais aussi la robustesse. Les matrices d’adjacence d’une trentaine de métros dans le monde sont en ligne dans un fichier xls. Histoire de gagner un peu de temps, le code pour créer une matrice d’adjacence peut être le suivant

loc="/data/Metro_Networks_Adjacency.xls" library(xlsx) E=read.xlsx(loc,"StPetersburg") n=nrow(E) nom=as.character(E[3:(n-2),1]) Adj=E[3:(n-2),(4:ncol(E)-1)] Adj[is.na(Adj)]=0 Adj=as.matrix(Adj) colnames(Adj)=rownames(Adj)=nom

On est ensuite prêt à manipuler le réseau,

library(igraph) iflo=graph_from_adjacency_matrix(Adj,mode = "undirected") plot(iflo)

On va utiliser les notions vues en cours, sur la centralité, mais surtout, on travaillera sur Quantifying the robustness of metro networks, inspiré de The complexity and robustness of metro networks. Plusieurs fonctions utiles sont déjà programmées dans R, comme l’assortativité.

# Networks with R

In order to practice with network data with R, we have been playing with the Padgett (1994) Florentine’s wedding dataset (discussed in the lecture). The dataset is available from

> library(network) > data(flo) > nflo=network(flo,directed=FALSE) > plot(nflo, displaylabels = TRUE, + boxed.labels = + FALSE)

The next step was to move from the network package to igraph. Since we have the adjacency matrix, we can use it

> library(igraph) > iflo=graph_from_adjacency_matrix(flo, + mode = "undirected") > plot(iflo)

The good thing is that a lot of functions are available, for instance we can get shortest paths, between two specific nodes. And we can give appropriate colors to the nodes that we’ll cross

> AP=all_shortest_paths(iflo, + from="Peruzzi", + to="Ginori") > L=AP$res[[1]] > V(iflo)$color="yellow" > V(iflo)$color[L[2:4]]="light blue" > V(iflo)$color[L[c(1,5)]]="blue" > plot(iflo)

We can also visualize edges, but I found it slightly more complicated (to extract edges from the output)

> liens=c(paste(as.character(L)[1:4], + "--", + as.character(L)[2:5],sep=""), + paste(as.character(L)[2:5], + "--", + as.character(L)[1:4],sep="")) > df=as.data.frame(ends(iflo,E(iflo))) > names(df)=c("src","target") > lstn=sort(unique(c(as.character(df[,1]),as.character(df[,2]),"Pucci"))) > Eliens=paste(as.numeric(factor(df[,1],levels=lstn)),"--", + as.numeric(factor(df[,2],levels=lstn)),sep="") > EU=unlist(lapply(Eliens,function(x) x%in%liens)) > E(iflo)$color=c("grey","black")[1+EU] > plot(iflo) But it works. It is also possible to use some D3js visualization > library( networkD3 ) > simpleNetwork (df) Then the next question was to add a vertice to the network. The most simple way to do it is probability through the adjacency matrix > flo2=flo > flo2["Pucci","Bischeri"]=1 > flo2["Bischeri","Pucci"]=1 > nflo2=network(flo2,directed=FALSE) > plot(nflo2, displaylabels = TRUE, + boxed.labels = + FALSE) Then, we’ve been playing with centrality measures. > plot(iflo,vertex.size=betweenness(iflo)) The goal was to see how related they were. Here, for all of them, “Medici” is the central node. But what about the others? > B=betweenness(iflo) > C=closeness(iflo) > D=degree(iflo) > E=eigen_centrality(iflo)$vector > base=data.frame(betw=B,close=C,deg=D,eig=E) > cor(base) betw close deg eig betw 1.0000000 0.5763487 0.8333763 0.6737162 close 0.5763487 1.0000000 0.7572778 0.7989789 deg 0.8333763 0.7572778 1.0000000 0.9404647 eig 0.6737162 0.7989789 0.9404647 1.0000000

Those measures are quite correlated. It is also possible to use a hierarchical graph to visualize how close those centrality measures can be

> H=hclust(dist(t(base)), + method="ward") > plot(H)

Instead of looking at values of centrality measures, it is possible to looks are ranks

> rbase=base > for(i in 1:4) rbase[,i]=rank(base[,i]) > H=hclust(dist(t(rbase)), + method="ward") > plot(H)

Here the eigenvector measure is very close to the degree of vertices.

Finally, it is possible to seek clusters (in the context of coalition here, in case a war should start between those families)

> kc <- fastgreedy.community ( iflo )

Here we have 3 classes (+1 for the node that is disconnected from the other families)

> V(iflo)\$color=c("yellow","orange", + "light blue")[membership ( kc )] > plot(iflo)

> plot(kc,iflo)


text:
Questions & Answers
Question
Answers

# How many cubic blocks of wood of side 20cm can be cut from a block of wood having dimensions of 2m, 80cm and 40cm?

Answer
Verified
129.3k+ views
Hint: This is the formula for the volume of a block of wood having different dimensions.
$Volume{\text{ }} = {\text{ }}length\; \times \;width\; \times height$
And volume of a block of wood having same dimensions is $Volume{\text{ }} = {\text{ }}length\; \times \;length\; \times length$
And after calculating the volumes of cubes and cuboid, one can check the number of cubes present in cuboids.

Complete step by step solution: When you look around the room in your house, you can see many objects that have the shape of a rectangular box. For example, furniture, books and TV are in the shape of a rectangular box. This shape is a cuboid. On the other hand, objects such as ice, dice and Rubik’s cube are examples of another 3D shape called a cube. As a matter of fact, a cube is a special type of cuboid in which all sides are squares and identical. Let us do it for cubes and cuboids. Here in this question, the bigger block of wood will be acting as a cuboid and smaller ones will be acting as cubes.
As we know, there are different lengths of all the sides of the given bigger wooden block, so it will be having the shape of a cuboid. Therefore,
$\begin{array}{*{20}{l}} {Length{\text{ }}of{\text{ }}Wooden{\text{ }}Block{\text{ }} = {\text{ }}2m{\text{ }} = {\text{ }}2{\text{ }} \times {\text{ }}100{\text{ }} = {\text{ }}200cm} \\ {Width{\text{ }}of{\text{ }}Wooden{\text{ }}Block{\text{ }} = {\text{ }}80cm} \\ {Height{\text{ }}of{\text{ }}Wooden{\text{ }}Block{\text{ }} = {\text{ }}40cm} \end{array}$
$\begin{array}{*{20}{l}} {Volume{\text{ }}of{\text{ }}block{\text{ }}of{\text{ }}wood{\text{ }} = {\text{ }}200cm{\text{ }} \times {\text{ }}80cm{\text{ }} \times {\text{ }}40cm{\text{ }} = {\text{ }}64000cm{\text{ }}cube} \\ {Volume{\text{ }}of{\text{ }}1{\text{ }}cubical{\text{ }}block{\text{ }} = {\text{ }}20cm{\text{ }} \times {\text{ }}20cm{\text{ }} \times {\text{ }}20cm{\text{ }} = {\text{ }}800cm{\text{ }}cube} \\ {No.{\text{ }}of{\text{ }}blocks{\text{ }} = 64000 \div 800{\text{ }} = {\text{ }}80{\text{ }}blocks} \end{array}$

Note: The total number of given smaller wooden blocks out of the bigger blocks can be calculated by dividing the volume of the bigger block with the volume of the smaller one.

format: plain


text:
It looks like the Mythbusters tend to ignore air resistance.

In a recent episode, they claimed to have demonstrated that a horizontally fired bullet and a bullet that is simply dropped fall to the ground in the same amount of time. They were wrong. What they actually demonstrated is that air resistance causes the fired bullet to hit the ground more slowly.

Their argument would apply perfectly in a vacuum, as on the surface of the Moon, but not here on the Earth, where the bullet’s motion is governed not just by the laws of gravity, but also by the laws of a non-conservative force, namely air resistance. (Why is it non-conservative? Some of the bullet’s kinetic energy is converted into heat, as it travels through the air at high speed. Unless we also include the thermodynamics of the air into our equations of motion, the equations will not conserve energy, as the amount of kinetic energy converted into heat will just appear “lost”.)

The bullet’s velocity, $v$, can be written as $v^2=v_h^2+v_v^2$, where $v_h$ is the horizontal and $v_v$ is the vertical component. The initial horizontal velocity is $v_0$. The initial vertical velocity is 0.

Air resistance is proportional to the square of the bullet’s velocity. To be precise, acceleration due to air resistance will be

$a_\mathrm{air}=\kappa v^2,$

where $\kappa$ is an unknown proportionality factor. (To be more precise, $\kappa=\frac{1}{2}cA\rho m$, where $c$ is the dimensionless drag coefficient, $A$ is the bullet’s cross-sectional area, $\rho$ is the density of the air, and $m$ is the bullet’s mass.) The direction of the acceleration will be opposite the direction of the bullet’s motion.

The total acceleration of the bullet will have two components: a vertical component $g$ due to gravity, and a component $a_\mathrm{air}$ opposite the direction of the bullet’s motion.

The resulting equations of motion can be written as:

\begin{align}\frac{dv_h}{dt}&=-\kappa vv_h,\\\frac{dv_v}{dt}&=g-\kappa vv_v.\end{align}

Right here we can see the culprit: air resistance not only slows the bullet down horizontally, it also reduces its downward acceleration.

This is a simple system of two differential equations in the two unknown functions $v_h(t)$ and $v_v(t)$. Its solution is not that simple, unfortunately. However, it can be greatly simplified if we notice that given that $v_v\ll v_h$, $v\simeq v_h$, and therefore, we get

\begin{align}\frac{dv_h}{dt}&=-\kappa v_h^2,\\\frac{dv_v}{dt}&=g-\kappa v_vv_h.\end{align}

This system is solved by

\begin{align}v_h&=\frac{1}{\kappa t+C_1},\\v_v&=v_h\left[\left(\frac{1}{2}\kappa t^2+C_1t\right)g+C_2\right].\end{align}

Given

$v_0=\frac{1}{\kappa t_0+C_1},$

we have

$C_1=\frac{1}{v_0-\kappa t_0},$

$v_h=\frac{v_0}{\kappa(t-t_0)v_0+1},$

or, if we set $t_0=0$,

$v_h=\frac{v_0}{\kappa tv_0+1}.$

Similarly, given $v_v(0)=0$, we get $C_2=0$, thus

$v_v=gt\frac{\kappa v_0t+2}{2\kappa v_0t+2}.$

The distance traveled horizontally ($s_h$) and vertically ($s_v$) between $t_0=0$ and $t_1$ can be obtained by simple integration of the respective velocities with respect to $t$ between 0 and $t_1$:

\begin{align}s_h&=\frac{\log(\kappa v_0t_1+1)}{\kappa},\\s_v&=g\frac{\kappa v_0t_1(\kappa v_0t_1+2)-2\log(\kappa v_0t_1+1)}{(2\kappa v_0)^2}.\end{align}

The claim by the Mythbusters was that the time it took for the fired bullet to hit the ground was only $\sim 40~{\rm ms}$ more than the time it took for a dropped bullet to fall, which is a negligible difference. But it is not! Taking $g\simeq 10~{\rm m}/{\rm s}^2$, it is easy to see that the time it takes for a bullet to fall from a height of $1~{\rm m}$, using the well-known formula $\frac{1}{2}gt^2$, is $447~{\rm ms}$; the difference measured by the Mythbusters is nearly 10% of this number!

Not only did the fired bullet take longer to hit the ground, the Mythbusters’ exquisite setup allows us to calculate the bullet’s initial velocity $v_0$ and drag coefficient $\kappa$. This is possible because the Mythbusters conveniently provided three pieces of information (I am using approximate numbers here): the length of the path that the bullet traveled $s_h\simeq 100~{\rm m}$, the height of the bullet at the time of firing ($s_v\simeq 1~{\rm m}$), and the time it took for the fired bullet to hit the ground. Actually, what they provided was the difference between the time for a fired vs. a dropped bullet to hit the ground, but we know what it is for the dropped bullet (and because it is never moving very rapidly, we can ignore air resistance in its case), so $t_1=447+40=487 {\rm ms}$. The solution is given by

\begin{align}\kappa&=0.0054~\mathrm{m}^{-1},\\v_0&=272.3~\mathrm{m}/\mathrm{s}.\end{align}

Given a bullet cross-sectional area of $A=2~{\rm cm}^2=2\times 10^{-4}~{\rm m}^2$, an approximate air density of $\rho=1~{\rm kg}/{\rm m}^3$, and a bullet mass of $m=20~{\rm g}=0.02~{\rm kg}$, the dimensionless drag coefficient for the bullet can be calculated as $c=2\kappa mA\rho=1.08$, which is not at all unreasonable for a tumbling bullet. Of course the actual values of $A$ and $m$ may differ from the ones I’m using here, resulting in a different value for the dimensionless drag coefficient $c$.

And now (May 28, 2012) I feel obliged to include a footnote, given the almost hysterical news coverage in the last couple of days about an exact solution to this same problem obtained by a 16-year old Dresden youth of Indian descent, Shouryya Ray.

The news coverage was breathless but incomplete: the youth, we were told, solved a problem that baffled Newton 350 years ago, but there was no description of the actual problem (other than a hint that it had to do with air resistance) and its solution, except for a press photo with him holding up a large sheet of paper with an equation.

I suspected that his variables $u$ and $v$ were in fact the horizontal and vertical components of a projectile, which I denoted with $v_h$ and $v_v$ above. So Ray's equation would read, in my notation:

$\frac{g^2}{2v_h^2}-\frac{\kappa g}{2}\left(\frac{v_v\sqrt{v_h^2+v_v^2}}{v_h^2}+{\rm arcsinh}\left|\frac{v_v}{v_h}\right|\right)={\rm const,}$

where I picked up an extra minus sign because I considered $g$ to be positive when pointing downwards.

But how did he find this neat form? Is it even correct? To find out, I tried a few ways to integrate the equations but to no avail. I admit I was ready to give up when I came across a solution posted on reddit. The derivation is actually quite simple. We begin with (I now use the overdot for the time derivative):

\begin{align}
\dot{v}_h=-\kappa\sqrt{v_h^2+v_v^2}v_h,\\
\dot{v}_v=g-\kappa\sqrt{v_h^2+v_v^2}v_v.
\end{align}

Now multiply the first equation by $\dot{v}_v$ and the second by $\dot{v}_h$:

\begin{align}
\dot{v}_h\dot{v}_v=-\kappa\sqrt{v_h^2+v_v^2}v_h\dot{v}_v,\\
\dot{v}_h\dot{v}_v=g\dot{v}_h-\kappa\sqrt{v_h^2+v_v^2}\dot{v}_hv_v.
\end{align}

Subtract one from the other and rearrange:

$g\dot{v}_h=\kappa\sqrt{v_h^2+v_v^2}(\dot{v}_hv_v-v_h\dot{v}_v).$

Now substitute $v_v=sv_h$:

$g\dot{v}_h=\kappa\sqrt{(s^2+1)v_h^2}(\dot{v}_hv_v-v_h\dot{v}_v).$

Divide both sides by $v_h^3$:

$g\frac{\dot{v}_h}{v_h^3}=-\kappa\sqrt{(s^2+1)}\dot{s}.$

In this form, both sides are full derivatives with respect to $t$ and thus both sides can be integrated, to yield

$-\frac{1}{2}\frac{g}{v_h^2}=-\kappa\left[\frac{1}{2}s\sqrt{s^2+1}+\frac{1}{2}{\rm arcsinh}~s\right]+C,$

where $C$ is an integration constant. Multiplying by $g$ and replacing $s$ with $v_v/v_h$, we get back Ray's equation.

This is an elegant result. Its practical utility may be limited, however, as the solution is implicit. Actually computing $v_h$ and $v_v$ as functions of time still requires numerical methods, so one might as well just solve the original differential equations numerically. Perhaps this explains why Ray's formula is not better known. But it has been discovered before. Apart from trivial notational differences, this very formulation is discussed by Parker (Am. J. Phys, 45, 7, 606, July 1977). But given its importance to ballistics, it should come as no surprise that this problem is discussed in depth in textbooks dating back to the middle of the 19th century, for instance in an 1860 book by Didion.


text:
# How to decrypt unusual Many Time Pad

How does one decrypt the many time pad:

($R$ is a random string, $C$ is ciphertext, $P$ is plaintext)

$R_0 \oplus R_1 = C_0$

$R_2 \oplus R_1 = C_1$

$P_0 \oplus R_1 = C_2$

assuming that $R_n$ are perfectly random pads, plaintext is writing and $C_n$ are the encrypted pads. The key $R_1$ is used three times but it doesn't seem like Crib Dragging (the solution I've seen described to decrypt many time pads) would work (edit: to retrieve $R_0$, $R_2$, and $P_0$ or even just $P_0$) since the other values that $R_1$ is added to are random. Is there something I'm missing?

• Hint: rewrite "random[0] xor random[1] = cipher[0]" as "(all-zero plaintext) xor (random[0] xor random[1]) = cipher[0]". is this really a many-time pad? Do cipher[0] and cipher[1] actually carry any information, or are they in fact statistically independent from random[1]? Write down the probabilities on paper to make sure. – Thomas May 3 '15 at 5:13
• Adding an all zero plaintext makes no sense. $R_0 and$R_2 are the information being carried by $C_0$ and $C_1$ and the goal is to retrieve them and $P_0$ or failing that retrieve just $P_0$. Since $R_0$, $R_1$, and $R_2$ are just theoretical values defined as random the probabilities should have a perfectly even distribution, $R_1$ should be statistically independent from $C_0$. – Daedalus6394 May 3 '15 at 20:58
• Both $C_0$ and $C_1$ serve no purpose at all. In order to decrypt and retrieve $P_0$, all you need is $R_1$. So if we assume the key contains $R_1$, I can decrypt even if I only get $C_2$ and not $C_1,C_2$. Sending random things alongside actual information does not mean it is any different than the original scheme. – tylo May 4 '15 at 11:15
• It is assumed in the one time pad that the key, in this case $R_1$ is perfectly secret and unknown. In this case that is extended to $R_n$. Note this is not a question about decryption by the intended receiver but decryption by an unintended third party. – Daedalus6394 May 4 '15 at 16:51

## 1 Answer

This cannot be used to compromise $P_0$. Suppose that, instead of generating $C_0$ and $C_1$ as $R_0\oplus R_1$ and $R_2\oplus R_1$, the attacker instead picked $C_0$ and $C_1$ at random from the set of bit strings the length of $C_2$. This situation is indistinguishable from yours; in both cases, $R_0$ and $R_2$ are independent and follow a uniform probability distribution, and the same is true for $C_0$ and $C_1$. If you can compromise $P_0$ in your case, you can compromise it if the sender had actually generated $C_0$ and $C_1$. But that would mean that sending two completely random strings alongside an OTP message somehow helps compromise the message. Obviously, that's not true; the attacker could otherwise create his own random strings and break the OTP that way.


text:
<img src="https://d5nxst8fruw4z.cloudfront.net/atrk.gif?account=iA1Pi1a8Dy00ym" style="display:none" height="1" width="1" alt="" />
Due to system maintenance, CK-12 will be unavailable on Friday,8/19/2016 from 6:00p.m to 10:00p.m. PT.

# 8.3: Multimedia Resources for Chapter 8

Difficulty Level: At Grade Created by: CK-12

Copy and distribute the lesson worksheets. Ask students to complete the worksheets alone or in pairs as a review of lesson content.

## The Periodic Table and Electron Configuration Worksheet

When Mendeleev organized the periodic table, he placed the elements in vertical columns according to their chemical behavior. That is, elements were placed in the same vertical columns because they behaved similarly in chemical reactions. All the alkali metals (Li,Na,K,Rb,Cs)\begin{align*}(Li, Na, K, Rb, Cs)\end{align*} react with water to produce heat, hydrogen gas, and the metal hydroxide in solution. Essentially, the only difference in the reactions is that the larger alkali metals react faster than the smaller ones. The vertical columns of elements are frequently referred to chemical “families” because of their similar chemical characteristics.

When quantum theory generated electron configurations which demonstrated that the elements in the same family have the same outer energy level electron configuration, the reason these elements behaved similarly became clear. Since chemical behavior is determined by outer energy level electron configuration, it was clear that elements that behaved similarly should have similar electron configuration.

The Electron Configuration of Family 1A Elements
Element Electron Configuration
Li 1s22s1\begin{align*}1s^22s^1\end{align*}
Na 1s22s2sp63s1\begin{align*}1s^22s^2sp^63s^1\end{align*}
K 1s22s2sp63s23p64s1\begin{align*}1s^22s^2sp^63s^23p^64s^1\end{align*}
Rb 1s22s2sp63s23p64s23d104p65s1\begin{align*}1s^22s^2sp^63s^23p^64s^23d^{10}4p^65s^1\end{align*}
Cs 1s22s2sp63s23p64s23d104p65s24d105p66s1\begin{align*}1s^22s^2sp^63s^23p^64s^23d^{10}4p^65s^24d^{10}5p^66s^1\end{align*}
The Electron Configuration of Family 7A Elements
Element Electron Configuration
F 1s22s22p5\begin{align*}1s^22s^22p^5\end{align*}
Cl 1s22s2sp63s23p5\begin{align*}1s^22s^2sp^63s^23p^5\end{align*}
Br 1s22s2sp63s23p64s23d104p5\begin{align*}1s^22s^2sp^63s^23p^64s^23d^{10}4p^5\end{align*}
I 1s22s2sp63s23p64s23d104p65s24d105p5\begin{align*}1s^22s^2sp^63s^23p^64s^23d^{10}4p^65s^24d^{10}5p^5\end{align*}

Exercises

1. If the outermost energy level electron configuration of an atom is ns2np1\begin{align*}ns^2np^1\end{align*},
1. to which family does it belong?____________
2. is the atom a metal, metalloid, non-metal, or a noble gas? ____________
3. how many valence electrons does it have? ____________
2. If the outermost energy level electron configuration of an atom is ns2np4\begin{align*}ns^2np^4\end{align*},
1. to which family does it belong? ____________
2. is the atom a metal, metalloid, non-metal, or a noble gas? ____________
3. how many valence electrons does it have? ____________
3. If the outermost energy level electron configuration of an atom is ns2np6\begin{align*}ns^2np^6\end{align*},
1. to which family does it belong? ____________
2. is the atom a metal, metalloid, non-metal, or a noble gas? ____________
3. how many valence electrons does it have? ____________
4. The electron configuration of an element is [Ar]4s23d3\begin{align*}[Ar]4s^23d^3\end{align*}.
1. What is the identity of the element? ____________
2. In what period does the element belong? ____________
3. In what group does the element belong? ____________
4. Is the element a main group element, a transition element, a lanthanide, or an actinide? ____________
5. Write the electron configuration of only the outermost energy level for an element that is in family 5A of the fifth period of the periodic table. ____________
6. Write the electron configuration of only the outermost energy level for an element that is in family 8A of the third period of the periodic table. ____________

• The worksheet answer keys are available upon request. Please send an email to teachers-requests@ck12.org to request the worksheet answer keys.

### Notes/Highlights Having trouble? Report an issue.

Color Highlighted Text Notes

Show Hide Details
Description
Tags:
Subjects:
Date Created:
Aug 18, 2012
Sep 03, 2015
Image Detail
Sizes: Medium | Original

CK.SCI.ENG.TE.2.Chemistry.8.3
Here


text:
News

# China Sets Renewable Energy Quota and Pens German Agreement

January 16, 2006 by Jeff Shepard

China has told power companies that 5 percent of their electricity will have to come from renewable energy sources by 2010, as the country tries to diversify away from fossil fuels to power its fast-growing economy, a news report said Friday. Zhang Guobao, vice-minister of China's top economic planning body, said the quota would increase to 10 percent by 2020, according to the China Daily newspaper. Wind and solar power are included but nuclear and hydropower are not, Zhang said.

In a related development, China and Germany will jointly invest about €50 million (US$60 million dollars) to build a wind power plant in eastern China. Five wind turbines with a power capacity of 5 megawatts each will be constructed off the coast of the eastern city of Qingdao, generating power when the city hosts the water events of the 2008 Beijing Olympics, the China Daily newspaper said. The two countries have signed a memorandum of understanding about the wind turbine project, as well as six other renewable energy projects, during the Sino-German Forum on Environment in Qingdao. The total investment for all seven projects was about €100 million (US$121 million).

The renewable energy quota applies to large power companies with an installed capacity of more than 5 gigawatts. Analysts estimate that China has 15 such companies generating about half of the country's power supply. China's larger power companies currently rely mainly on coal and other fossil fuels with renewable energy making up a small percentage.

Zhang said China wanted to reduce its reliance on polluting fossil fuels and secure a reliable energy supply for the booming economy, according to the paper. Financial incentives are also offered to power grids that use energy from renewable sources, Zhang said.

Two of China's top power companies, Datang International Power Generation Co. Ltd. and China Huaneng Group, have both invested in wind turbines. Datang said it plans to cut its coal-fired power generation from the current 99 percent to 75 percent by 2014, according to the report.

Total installed capacity of all of China's power plants last year reached 508 gigawatts, up nearly 15 percent from the previous year. It is expected to exceed 1,000 gigawatts within 15 years.


text:
#### Archived

This topic is now archived and is closed to further replies.

# DirectX + Radeon 7500 = No Fog?

## Recommended Posts

Goodlife    122
Hi all, I''ve written a game (written on a Radeon 9500). It tests good on several NVidia cards, and on a Radeon 8500. Then I ran it on a mobile laptop with a Radeon 7500, and I have a funny problem: Instead of fogging (using pixel fog) the display simply cuts off at the edge of visibility, just as though I had fog turned off. I downloaded the latest drivers, double checked to make sure my projection matrix was W compliant, and still no goods. The code I use to set up the fog is here: mDevice->SetRenderState(D3DRS_FOGENABLE,true); Device->SetRenderState(D3DRS_FOGTABLEMODE, D3DFOG_LINEAR); // theNearValue=2000 // theFarValue=4000 (edge of view range) Device->SetRenderState(D3DRS_FOGSTART,*(DWORD *)(&theNearValue)); Device->SetRenderState(D3DRS_FOGEND,*(DWORD *)(&theFarValue)); // theColor=0xFFFFFFFF; Device->SetRenderState(D3DRS_FOGCOLOR,theColor A bit of reading turned up that you''re ''supposed'' to use a value from 0.0 to 1.0 for the FOGSTART/FOGEND values. That''s in DirectX9, not sure if it''s the same in DirectX8.1, which is what I''m using. Anyway, I tried it with a 0.0 thru 1.0 value setup, and still no fog. Does anyone know what the problem might be here? I know the Radeon 7500 card supports fog, because I play other games with fog on it. I assume there''s some sort of default that is different on this card, or something I need to set up. However, all the directx samples I''ve seen show only these renderstates getting set. Thanks in advance for any help

##### Share on other sites
Chaucer    122
I''m not sure what your problem is but I have a radeon 7200 and I''m able to render vertex fog like this:

  m_pd3dDevice->SetRenderState(D3DRS_FOGENABLE, TRUE);	m_pd3dDevice->SetRenderState(D3DRS_FOGCOLOR,  fogColor);    m_pd3dDevice->SetRenderState(D3DRS_FOGSTART,   *(DWORD *) (&fogStart));    m_pd3dDevice->SetRenderState(D3DRS_FOGEND,     *(DWORD *) (&fogEnd));    m_pd3dDevice->SetRenderState(D3DRS_FOGDENSITY, *(DWORD *) (&fogDensity));	m_pd3dDevice->SetRenderState(D3DRS_FOGVERTEXMODE, fogMode);		m_pd3dDevice->SetRenderState(D3DRS_RANGEFOGENABLE, TRUE);

##### Share on other sites
- DX8.1 docs on fogstart:
Depth at which pixel or vertex fog effects begin for linear fog mode. The default value is 0.0f. Depth is specified in world space for vertex fog, and either device space [0.0, 1.0] or world space for pixel fog. For pixel fog, these values are in device space when the system uses z for fog calculations, and world space when the system is using eye-relative fog (w-fog). For more information, see Fog Parameters and Eye-Relative vs. Z-Based Depth.
Values for this render state are floating-point values. Because the IDirect3DDevice8::SetRenderState method accepts DWORD values, your application must cast a variable that contains the value, as shown in the following code example.

pd3dDevice8->SetRenderState(D3DRS_FOGSTART, *((DWORD*) (&fFogStart)));

- If that doesn''t fix the problem, make sure the card supports it (it should, of course, but maybe the drivers don''t expose that).

- If you can try the MFCFog sample on the card, do so.

To Chaucer: He''s using pixel (table) fog, not vertex fog.

Cheers,

##### Share on other sites
Infinisearch    2971
I can''t help you with your problem but i can suggest going to www.rage3d.com they are an ati related hardware site and have a programming forum where people who work for ati frequent. Well at least last time i was there. good luck.

##### Share on other sites
JoeyBlow2    100
They don''t support pixel/table fog. It doesn''t have the CAP... Check the cap then go over to vertex fog (unless that doesn''t work good)

##### Share on other sites
Goodlife    122
Hi Joey,
Thanks, that was it.

Argh. My meshes aren''t very tesselated (the game doesn''t need it) and vertex fog looks HORRIBLE. What a mess. Oh well, I''ll figure out something.


text:
# What could be the possible future applications for HHL algorithm?

Note on the vocabulary: the word "hamiltonian" is used in this question to speak about hermitian matrices.

The HHL algorithm seems to be an active subject of research in the field of quantum computing, mostly because it solve a very important problem which is finding the solution of a linear system of equations.

According to the original paper Quantum algorithm for solving linear systems of equations (Harrow, Hassidim & Lloyd, 2009) and some questions asked on this site

the HHL algorithm is limited to some specific cases. Here is a summary (that may be incomplete!) of the characteristics of the HHL algorithm:

# HHL algorithm

The HHL algorithm solves a linear system of equation $$A \vert x \rangle = \vert b \rangle$$ with the following limitations:

## Limitations on $\vert x \rangle$ (output):

• $\vert x \rangle$ cannot be recovered fully by measurement. The only information we can recover from $\vert x \rangle$ is a "general information" ("expectation value" is the term employed in the original HHL paper) such as $$\langle x\vert M\vert x \rangle$$

Question: Taking into account all of these limitations and imagining we are in 2050 (or maybe in 2025, who knows?) with fault-tolerant large-scale quantum chips (i.e. we are not limited by the hardware), what real-world problems could HHL algorithm solve (including problems where HHL is only used as a subroutine)?

I am aware of the paper Concrete resource analysis of the quantum linear system algorithm used to compute the electromagnetic scattering cross section of a 2D target (Scherer, Valiron, Mau, Alexander, van den Berg & Chapuran, 2016) and of the corresponding implementation in the Quipper programming language and I am searching for other real-world examples where HHL would be applicable in practice. I do not require a published paper, not even an unpublished paper, I just want to have some examples of real-world use-cases.

EDIT:

Even if I am interested in every use-case, I would prefer some examples where HHL is directly used, i.e. not used as a subroutine of an other algorithm.

I am even more interested in examples of linear systems resulting of the discretisation of a differential operator that could be solved with HHL.

But let me emphasise one more time I'm interested by every use-case (subroutines or not) you know about.

• You mention that you want some examples where HHL is "directly used". I am not very clear on what you mean by that. I do know some algorithms (which can potentially have practical uses) in which HHL is one of the primary steps, but surely not the only step. Would something like recognizing genetic sequences using HHL as one of the primary steps (subject to all the constraints you mentioned), be a suitable answer? The other primary steps mainly involve Hamiltonian simulation and state preparation. – Sanchayan Dutta Jul 11 '18 at 15:25
• I would prefer some examples where HHL is directly used. It means that the problem can be directly formulated as a linear system of equation to solve. This is the case when solving differential equations: we discretise the equation and solve the discretised problem which is most of the time a sparse linear system. But other examples are welcomed. – Nelimee Jul 11 '18 at 15:29

A couple years ago it was shown in Quantum algorithms and the finite element method by Montanaro and Pallister that the HHL algorithm could be applied to the Finite Element Method (FEM) which is a "technique for efficiently finding numerical approximations to the solutions of boundary value problems (BVPs) for partial differential equations, based on discretizing the parameter space via a finite mesh".

They showed that within this context HHL could be used to achieve (perhaps at most) a polynomial speedup over the standard classical algorithm (the "conjugate gradient method").

With respect to real-world use-cases, they state that

"One example application is any dynamical problem involving $n$ bodies, which implies solving a PDE defined over a configuration space of dimension 2n. Also, there may be a significant advantage for problems in mathematical finance; for example, pricing multiasset options requires solving the Black-Scholes equation over a domain with dimension given by the number of assets"

This opens up a whole area of potential use-cases for HHL (assuming conditions on the sparsity of $A$ can be satisfied).

• I didn't downvote, but the question asks for specific instances where the matrix and vector satisfy the restrictions of the HHL algorithm. For example you say "this opens up a whole area of potential use-cases for HHL" but what is a useful example where A and b satisfy the HHL conditions? – user1271772 Jul 19 '18 at 18:42
• Where did it ask for a specific instance? As far as I can see they only asked for a real-world problem to solve, which I gave by the way of FEM and further an instance where FEM is applied, I.e Black-Sholes – SLesslyTall Jul 19 '18 at 18:48
• "Question: Taking into account all of these limitations..." means you have to take into account that A and b cannot just be any matrix and any vector. Otherwise the question is too easy. "Black-Sholes differential equation" is one application, so is any other differential equation. See the two deleted answers to this question and the reasons they deleted those answers (if you're able to). – user1271772 Jul 19 '18 at 18:52
• The whole point of the paper that I linked was that they ensured that both A and b used by FEM were sufficient to satisfy the HHL conditions to ensure a quantum speedup--not just any old system of simultaneous equations. The reason Black-Sholes is referenced is that it is is solvable using FEM as described by the paper. – SLesslyTall Jul 19 '18 at 18:57
• I don't have the time to read the article right now, but from the abstract it seems that the paper is interesting. I read quickly the section III and could not find any references to the eigenvalues of $M$ but the other points are either trivial ($M$ is hermitian), covered in other articles (hamiltonian simulation of $s$-sparse matrix with $s=3$) or covered in the article (state preparation). I did not know about this article, and the application is of particular interest for me (finite-elements is closely linked with PDEs). You have my upvote (and probably the bounty) :) – Nelimee Jul 20 '18 at 12:09

Quantum computers are devices that harness quantum mechanics to perform computations in ways that classical computers cannot. For certain problems, quantum algorithms supply exponential speedups over their classical counterparts, the most famous example being Shor's factoring algorithm. Few such exponential speedups are known, and those that are (such as the use of quantum computers to simulate other quantum systems) have so far found limited use outside the domain of quantum mechanics. This algorithm provides an exponentially faster method of estimating features of the solution of a set of linear equations, which is a problem ubiquitous in science and engineering, both on its own and as a subroutine in more complex problems.

Electromagnetic scattering

Clader et al. provided a preconditioned version of the linear systems algorithm that provided two advances. First, they demonstrated how a preconditioner could be included within the quantum algorithm. This expands the class of problems that can achieve the promised exponential speedup since the scaling of HHL and the best classical algorithms are both polynomials in the condition number. The second advance was the demonstration of how to use HHL to solve for the radar cross-section of a complex shape. This was one of the first end-to-end examples of how to use HHL to solve a concrete problem exponentially faster than the best known classical algorithm. [9]

Linear differential equation solving

Dominic Berry proposed a new algorithm for solving linear time-dependent differential equations as an extension of the quantum algorithm for solving linear systems of equations. Berry provides an efficient algorithm for solving the full-time evolution under sparse linear differential equations on a quantum computer.[10]

Least-squares fitting

Wiebe et al. provide a new quantum algorithm to determine the quality of a least-squares fit in which a continuous function is used to approximate a set of discrete points by extending the quantum algorithm for linear systems of equations. As the amount of discrete points increases, the time required to produce a least-squares fit using even a quantum computer running a quantum state tomography algorithm becomes very large. Wiebe et al. find that in many cases, their algorithm can efficiently find a concise approximation of the data points, eliminating the need for the higher-complexity tomography algorithm.[11]

Machine learning and big data analysis

Machine learning is the study of systems that can identify trends in data. Tasks in machine learning frequently involve manipulating and classifying a large volume of data in high-dimensional vector spaces. The runtime of classical machine learning algorithms is limited by a polynomial dependence on both the volume of data and the dimensions of the space. Quantum computers are capable of manipulating high-dimensional vectors using tensor product spaces are thus the perfect platform for machine learning algorithms.[12]

The quantum algorithm for linear systems of equations has been applied to a support vector machine, which is an optimized linear or non-linear binary classifier. A support vector machine can be used for supervised machine learning, in which training set of already classified data is available, or unsupervised machine learning, in which all data given to the system is unclassified. Rebentrost et al. show that a quantum support vector machine can be used for big data classification and achieve an exponential speedup over classical computers.[13]

Rebentrost et al. also recently used the HHL09 algorithm in their A Quantum Hopfield Neural Network (2018) paper, for optimization of the Hopfield network's energy function.

Basically, if the Lagrangian (which is used to optimize the network energy $E = -\frac{1}{2}x^{T}Wx + \theta^Tx$ given the constraint $Px - x^{\text{(inc)}} = 0$) is:

$$\mathcal{L} = -\frac{1}{2}x^{T}Wx + \theta^Tx - \lambda^T (Px - x^{\text{(inc)}}) + \frac{\gamma}{2}x^T x$$ then the optimization equations $\frac{\partial \mathcal{L}}{\partial x} = 0$ and $\frac{\partial \mathcal{L}}{\partial \lambda} = 0$ can be written in the form $A \mathbf{v} = \mathbf{w}$. Note that the $\gamma$ in the expression is the regularization parameter. We need to find $\mathbf{v}$ which extremizes network energy subject to the constraint $Px = x^{(\text{inc})}$ and thus, we need a matrix inversion technique. In the paper they've done exactly that and for the matrix inversion they utilized the HHL09 algorithm. See page 4 of the paper.

In short, I believe that once we have quantum computers with a sufficiently large number of qubits and decoherence time, the HHL algorithm is going to be one of the most useful subroutines for any quantum machine learning algorithm (since almost all machine learning and neural network algorithms involve some form of "gradient descent" or "optimization").

• But the question asks for specific examples given the restrictions on A and b, so even though you can solve matrix-valued differential equations, the matrices are restricted, and the question is what specific applications are there that satisfy these restrictions? So what useful differential equation satisfies the mentioned restrictions? – user1271772 Jul 11 '18 at 12:20
• Same goes for machine learning, what's a specific machine learning matrix and vector that satisfies the restrictions mentioned in the question? The quesiton is pointing out that HHL cannot be applied to all possible matrices and vectors, so are there any specific examples of problems that do satisfy these restrictions, that HHL can be applied to? – user1271772 Jul 11 '18 at 12:24
• @user1271772 In the Hopfield network paper I linked, the matrix is necessarily Hermitian since $w_{ij}=w_{ji}$. Also, the restriction on spectral norm has been considered in that paper (for the weighing matrix $||W||<1$) As for state preparation, that is still a problem, yes. Arbitrary states cannot be prepared efficiently. It has been addressed in the first part of the paper: "This could in principle be achieved using the developing techniques of quantum random access memory (qRAM) [33] or efficient quantum state preparation, for which restricted, oracle based, results exist [34]" – Sanchayan Dutta Jul 11 '18 at 12:33
• Also, see the related post: Is it possible to speed up the generation of the weighting matrix using a quantum algorithm? – Sanchayan Dutta Jul 11 '18 at 12:35
• The Hopfield example was not my problem. You copied and pasted a bunch of examples from Wikipedia (such as differential equation solving), which anyone could have done including the asker, but the reason I didn't (and presumably the reason why the asker asked this question in the first place) is because the asker seems to want (for example) specific differential equations where the matrices and vectors satisfy the restrictions imposed by algorithms of the HHL family. – user1271772 Jul 11 '18 at 12:44


text:
# Is it possible to find indefinite integral of $\int \frac{1} {{\sin(x)+\sec^2(x)}}\mathrm{d}x$?

I am in standard $XI$ (i.e.11) and newbie in learning topic of integration. My friend asked me to find indefinite integral of the example shown below

$$y=\int \frac{1} {{\sin(x)+\sec^2(x)}} \, \mathrm{d}x \tag 1$$

What I tried until now is I done substitution as $m=\sin(x)$

$$\frac{\textrm{d}m}{\textrm{d}x}=\cos(x)$$

Now, converting equation $(1)$ in terms of $m$ becomes as

$$y=\int \frac{(1-m^2)^{1/2}} {{1+m(1-m^2)}} \, \mathrm{d}m$$

but as you can see that it became more complicated than original equation $(1)$

So can anybody help me to compute integration of $\int \frac{1} {\sin(x)+\sec^2(x)} \, \mathrm{d}x$ ?

## 1 Answer

With the usual tangent half-angle substitution, $x=2\arctan t$, the integral becomes: $$I = 2\int \frac{(1-t^2)^2}{1+2t+3t^2-4t^3+3t^4+2t^5+t^6}\,dt\tag{1}$$ So, assuming we know the roots of the polynomial $p(t)=1+2t+3t^2-4t^3+3t^4+2t^5+t^6$, we can solve the above integral through partial fraction decomposition. That polynomial is palyndromic, so if $\zeta$ is a root, $\frac{1}{\zeta}$ is a root, too, and the original problem boils down to finding the roots of a third-degree polynomial.

For instance, by replacing $t+\frac{1}{t}$ with $u$, then $u$ with $2v$, we get: $$I = 2\int\frac{\sqrt{u^2-4}}{u^3+2u^2-8}\,du = \int \frac{\sqrt{v^2-1}}{v^3+v^2-1}\,dv\tag{2}$$ and by replacing $v$ with $\cosh z$ we have: $$I = \int \frac{\sinh^2 z}{\cosh^3 z+\cosh^2 z-1}\,dz. \tag{3}$$ Anyway, since the discriminant of $v^3+v^2-1$ is $-23$, the closed form of $(1)$ is not nice at all.

• Should it not be $\sinh(z)$ and not $\sinh^2(z)$? Upvote anyway! – johannesvalks Jul 11 '15 at 1:45
• @jack what in the answer written is high level mathematics. I am beginner in learning topic of integration and I would need explanation after each step. I am not understanding the given answer after equation (1) from your answer. can you write the answer in detail or simplify a little bit more ? – devraj Jul 11 '15 at 5:40
• A minor shortcut would come from not using the Weierstrass substitution, but rather rewriting the integrand in $\sin$ and $\cos$ as $\frac{\cos^2(x)}{\sin(x)\cos^2(x)+1}=\frac{\cos(x)\cdot\cos(x)}{\sin(x)(1-{\sin^2} (x))+1}$ and substituting $u=\sin(x)$. Still leads to the same kind of thing as in (2) of course. – alex.jordan Jul 11 '15 at 7:01
• @devraj: the substitutions are straightforward; the point is that our integral equals the integral of a rational function $\frac{p(x)}{q(x)}$ with $q(x)$ being a sixth-degree palyndromic polynomial, so $I$ cannot be simplified much more than in $(3)$. – Jack D'Aurizio Jul 11 '15 at 11:47


text:
In this update we're sharing the latest news on the development of the Sigmadex Protocol and a recap of this months progress. We would like to thank our community for their continued enthusiasm and are thrilled to be joined by such a supportive group of people. Together, we are making this powerful game-theory infused protocol come to life!

August has been an exciting month as our team is getting closer to achieving a number of major goals we had set out to achieve early on. During the conception phase of Sigmadex our team discussed the long term future and what sort of requirements would be necessary to ensure our longevity. While we cannot disclose exact details at this time, we can share that Sigmadex is pursuing multiple reputable launchpads for our TGE and have participated in talks with the necessary parties to bring this goal to fruition. While we have yet to reach this monumental goal, we are well on our way and will be sharing guaranteed information once it becomes available.

We are entertaining multiple reputable centralized exchange partners for the Sigmadex TGE for two main reasons. The first being for adequate token liquidity and volume and the second reason is low barrier to entry for newcomers to obtain SDEX tokens.

While some of our supporters have been eager to see us launch quickly, we cannot hastily launch Sigmadex simply for the sake of launching. We are undertaking the necessary strategic plans to align ourselves prior to launch. There are numerous examples of projects that launched their projects only to find their token values plummet in a short period of time shortly thereafter. This is the exact scenario we want to avoid. Our team is working feverishly to ensure the best possible launch for Sigmadex and truly appreciate our community and their understanding during this critical phase.

### Moving Closer to TGE

Our community is eager to know the official date of the Sigmadex TGE. Working to combine a centralized listing on a tier 1 launchpad for our TGE is the main bottleneck in establishing a hard date for this event. Working together with numerous partners to establish a date for our TGE means relying on external factors that are out of our hands. While our team can control the development of the Sigmadex protocol we still need to collaborate with other external parties to bring Sigmadex to market and establish an exact date for our TGE. We are working closely with our partners and are in constant communication to determine an exact date for this event and will be sharing this information as soon as it can be made public.

## Refining Penalty Formula

If you are familiar with Sigmadex you know that the reward/penalty system is an important mechanism for the protocol. This system actively drives the ecosystem; penalties are a major part of the cyclical momentum of the Sigmadex protocol since they are distributed to users who go full term on their liquidity stake. Our development team has been hard at work creating and testing the mechanics of this penalty formula. This formula is proposed to use Net Present Value Analysis with penalties being calculated in perpetuity and compared to the remaining time in the initial stake. This method allows a more efficient and accurate penalty calculation for Sigmadex users. With so much importance on this particular mechanism in-depth testing is required to ensure its success.

The penalty ratio formula for rebalancing the market after unstaking liquidity as a function of time elapsed, total time originally staked, and time horizon can be expressed as:

$$p(t)= \frac{T-t}{T+T_{long} -t}$$

where,
p(t) : penalty ratio for original staking amount
t : percent of time elapsed since staking began
T : total amount of time originally staked
$T_{long}$ : horizon period for rebalancing the market pool

Selecting the time horizon is a matter of analyzing the discount rate (risk-weighted interest rate available to the liquidity provider in the DeFi community) and the volatility horizon used in the Dynamic Automatic Market Maker algorithm. These relations are currently being studied to optimize the time horizon used for calculating the penalty for unstaking liquidity. As our development team works to refine this formula, changes may still be made in the near future.

## NFT Farming Event Audits

In recent months we announced a Public NFT Farming Event that would include a limited number of powerful NFTs for use with the Sigmadex protocol. These powerful tokens have the ability to override the standard parameters of the penalty/reward system and can provide holders with incredibly valuable one-off benefits.

To make sure these NFTs work as intended we have taken the steps over the last month to carefully begin auditing the function and testing of these various tokens. Since these NFTs work largely outside of the protocols dictated parameters extensive auditing and testing is required.

## Market Conditions

Market conditions are improving since the brief downturn in early May with signs of stable and steady growth so far. The impact of the downward market allowed the Sigmadex team to focus on development, establishing new partnerships and continuing toward cementing factors regarding our TGE such as potential launchpads and liquidity partners. Since postponing our original timeframe for our TGE we have sparked a number of new relationships within the space, expanded our internal team and have on-boarded a number of prominent and experienced advisors to aid us in developing Sigmadex to it's fullest potential and bringing it to market. Paving the right path to ensure the best possible success for Sigmadex is critical.

While the market has stabilized our team remains adament that Sigmadex is immune to unfavorable market conditions. The market downturn did not cause any negative effects to our project rather it gave us the opportunity to grow substantially in every facet. We are building Sigmadex as a disrupter to the industry and we want our community to know that we will do everything within our power to ensure the longevity of this project. Crypto is constantly evolving and reinventing itself - This leaves opportunities for Sigmadex to continually improve, collaborate and add necessary features and functionality as time goes by.


text:
Differential and Integral Equations

Finite energy weak solutions to some Dirichlet problems with very singular drift

Lucio Boccardo

Abstract

In this paper, the boundary problems (1.1) and (3.1) are studied. The main results are the existence of a bounded weak solution of (1.1) under the minimal assumption (1.3) on $E$, and of the quasilinear problem (Hamilton-Jacobi equation) (3.1).

Article information

Source
Differential Integral Equations, Volume 32, Number 7/8 (2019), 409-422.

Dates
First available in Project Euclid: 2 May 2019

Permanent link to this document
https://projecteuclid.org/euclid.die/1556762423

Mathematical Reviews number (MathSciNet)
MR3945762

Citation

Boccardo, Lucio. Finite energy weak solutions to some Dirichlet problems with very singular drift. Differential Integral Equations 32 (2019), no. 7/8, 409--422. https://projecteuclid.org/euclid.die/1556762423


text:
# All Questions

3 views

20 views

### the average speed of the car

A car started its journey from town A to town B with a constant acceleration. From the time the speed of 80 km/h was reached, for half of the total travel time the car continued at this speed. The ...
8 views

### Is a map from a smooth curve, bijective away from the singularities, already a normalization?

Over $\mathbb{C}$, I am considering a map between projective curves $f: C \rightarrow C'$, where $C$ is smooth. Suppose that $f$ is surjective as well as bijective away from the singularities of ...
9 views

### Does “pseudo-independent implies independent” imply that $R$ is a field?

(All my rings are unital.) Suppose $R$ is a commutative ring and that $M$ is an $R$-module. Definition. Call a subset $X \subseteq M$ pseudo-independent iff for all proper subsets $Y$ of $X,$ the ...
8 views

### Proving max-mod principle by contradiction

This is a homework exercise I have to make which I am kind of stuck on. First let $U$ be open and connected, $\overline{D}$ be the closure of the disk $D$ contained in $U$ and ...
17 views

### What is the maximal number of residents of the city

The population of the city Roachrach consists of cockroaches and cucarachas. The number of residents of the city is no greater than 2,000,000. Each cockroach is acquainted with 1,000 cucarachas, and ...
12 views

### Perfect pairing induces isomorphism of tensor products

Let $M, N$ be $R$-modules and $(\cdot, \cdot): M \times N \to R$ be a perfect pairing. Wikipedia sais that this means that the map $\varphi: M \to \text{Hom}_R(N, R), m \mapsto (n \mapsto (m, n))$ is ...
36 views

### Proving the equation has no root.

How to show that for $a\in \mathbb R$, the equation $x^2+12a^2+4ax-8a+8=0$ has no root?
34 views

### $\frac{1}{{1 + \left| {\left\| A \right\|} \right|}} \le \left| {\left\| {{{(I - A)}^{ - 1}}} \right\|} \right|$

Let a matrix norm $\left| {\left\| . \right\|} \right|$ have the property that $\left| {\left\| I \right\|} \right| = 1$ and $\left| {\left\| A \right\|} \right| < 1$. Why does the following ...
26 views

### How to find one matrix, which is subject to $B^3 = A$. How much is such matrices?

Here I have a problem with row echelon form. $$A := \begin{bmatrix}-6 & 3 & 7 \\ 0 & -1 & 0 \\ -14 & 12 & 15\end{bmatrix}$$
16 views

### Cyclic projective module

Let $R$ be an integral domain. If $M$ is a cyclic $R$-module which is also projective, then must there necessarily be an isomorphism of $R$-modules $M \cong R$?
8 views

### Complexity of Newton iteration problem for a d-dimensional problem

If we assume that we have $f:\mathbb{R}^{d} \rightarrow \mathbb{R}^{d}$ and we want to use the Newton iteration method to solve $f(x)=0_{\mathbb{R}^{d} }$. Is there any theorem regarding the ...
26 views

### Is there an established notation for this “replacement” operation?

If $S$ is a set, define $$(x \to y) \cdot S := \begin{cases} (S \setminus \{x\}) \cup \{y\} & \text{ if } x \in S \text{ and } y \not \in S; \\ S & \text{ otherwise.} \end{cases}$$ In other ...
12 views

### Combinations for pairing groups

I have a little bit of a complex question and I don't know anything about combinatorics, but I'm working on software problem and I'm trying to figure out how my algorithm will scale. I'm having to ask ...
34 views

### How to prove that one guy in all groups [on hold]

I dont know how even think about it. Anyone? thanks
48 views

### $3^x + 4^y = 5^z$

This is an advanced high-school problem. Find all natural $x,y$, and $z$ such that $3^x + 4^y = 5^z$. The only obvious solution I can see is $x=y=z=2$. Are there any other solutions?
24 views

### What exactly is wrong with this argument (Lucas-Penrose fallacy)

Argument "For every computer system, there is a sentence which is undecidable for the computer, but the human sees that it is true, therefore proving the sentence via some non-algorithmic method." ...
11 views

### Prove that $\sum_{n=0}^\infty e^{-nz}$ is analytic in the right half plane $\text{Re}(z)>0$

Consider$$\sum_{n=0}^\infty e^{-nz}$$ Using Weierstrass theorem, prove that the series is analytic in $\text{Re}(z)>0$. I know that $f$ is analytic if it satisfies Cauchy–Riemann equations. Could ...
4 views

### How to deal with a barrier function when constrained variables reach their bounds?

I am implementing an algorithm of Dang and Xu's, Non-convex Quadratic Programming Problem with Box Constraints'' and I'm hoping that somebody could verify what I'm doing. Their algorithm minimizes ...
6 views

### anlyalytic paths through convergent cauchy sequence II

Assume we have a Cauchy sequence $\{\vec{a_i}:i\in\mathbb{N}\}$ converging to $\vec{0}$ in $\mathbb{C}^n$ such that $|\vec{a_i}|<|\vec{a_j}|$ whenever $i>j$. Can we find an analytic path ...
12 views

### Invertibility for a matrix that I don't know

I would like to know why $(e^{-At}-I)^{-1}$ is invertible when matrix A is Hurwitz.
3 views

### Direct Product of Chernikov Groups is Chernigov group?

A group $G$ is said to be Chernikov if it contains a normal subgroup N such that $G/N$ is finite and $N$ is direct product of finitely many Prufer groups. The problem is the following: If $G$ is a ...
45 views

### A logic problem about set theory

In a group of n people, subgroups with common interest are formed (football,tennis,snooker). The number of subgroups equals $2^{n-1}$. Any 3 subgroups have a common member. Prove that there is a ...
29 views

### Groups with finite automorphism groups.

An easy argument shows that for any finite group $G$ the cardinal of $Aut(G)$ is less than $(|G|-1)!$. In particular the automorphisms group of a finite group is finite. Basically my question is about ...
24 views

### minimal volume of pyramid x0,y0,z0

I dont know how should i even start. I tried to think about something but get nothing. can someone help me please? Thanks
30 views

### Why $p\{N>n\}=p\{X_1+…+X_n\leq x\}$.

Let $(X_k)$ a sequence iid of random variable uniform on $[0,1]$. Let $x\in]0,1[$ and $N=\min\{n\geq 1\mid X_1+...+X_n>x\}$. Why $$p\{N>n\}=p\{X_1+...+X_n\leq x\} \ \ ?$$
35 views

### a linear differential equation with periodic coefficients

Let $$y' = a(x) y + b(x)$$ be a linear differential equation with continuous, periodic coefficients $a, b: \mathbb{R} \to \mathbb{R}$ that both have a period of $T > 0$. Also, we assume that ...
32 views

### tetrahedron problem with center and reflections [on hold]

I can i solve it? help me please i dont know how to do it. thanks
Prove that for all prime numbers $p>2000$ the sum $1 + 2\cdot2000 + 3\cdot2000 + … + (p-1)\cdot2000$ is divisible by $p$.


text:
# Least Squares and Fourier Series

I have a little bit of problem figuring out the relation between Fourier series and Least Squares.

As far as I understand, LS is a way of minimizing the quadratic error between a measured value $y_i$ and a linear combination of functions evaluated at points $x_i$. So, in order to solve for least squares we have to have the $(x_i,y_i)$ and the functions $f_j(\cdot)$ that are going to be linearly combined.

If I would like to approximate a set of points $(x_i,y_i)$ using Fourier basis functions, I could assemble a matrix as

$$\begin{pmatrix} 1 & \cos(\omega x_1) & \sin(\omega x_1) \\ 1 & \cos(\omega x_2) & \sin(\omega x_2) \\ \vdots & \vdots & \vdots\\ 1 & \cos(\omega x_n) & \sin(\omega x_n) \\ \end{pmatrix} \begin{pmatrix} A_0 \\ A_1 \\ A_2 \\ \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{pmatrix}$$

To obtain the values of the coefficients $A_0,A_1,A_2$ using least-squares I can solve for $x = (A^TA)^{-1}A^Tb$

Now, my doubt arises since it seems to me that the coefficients $A_0,A_1,A_2$ computed with the least-squares formula can be found with the Fourier transform. Are the coefficients $A_0,A_1,A_2$ the same as the Fourier coefficients $C_k$ ?

EDIT

clc; close all;

N = 10;
x = 1:N;
%y = rand(1,N);
y = [0.7099 0.6745 0.7374 0.1144 0.7732 0.0021 0.0357 0.4804 0.8107 0.0410];

% Polynomial Fit
A = [(x.*x)' x' ones(N,1)];
c = A\y'

% Fourier Fit - K+1 Coefficients
A = [];
K = 8;
for k = 0:K
A = [A exp(-1j*2*pi*k*x/N)'];
end
c = A\y'
c = (fft(y)/length(y))'
% Test the approximation on a bigger range
x_nodes = linspace(0,10,1000);
f = c(1);
for i = 2:K
f = f + c(i)*exp(1j*2*pi*i*x_nodes/N);
end

If I try to use the complex form of the Fourier Series $f = C_0 + \sum_k C_ke^{-j2\pi kx/L}$ and follow the least-squares approach and compare the results to the coefficients of the FFT, the coefficients seem to be different. What am I missing?

• It's probably better to make a new question than add an edit with another question, but ... For real valued functions: $C_0=A_0$ and $C_k=(A_k-iB_k)/2$ for $k\ge 1$. Also, $C_{-k}=Conj(C_k)$. In the case of an even number of points the extra (highest wave number) coefficient is purely real and given by $C_{k_{\max}}=A_{k_{\max}}$. – Doug Lipinski May 20 '14 at 16:30
• Also, for Fourier methods, your interval length should be dx*(N+1) due to the implied periodicity. – Doug Lipinski May 20 '14 at 16:35

In a word, yes.

A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit of a 4 term Fourier series to a set of 20 data points.

On the other hand, the discrete Fourier transform of a set of points always gives the same number of Fourier coefficients as input points. Thus, if you start with 20 points you will get 20 Fourier coefficients. Additionally, the first N Fourier coefficients are exactly the same as a least squares fit of a Fourier series with only N terms.

The key here is that the Fourier basis is an orthogonal basis on a given interval. The math works out so that the least squares best fit based of a lower order Fourier series is exactly equivalent to the truncated FFT. See this PDF for more details.

To verify, you can easily try this yourself. Pick a set of data points and compute the DFT and the least squares best fit and compare them. Here's an example in MATLAB:

    clear

%interval length and number of points
L = 10;
N = 15;

%data points
x = (linspace(0,L,N+1))'; x(end)=[];
y = x.^2; %function values

% FFT:
yk = fft(y)/length(x); %divide by length(x) for normalization

% pick out the first half of the coefficients
% (coefficients are symmetric for real valued functions)
yk = yk(1:ceil((length(x)+1)/2));

% Transform coefficients to the form:
%       y = A0 + sum( Ak*cos(...) + Bk*sin(...) )
A0 = yk(1); %constant term
A = 2*real(yk(2:end)); %cosine coefficients
B = -2*imag(yk(2:end)); %sine coefficients

%account for difference in highest wave number for odd/even numbers of points
if mod(numel(x),2)==0
A(end) = A(end)/2;
end

%reconstruct the original function for verification
yy=A0;
for k=1:length(A)
yy = yy + A(k)*cos(2*pi/L*k*x) + B(k)*sin(2*pi/L*k*x);
end

%least squares fit, unknowns C = [A0 A(1) B(1) A(2) B(2) ... ]'

%number of sin() and cos() coefficients to use:
N = 4;

%set up matrix:
M = zeros(length(x),1+2*N);
M(:,1) = 1;
for k=1:N
M(:,2*k) = cos(2*pi/L*k*x);
M(:,2*k+1) = sin(2*pi/L*k*x);
end

C = M\y; %equivalent to (M'*M)\(M'*y)

% least squares coefficients:
A0L2 = C(1);
AL2 = C(2:2:end);
BL2 = C(3:2:end);

%reconstruct the function values
yyL2=A0L2;
for k=1:length(AL2)
yyL2 = yyL2 + AL2(k)*cos(2*pi/L*k*x) + BL2(k)*sin(2*pi/L*k*x);
end

% Plot:

figure(1)
clf
hold on
plot(x,y,'ko')
plot(x,yy,'b-')
plot(x,yyL2,'r--')
legend('data','Fourier Series','Least Squares','location','northwest')
title('Comparison of function values')
xlabel('x')
ylabel('y')

figure(2)
clf
hold on
plot(A,'bo')
plot(B,'bs')
plot(AL2,'r.')
plot(BL2,'rx')
plot(0,A0,'bo')
plot(0,A0L2,'r.')
legend('Fourier Series, A_k','Fourier Series, B_k','Least Squares, A_k','Least Squares, B_k','location','north')
title('Comparison of coefficient''s values')
xlabel('wave number, k')
ylabel('coefficients')

This code produces the following plots. The second plot in particular shows that the coefficients of a truncated Fourier series are the same as the least squares best fit.

• In the second plot of the LS and FS coefficients, what we see is the complex conjugate of the coefficients? – BRabbit27 May 20 '14 at 9:52
• No, these are the coefficients for the series of the form $f = A_0+\sum_k\left( A_k \cos\left(\dfrac{2\pi}{L}kx\right) + B_k \sin\left(\dfrac{2\pi}{L}kx\right) \right)$ so they are real valued. The x-axis is wave number and there are two coefficients ($A_k$ and $B_k$) for each wave number for a given method. I'll update the plot soon to make that more clear. – Doug Lipinski May 20 '14 at 14:05
• Oh, I see, I understand now. I tried to get the coefficients using the complex form of the Fourier Series $f = A_0+\sum_k C_ke^{j2\pi kx/L}$ but the coefficients doesn't seem to be the same than those from the FFT. What am I doing wrong? – BRabbit27 May 20 '14 at 14:52
• Take a look at the code I posted and at @JLC's answer. There is a relatively simple relationship between the coefficients with $A_0=C_0$, $A_k=2Re(C_k)$ for $k\ge 1$, and $B_k=-2Imag(C_k)$ for $k\ge 1$. Also, if the number of points is even then the coefficients for the highest wave number are slightly different with $A_{k_{\max}}=Re(C_{k_{\max}})$ and $B_{k_{\max}}=0$. – Doug Lipinski May 20 '14 at 15:39
• New link of the PDF file: pdf – user153245 Oct 23 '17 at 7:34

They're related, but not the same, for a square system, in which the least squares fit becomes an interpolation. You have the relationship that

$c_k = \pm (a_k+ib_k)/2$ if $k\neq 0$

and $c_0 = a_0$ for $k=0$. You can do this by noting that the Fourier expansion over a period of $\tau$ gives

$f(t) = a_0 + \sum_{k=1}^\infty a_k \cos(\frac{2\pi k}{\tau}t) + \sum_{k=1}^\infty b_k \sin(\frac{2\pi k}{\tau}t)$.

Rewriting using Euler's formula gives

$f(t) = \sum_{k=1}^\infty c_k e^{i(2\pi k/\tau)t}$

with the above relationship between $c$ and $a,b$ above.

I'm not sure if the analogy holds for least squares.

• Sorry, I didn't get the idea of your answer. Could you elaborate a bit more? The goal of this is to understand this paper where the authors "project" the non-integrable gradient field on to a set of Fourier basis function so as to have an integrable gradient field. There's a part where they propose the minimization function and then use the Fourier basis functions to achieve that but I don't really get the "projection" step. – BRabbit27 May 19 '14 at 19:12
• Sorry for the confusion - I was just giving the connection between the $A_0,A_1,A_2$ coefficients and the $C_k$ coefficients from the Fourier transform. There's a relation between these coefficients for the case when you have an equal number of samples and frequencies. However, I don't believe this is the case when you have more frequencies than samples (which involves solving a least squares problem instead of just a linear system). – Jesse Chan May 19 '14 at 19:32
• What do you mean by "equal number of samples and frequencies"? Could you give me an example were I have more frequencies than samples? – BRabbit27 May 19 '14 at 19:50
• Well, in your problem you appear to have more samples than freqencies (more frequencies than samples is underdetermined). Doug's answer addresses the case when you have more samples than frequencies, and shows that it nails the first $N$ coefficients. I had forgotten that this comes out from the Fourier matrix having orthogonal columns! – Jesse Chan May 20 '14 at 3:20