text:
## Review for Second part (section 007)

Note from the instructor: This review exam is not very representative of the second part of the course. It touches but a small portion of the material taught in this important part. I don’t feel this is a good study guide.

Things I miss in this exam, and will surely be part of the final: (1) Product rule. (2) Quotient rule. (3) Marginal Cost, Revenue and Profit. (4) Concavity and inflection points. (5) First and second derivative tests. (6) Interpretation of the derivative (other than slopes). (7) Optimization story problems.

1. Find the inflection points of the following functions
1. $y=2x^3-3x^2-36x+5$
2. $x^4+x^3-3x^2+2$
2. Compute the following slopes:
1. $f(x)=6x^2-4x+36$ at $(3,6)$
2. $f(x)=x^3-8x^2+2x+40$ at $(2,7)$
3. Find the intervals of increase and decrease for the following functions:
1. $f(x)=x^3-9x^2-48x+52$
2. $f(x)=x^2-11x+3$
4. Compute the derivative of the following functions
1. $6e^{5x}+e^{-x^2}$
2. $\sqrt{1-x^2}$
5. Find the absolute maximum and minimum of the function $f(x)=(x-3)(x+2)$ on the interval $[-5,2]$


text:
# Cubic Inch To Imperial Pint Calculator

Cubic Inch (in3):

Imperial Pint (uk) :

The Cubic Inch to Imperial Pint(uk) Calculator an online tool which shows Cubic Inch to Imperial Pint(uk) for the given input. Byju's Cubic Inch to Imperial Pint(uk) Calculator is a tool
which makes calculations very simple and interesting. If an input is given then it can easily show the result for the given number.

Solve for x:
2x+5=3(x+2)

format: plain


text:
# Nodes and weight of Gauss Quadrature

• MHB
Gold Member
MHB
Hey! :giggle:

Calculate the node $x_0$ and the weight $a_0$ of Gauss Quadrature so that $$\int_0^1w(x)f(x)\, dx\approx I_0(f)=a_0f(x_0)$$ where $w(x)=1+\sqrt{x}$.

I have done the following:

The Gauss quadrature formula with $(n + 1)=1$ node (i.e. $n=0$) integrates polynomials of degree $2n + 1=1$ exactly.
\begin{align*}\int_0^1w(x)\cdot 1\, dx=a_0 &\Rightarrow \int_0^1\left (1+\sqrt{x}\right )\cdot 1\, dx=a_0 \\ & \Rightarrow \int_0^1\left (1+\sqrt{x}\right )\, dx=a_0 \\ & \Rightarrow \int_0^1\left (1+x^{\frac{1}{2}}\right )\, dx=a_0 \\ & \Rightarrow \left [x+2x^{\frac{1}{2}+1}\right ]_0^1=a_0 \\ & \Rightarrow \left [x+2x^{\frac{3}{2}}\right ]_0^1=a_0 \\ & \Rightarrow 1+2 =a_0 \\ & \Rightarrow a_0 =3 \\ \int_0^1w(x)\cdot x\, dx=a_0\cdot x_0 &\Rightarrow \int_0^1\left (1+\sqrt{x}\right )\cdot x\, dx=3x_0 \\ & \Rightarrow \int_0^1\left (x+x\sqrt{x}\right )\, dx=3x_0 \\ & \Rightarrow \int_0^1\left (x+x^{\frac{3}{2}}\right )\, dx=3x_0 \\ & \Rightarrow \left [\frac{x^2}{2}+\frac{2}{3}x^{\frac{3}{2}+1}\right ]_0^1=3x_0 \\ & \Rightarrow \left [\frac{x^2}{2}+\frac{2}{3}x^{\frac{5}{2}}\right ]_0^1=3x_0 \\ & \Rightarrow \frac{1}{2}+\frac{2}{3} =3x_0 \\ & \Rightarrow \frac{7}{6} =3x_0 \\ & \Rightarrow x_0=\frac{7}{18} \end{align*}

Is that correct? Or if we say that it integrates polynomials of degree $1$ do we mean that $w(x)f(x)=1$ and $w(x)f(x)=x$ instead of $f(x)=1$ and $f(x)=x$ ?

:unsure:

Homework Helper
MHB
Or if we say that it integrates polynomials of degree $1$ do we mean that $w(x)f(x)=1$ and $w(x)f(x)=x$ instead of $f(x)=1$ and $f(x)=x$ ?

Hey mathmari!

I believe it is indeed intended that $f(x)=1$ and $f(x)=x$. (Nod)

$$\int_0^1\left (1+x^{\frac{1}{2}}\right )\, dx=a_0 \Rightarrow \left [x+2x^{\frac{1}{2}+1}\right ]_0^1=a_0$$

If we take the derivative of $2x^{\frac{1}{2}+1}$, then we don't get $x^{\frac{1}{2}}$ do we?

$$\int_0^1\left (x+x^{\frac{3}{2}}\right )\, dx=3x_0 \Rightarrow \left [\frac{x^2}{2}+\frac{2}{3}x^{\frac{3}{2}+1}\right ]_0^1$$

If we take the derivative of $\frac{2}{3}x^{\frac{3}{2}+1}$, then we don't get $x^{\frac{3}{2}}$ do we?

Gold Member
MHB
If we take the derivative of $2x^{\frac{1}{2}+1}$, then we don't get $x^{\frac{1}{2}}$ do we?

If we take the derivative of $\frac{2}{3}x^{\frac{3}{2}+1}$, then we don't get $x^{\frac{3}{2}}$ do we?

Oh yes...

It should be:
\begin{align*}\int_0^1w(x)\cdot 1\, dx=a_0 &\Rightarrow \int_0^1\left (1+\sqrt{x}\right )\cdot 1\, dx=a_0 \\ & \Rightarrow \int_0^1\left (1+\sqrt{x}\right )\, dx=a_0 \\ & \Rightarrow \int_0^1\left (1+x^{\frac{1}{2}}\right )\, dx=a_0 \\ & \Rightarrow \left [x+\frac{1}{\frac{1}{2}+1}x^{\frac{1}{2}+1}\right ]_0^1=a_0 \\ & \Rightarrow \left [x+\frac{2}{3}x^{\frac{3}{2}}\right ]_0^1=a_0 \\ & \Rightarrow 1+\frac{2}{3} =a_0 \\ & \Rightarrow a_0 =\frac{5}{3}\approx 1.6667 \\ \int_0^1w(x)\cdot x\, dx=a_0\cdot x_0 &\Rightarrow \int_0^1\left (1+\sqrt{x}\right )\cdot x\, dx=\frac{5}{3}x_0 \\ & \Rightarrow \int_0^1\left (x+x\sqrt{x}\right )\, dx=\frac{5}{3}x_0 \\ & \Rightarrow \int_0^1\left (x+x^{\frac{3}{2}}\right )\, dx=\frac{5}{3}x_0 \\ & \Rightarrow \left [\frac{x^2}{2}+\frac{1}{\frac{3}{2}+1}x^{\frac{3}{2}+1}\right ]_0^1=\frac{5}{3}x_0 \\ & \Rightarrow \left [\frac{x^2}{2}+\frac{2}{5}x^{\frac{5}{2}}\right ]_0^1=\frac{5}{3}x_0 \\ & \Rightarrow \frac{1}{2}+\frac{2}{5} =\frac{5}{3}x_0 \\ & \Rightarrow \frac{9}{10} =\frac{5}{3}x_0 \\ & \Rightarrow x_0=\frac{27}{50}=0,54 \end{align*}
:geek:

Homework Helper
MHB
It should be:
:geek:
It looks correct to me now. (Nod)

Gold Member
MHB
It looks correct to me now. (Nod)

Great! Thank you! 👏


text:
# Computational Biomechanics

Dr. Kewei Li

## Implementation of Material Model in FEAPpv

In FEAPpv as well as the full version of FEAP, you can write a Fortran subroutine to add a new constitutive model to the main program. Not only new constitutive models, you can also add a new element subroutine or even a new solver. To get started, you can study the built-in material models in FEAPpv. There are many examples of material models in FEAPpv. For example, the Neo-Hookean material model has been implemented in FEAPpv. You can download a copy of the FEAPpv source code online. The Neo-Hookean model as well as other material models are implemented in the subroutine elements\material.f.

The first step to implement your own material model in FEAPpv is to derive the Cauchy stress and elasticity tensors of your material model. Then, you can write a Fortran subroutine by using the material model interface available in FEAPpv. Once your user material model is ready, you can compile it together with the main FEAPpv program. The instruction to compile and build the FEAPpv executable file is described in the FEAPpv readme file. If you encountered an error saying that some .dll files such as libifcoremd.dll are missing while running the executable file for the first time on a Windows machine, please install the Intel Fortran Compiler Redistributable Libraries for your computer first.

As an example of the user material subroutine, the implementation of the Neo-Hookean hyperelastic material model in FEAPpv have been uploaded to the course repository together with a FEAPpv input file for the uniaxial extension test of a 2D rectangular model. Due to the limitation of FEAPpv, only the 2D mixed element is “available” (only partially, read below), it is not possible to test those user subroutines with 3D mixed elements in FEAPpv. Since we have already studied the analytical solution of a rectangular specimen under uniaxial test,

where $\lambda$ is the stretch in the loading direction. The following figure shows the initial (red frame) and deformed configurations of the rectangular model under uniaxial extension test up to a stretch of $\lambda = 1.5$.

Here is a comparison of the computational solution by using FEAPpv with user subroutines and the analytical solution of the uniaxial extension test for $\mu = 10$ kPa,

As shown, there is a small difference between the computational and analytical solutions. In fact, the total area of the deformed shape of the model is not equal to the total area of the model in the initial configuration. Thus, the material becomes compressible already. This is because the mixed 2D element in FEAPpv has not been completely implemented yet. The developers of FEAPpv are currently working to finish this implementation, and you may find it in a future release.

To make sure the user material subroutine is correct, a simple test with the 3D mixed elements in the full version of FEAP 8.5 (not FEAPpv because the 3D mixed element is not available in FEAPpv) shows exact match between the computational and analytical solutions.

The following figure shows the initial (red frame) and deformed configurations of the 3D model under uniaxial extension test up to a stretch of $\lambda = 1.5$ with FEAP 8.5.

Note that the thickness of the model becomes smaller in the deformed configuration too.

Go Back


text:
# Midpoint and Distance Formula on the Coordinate Plane

In the Cartesian Coordinate System, we take two different ponts and . These two points form a line segment with endpoints and .

We know that each line segment has the following two properties:

• Length (the distance between the endpoints)
• Midpoint (the point in the center of the line segment)

In this lesson, we will see how to calculate the length of the segment and how to find the coordinates of the midpoint, which is very easy by applying the rules that we will see shortly.

## Distance Formula

To measure the distance between two points in an orthogonal coordinate system, we use a formula that is very important in geometry. This formula is actually extracted from the Pythagorean Theorem as it is written below:

In which:

• is the length of the hypotenuse of a right triangle (the longest side).
• and the length of the two others sides (the shortest sides).

As a reminder, the Pythagorean Theorem only applies to a right-angled triangle

### Note:

Let and two points. The distance between and is denoted by , and we have:

But how did we get this formula? Let's first look at the figure below.

• The number is the length of the hypotenuse (Distance between and )
• The number is the length of the horizontal leg.
• The number is the length of the vertical leg.

So, applying the Pythagorean theorem, we have:

Remember that the distance is always positive.

## Midpoint Formula

Using geometry tools (for example, Compass), we can easily determine the midpoint of a line segment. But isn't there another way, without using geometry tools, that we can determine the midpoint?

Yes, of course, there is an arithmetic method that depends on knowing the coordinates of the two endpoints of the line segment in an orthogonal coordinate system, using a formula to calculate the coordinates of the center.

We take two points and in an orthogonal coordinate system, and take the point as the midpoint of the line segment .

### Note:

We say that is the midpoint of the line segment if the following two conditions are met together:

• There is the same distance between the two endpoints of the line segment, i.e.
• The point belongs to the line segment .

To calculate the midpoint coordinates, we use the following property:

### Note:

Let and two points and the midpoint of the line segment . We have:

The following figure shows the coordinates of the midpoint of a segment in a coordinate plane.

Let's take some examples to calculate distance and midpoint coordinates.

Example1:

Consider the points and .

1. Calculate the distance of .
2. Calculate the coordinates of the midpoint of the segment .
3. Deduce the lengths of the two segments and .

#### Solution:

1) Plotting the two points and in the coordinate plane.

Let and ; then we substitute in the distance formula:

2) We calculate the coordinates of using the midpoint formula:

3) is the midpoint of the line segment , This means that there is the same distance from both endpoints of the segment

Therefore,

Example2:

Consider the points and .

1. Calculate the distance of .
2. Calculate the coordinates of the midpoint of the segment .
3. Deduce the lengths of the two segments and .

#### Solution:

1) Pltting the two points and in the coordinate plane

.

Let and ; then we substitute in the distance formula:

2) We calculate the coordinates of using the midpoint formula:

3) is the midpoint of the line segment , This means that there is the same distance from both endpoints of the segment.

Therefore,

Example3:

Consider the points and .

1. Calculate the distance between and the origin.
2. Calculate the distance between and the origin.

#### Solution:

1) The origin of the coordinate plane is the point .

Let and ; then we substitute in the distance formula:

2) In the same way we calculate the distance between and the origin.

Let and ; then we substitute in the distance formula:

Example4:

Consider the points and .

1. Determine the coordinates of if is the midpoint of the line segment
2. Calculate the distance of the line segment

#### Solution:

1) Unlike examples 1 and 2, in this example we have the coordinates of the midpoint but we don't know the coordinates of the endpoint of the line segment . To solve this problem, it's very easy, we just apply the Midpoint Formula but this time differently.

We have and let

Since is the midpoint of then:

The coordinates on both sides of the equation are equal. It means that the -coordinate on the left side is equal to the -coordinate on the right side.

The mathematical expression is:

And the -coordinate on the left side is equal to the -coordinate on the right side.

The mathematical expression is:

Here we get two simple equations:

• The unknown in equation (1) is
• The unknown in equation (2) is

Let's solve equation (1), we have

Let's solve equation (2), we have:

Finally, the coordinates of the endpoint of the line segment are

2) Let's calculate the distance of the line segment

Let and

Therefore, the length of the line segment is:


text:
# When is it inappropriate to control for a variable?

I can think of at least one naive example. Suppose I want to study the relationship between X and Z. I also suspect that Y influences Z, so I control for Y. However, as it turns out, unbeknownst to me, X causes Y, and Y causes Z. Therefore, by controlling for Y, I "cover up" the relationship between X and Z, since X is independent of Z given Y.

Now, in the previous example, it may be the case that the relationships I should be studying are the ones between X and Y, and Y and Z. However, if I knew such things a priori, I wouldn't be doing science in the first place. The study that I DID do now suggests that there is no relationship between X and Z, which is not the case.... X and Z ARE related.

This is illustrated in the following dependence-diagram. In the right scenario, Z depends on X and Y and X and Y are independent. We rightly control for Y to determine the relationship between X and Z. In the left scenario Z depends on Y which depends on X. X and Z are independent given Y, so the relationship between X and Z is "covered up" by controlling for Y.

My question is basically "When is it appropriate to control for variable Y and when not?"... It may be difficult or impossible to fully investigate the relationship between X and Y, but, for instance, controlling Y at a given level is an option. How do we decide before conducting our study, and what are common pitfalls of controlling too much or too little?

Citations appreciated.

• To give an example, your exact situation comes up in assessing the impact of racial discrimination. Let $X$ be race. Let $Z$ be wages. Let $Y$ be education. Clearly education impacts wages, so you'd want to control for that, BUT if racial discrimination causes a racial minority to receive worse education then controlling for education would cover up that discrimination. Eg. see Neal and Johnson (1996). As Alexis's answer points out, you need to get into the particulars of your problem. There isn't a simple button to press that solves everything. – Matthew Gunn Jan 30 '18 at 17:38
• I had no hopes for a simple button to press. Indeed, I would have been quite disappointed if my question had turned out to have a trivial answer. :) – Scott Jan 30 '18 at 17:46
• @Repmat Yes. However, IV estimation's 4 assumptions are seldom met, and even when they are, the strength of the associations involved can make IV estimation give biased results. See, for example, Hernán and Robins Causal Estimation (full cite and link in my answer), chapter 16: Instrumental variables estimation. – Alexis Jan 30 '18 at 19:44
• @Alexis Naturally, IVs are hard to come by... "No free lunch" and all, however when you do, you will often know so with little or no doubt. – Repmat Jan 30 '18 at 20:19
• @Repmat... the assumptions for valid IV estimates require more than that pictured DAG... they really are quite frail. – Alexis Jan 30 '18 at 20:54

Conditioning (i.e. adjusting) the probabilities of some outcome given some predictor on third variables is widely practiced, but as you rightly point out, may actually introduce bias into the resulting estimate as a representation of causal effects. This can even happen with "classical" definitions of a potential causal confounder, because both the confounder itself, and the predictor of interest may each have further causal confounders upstream. In the DAG below, for example, $L$ is a classic confounder of the causal effect of $E$ on $D$, because (1) it causes and is therefore associated with $E$, and (2) is associated with $D$ since it is associated with $U_{2}$ which is associated with $D$. However, either conditioning or stratifying $P(D|E)$ on $L$ (a "collider") will produce biased causal estimates of the effect of $E$ on $D$ because $L$ is confounded with $D$ by the unmeasured variable $U_{2}$, and $L$ is confounded with $E$ by the unmeasured variable $U_{1}$.

Understanding which variables to condition or stratify one's analysis on to provide an unbiased causal estimate requires careful consideration of the possible DAGs using the criteria for causal effect identifiability—no common causes that are not blocked by backdoor paths—described by Pearl, Robins, and others. There are no shortcuts. Learn common confounding patterns. Learn common selection bias patterns. Practice.

References

Greenland, S., Pearl, J., and Robins, J. M. (1999). Causal diagrams for epidemiologic research. Epidemiology, 10(1):37–48.

Hernán, M. A. and Robins, J. M. (2018). Causal Inference. Chapman & Hall/CRC, Boca Raton, FL

Maldonado, G. and Greenland, S. (2002). Estimating causal effects. International Journal of Epidemiology, 31(2):422–438.

Pearl, J. (2000). Causality: Models, Reasoning, and Inference. Cambridge University Press.

When is it appropriate to control for variable Y and when not?

is the "back-door criterion".

Judea Pearl's Structural Causal Model can tell you definitively which variables are sufficient (and when it's necessary) for conditioning, to infer the causal impact of one variable on another. Namely, this is answered using the back-door criterion, which is described in page 19 of this review paper by Pearl.

The major caveat is that it requires you to know the causal relationship between the variables (in the form of directional arrows in a graph). There's no way around that. This is where the difficulty and possible subjectivity can come into play. Pearl's structural causal model only allows you to know how to answer the right questions given a causal model (i.e. directed graph), which set of causal models are possible given a data distribution, or how to look for causal structure by performing the right experiment. It doesn't tell you how to find the right causal structure given only the data distribution. In fact, it claims that this is impossible without using external knowledge/intuition about the meaning of the variables.

The back-door criteria can be stated as follows:

To find the causal impact of $X$ on $Y,$ a set of variable nodes $S$ is sufficient to be conditioned on as long as it satisfies both of the following criteria:

1) No elements in $S$ is a descendant of $X$

2) $S$ blocks all "back-door" paths between $X$ and $Y$

Here, a "back-door" path is simply a path of arrows that begin at $Y$ and end with an arrow pointing at $X.$ (The direction that all other arrows point is not important.) And "blocking" is, itself, a criterion that has a specific meaning, which is given in page 11 of the above link. This is the same criterion that you would read when learning about "D-separation". I personally found that chapter 8 of Bishop's Pattern Recognition and Machine Learning describes the concept of blocking in D-separation far better than the Pearl source I linked above. But it goes like this:

A set of nodes, $S,$ blocks a path between $X$ and $Y$ if it satisfies at least one of the following criteria:

1) One of the nodes in the path, that is also in $S,$ emits at least one arrow on the path (i.e. the arrow is pointing away from the node)

2) A node that is neither in $S$ nor an ancestor of a node in $S$ has two arrows in the path "colliding" towards it (i.e. meeting it head-to-head)

This is an or criterion, unlike the general back-door criterion which is an and criterion.

To be clear about the back-door criterion, what it tells you is that, for a given causal model, when conditioning on a sufficient variable, you can learn the causal impact from the probability distribution of the data. (As we know, the joint distribution alone isn't sufficient for finding causal behavior because multiple causal structures can be responsible for the same distribution. This is why the causal model is required as well.) The distribution can be estimated using ordinary statistical/machine learning methods on the observational data. So as long as you know that the causal structure allows for conditioning on a variable (or set of variables), your estimate of the causal impact of one variable on another is as good as your estimate of the distribution of the data, which you obtain through statistical methods.

Here is what we find when we apply the back-door criterion to your two diagrams:

In neither case does there exist a back-door path from $Z$ to $X.$ So it's true that $Y$ blocks "all" back-door paths, because there aren't any. However, in the left diagram, $Y$ is a direct descendant of $X,$ while in the right diagram it is not. Therefore $Y$ follows the back-door criterion in the right diagram, but not the left. These are unsurprising results.

What is surprising, however, is that in the right diagram, as long as it is the complete picture, you need not condition on $Y$ to get the full causal impact of $X$ on $Z$. (Said another way, the null set satisfies the back-door criteria, and is thus, sufficient for conditioning.) Intuitively this is true because the value of $X$ is not associated with that of $Y$ so for sufficient data you can simply average over the values of $Y$ to marginalize the effect of $Y$ on $Z.$ One objection to this point can be that the data is limited, so you don't have a representative distribution of $Y$ values. But recall that the back-door criterion assumes you have the probability distribution of the data. In that case you can analytically marginalize $Y.$ Marginalization over a finite data set is just an estimation. Also, note that it's highly unlikely that this is the complete picture. There are likely external factors that impact $X.$ If those factors are also associated with $Y$ in any way, then more work must be done to see if $Y$ must be conditioned on, or if it's even sufficient. If you draw another arrow pointing from $Y$ to $X$ then $Y$ becomes necessary to control.

Those are, of course, very simple examples where intuition is enough to know when $Y$ can or cannot be controlled for. But here's a couple more examples where it's not obvious by looking at the diagram, and you can use the back-door criteria. For the following diagram we ask if it is sufficient to control for $Y$ when determining the causal impact of $X$ on $Z.$

The first thing to note is that, in both cases, $Y$ is not a descendant of $X.$ So it passes that criterion. The next thing to note is that, in both cases, there are several backdoor paths from $Z$ to $X.$ Two in the left diagram and three in the right.

In the left diagram the backdoor paths are $Z \leftarrow Y \rightarrow X$ and $Z \leftarrow W \rightarrow B \leftarrow A \rightarrow X. \hspace{1mm}$ $Y$ blocks the first path because it is an arrow-emitting node that is directly in the path. $Y$ also blocks the second path because it is neither $B,$ nor is it a descendant of $B,$ which is the only arrow colliding node in the path. Therefore $Y$ is a sufficient set for conditioning. (Note, unlike in your right diagram, the null set is not sufficient for conditioning, because it does not block the path $Z \leftarrow Y \rightarrow X$.)

In the right diagram the backdoor paths are the same two as in the left, plus the path $Z \leftarrow W \rightarrow B \rightarrow Y \rightarrow X. \hspace{1mm}$ $Y$ does block this path, because it is an arrow emitting node in the path. It also blocks the path $Z \leftarrow Y \rightarrow X$ for the same reason as the left diagram. However, it does not block the path $Z \leftarrow W \rightarrow B \leftarrow A \rightarrow X,$ because it is a direct descendant of the collider node $B.$ Therefore it is not sufficient for conditioning.

It is pretty unintuitive to see why $Y$ is sufficient for conditioning on the left diagram, because of the exogenous variables $A$ and $W$ that affect $X$ and $Z$ respectively. However, suppose there was no $B.$ In that case, there would be no spurious relationship between $X$ and $Z$ due to these exogenous variables so they're not of concern. The existence of $B,$ however, puts that in question. If $B$ is allowed to take whatever value it naturally takes given $A$ and $W$, it would not be a problem because it does not have any impact on the important variables, or the exogenous variables determining them. However, if $B$ (or any of its descendants) is controlled then it actually renders $A$ and $W$ dependent, which creates the spurious relationship beteween $X$ and $Z$ that we don't want. As mentioned in the linked source, this is an example of Berkson's paradox, where an observation of a variable caused by two independent sources makes those sources dependent (e.g. the result of two independent coin flips is rendered dependent upon the observation of the number of total heads flipped).

As I mentioned before the use of the back-door criterion requires that you know the causal model (i.e. the "correct" diagram of arrows between the variables). But the Structural Causal Model, in my opinion, also gives the best and most formal way to search for such a model, or to know when the search is futile. It also has the wonderful side effect of rendering terms like "confounding", "mediation", and "spurious" (all of which confuse me) obsolete. Just show me the picture and I'll tell you which circles ought to be controlled.

• Nice. I was debating whether to add Pearl's Causality to the references section of my answer.... and now have done so. :) – Alexis Jan 30 '18 at 20:09

The following might or might not be appropriate to your case: if X is a treatment, then you might be able to go around your problem by using propensity score matching in which you would still keep the variable Y when you do the matching. In other words, you balance the covariates (Y is one of such covariates) that predict receiving the treatment X.
Note how there is no reference to the outcome variable Z in the above. You can also check how balanced your observations are (by generating a before and after matching balance table), which might give you insights into how much of X is determined by Y.


text:
Page updated: May 7, 2021
Author: Curtis Mobley
View PDF

# Counting Photons

[Bryan Monosmith, Jeremy Werdell, and Curtis Mobley contributed to this page.]

The foundation of ocean color remote sensing is sunlight that has entered the ocean, been transformed through absorption and scattering by the myriad constituents of the water body, and then been scattered out of the ocean and into a detector. It is worthwhile to consider how many photons going through this process are actually available for detection by a satellite sensor. An order-of-magnitude calculation suﬃces to identify many of the engineering constraints on the design of an ocean color sensor.

### Radiance from the Sea Surface

Figure 1 shows that the downwelling spectral plane irradiance ${E}_{d}$ onto the sea surface on a clear day is of order $1\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}$ at visible wavelengths. Most (90-98%, depending on sun zenith angle and wind speed; e.g. Light and Water Figs. 4.11-4.13) of this irradiance enters the ocean. In-water irradiance reﬂectances $R={E}_{u}∕{E}_{d}$ are typically in the 0.01 to 0.05 range, depending on the water constituents and wavelength ($R$ can reach 0.1 in very turbid, highly scattering waters).

Suppose that $R=0.03$ of downwelling plane irradiance of magnitude $1\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}$ entering the ocean is backscattered into upward directions. If this upwardly scattered light is isotropically scattered into the $2\pi \phantom{\rule{2.6108pt}{0ex}}sr$ of the upward hemisphere, then the upwelling isotropic radiance just below the sea surface would be

 ${L}_{u}=\frac{0.03}{2\pi \phantom{\rule{2.6108pt}{0ex}}sr}\left(1\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}\right)\phantom{\rule{2.6108pt}{0ex}}.$

For ocean waters the upwelling radiance is not isotropic. Just below the sea surface, upwelling radiance within 30 deg of the zenith (the angles relevant to remote sensing) is less than the radiance in more nearly horizontal directions by a factor of two or three, as shown in Fig. 2. Thus the above value for isotropic scattering can be reduced by a factor of roughly $\frac{1}{2}$ to estimate the underwater radiance in near-zenith directions.

A beam of this radiance just below the sea surface will be reduced by a factor of $t∕{n}^{2}$ when passing through the sea surface. Here $t$ is the water-to-air radiance transmittance, which is close to 1 for directions relevant to remote sensing, namely directions within a few tens of degrees of the zenith. $n\approx 1.34$ is the water index of refraction. The water-leaving radiance is then roughly

 ${L}_{w}=\frac{0.03}{2{n}^{2}2\pi }\approx 0.0013\phantom{\rule{1em}{0ex}}\left[W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}s{r}^{-1}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}\right]$

Atmospheric transmittance is 0.7 to 0.95 at visible wavelengths, so most of the water-leaving radiance will be transmitted to the top of the atmosphere (TOA), where it can be detected by a satellite. However, along the way, atmospheric scattering of solar radiation into the beam will add typically 10 to 20 times as much radiance to the beam. The TOA radiance seen by a satellite would then be of order $10×{L}_{w}\approx 0.01\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}s{r}^{-1}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}$. Typical open-ocean TOA radiances for the MODIS sensor are in the range of 0.01 (red wavelengths) to 0.08 (blue wavelengths) $W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}s{r}^{-1}\phantom{\rule{0.3em}{0ex}}n{m}^{-1}$.

### Photons Detected

Now that we have the radiance detected at the satellite, we can compute the numbers of photons collected and consider the related engineering matters. We start by computing the number of detected photons that come from a $1\phantom{\rule{2.6108pt}{0ex}}{m}^{2}$ patch of sea surface in 1 second of observation time. For speciﬁc numbers, the MODIS sensors have the following characteristics:

 Physical quantity value altitude 705 km = $7.1×1{0}^{5}\phantom{\rule{2.6108pt}{0ex}}m$ oﬀ-nadir viewing angle 20 deg slant range 750 km = $7.5×1{0}^{5}\phantom{\rule{2.6108pt}{0ex}}m$ sensor aperture radius 89 mm = 0.09 m band width $\Delta \lambda =10\phantom{\rule{2.6108pt}{0ex}}nm$ optical eﬃciency OE = 0.6 quantum eﬃciency QE = 0.9

Table 1: MODIS sensor characteristics.

The slant range is the distance from the satellite to the observed point on the ocean surface along the line of sight, which is taken here to be 20 deg oﬀ nadir. The optical eﬃciency is the fraction of light incident onto the sensor fore-optics that eventually reaches sensor material itself; the losses are due to reﬂections from lens surfaces or diﬀraction gratings, absorption by ﬁlters, etc.

The quantum eﬃciency is the fraction of photons reaching the sensor material that actually results in a signal, e.g., by the generation of a photo-electron. Figure 3 shows the quantum eﬃciencies of typical sensors. The ”standard CCD” curve is for CCDs like those used in consumer-grade video cameras. “Scientiﬁc CCD” is for a much more expensive “science grade” CCD. The MODIS curve shows what is achievable if your budget is almost unlimited. The bottom curves are for the human eye. The human eye has a maximum QE of only a few percent for color (photopic) vision and 10% for night-time (scotopic) vision. It thus seems that evolution has given us eyes that are adequate for ﬁnding food and avoiding tigers, and even for reading websites, but which are shockingly ineﬃcient from an optical engineering standpoint. (This ﬁgure gives another refutation of “intelligent design” for the human eye. The intelligent design in this ﬁgure was done by the physicists and optical engineers working on the CCD and MODIS sensors.)

The solid angle of the sensor as seen from the Earth’s surface is

 ${\Omega }_{aperture}=\frac{\pi {\left(aperture\phantom{\rule{2.6108pt}{0ex}}radius\right)}^{2}}{{\left(slant\phantom{\rule{2.6108pt}{0ex}}range\right)}^{2}}=4.5×1{0}^{-14}\phantom{\rule{2.6108pt}{0ex}}sr\phantom{\rule{2.6108pt}{0ex}}.$

The power detected by the sensor coming from a square meter of the ocean surface is

$\begin{array}{llll}\hfill {P}_{detector}=& L\phantom{\rule{2.6108pt}{0ex}}{\Omega }_{aperture}\phantom{\rule{2.6108pt}{0ex}}Are{a}_{surface}\phantom{\rule{2.6108pt}{0ex}}OE\phantom{\rule{2.6108pt}{0ex}}\Delta \lambda \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \left(0.01\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2.6108pt}{0ex}}{m}^{-2}\phantom{\rule{2.6108pt}{0ex}}s{r}^{-1}\phantom{\rule{2.6108pt}{0ex}}n{m}^{-1}\right)\phantom{\rule{2.6108pt}{0ex}}\left(4.5×1{0}^{-14}\phantom{\rule{2.6108pt}{0ex}}sr\right)\phantom{\rule{2.6108pt}{0ex}}\left(1\phantom{\rule{2.6108pt}{0ex}}{m}^{2}\right)\phantom{\rule{2.6108pt}{0ex}}\left(0.6\right)\phantom{\rule{2.6108pt}{0ex}}\left(10\phantom{\rule{2.6108pt}{0ex}}nm\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill \approx & \phantom{\rule{2.6108pt}{0ex}}2.7×1{0}^{-15}\phantom{\rule{2.6108pt}{0ex}}W\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ [This equation has been written for viewing the satellite from the sea surface. We could also take the viewpoint of the satellite viewing the earth. In that case, the relevant solid angle would be that of the $1\phantom{\rule{2.6108pt}{0ex}}{m}^{2}$ pixel as seen from the satellite, and the area factor would be that of the sensor aperture. The throughput, or etendue, of an optical system is
 ${\Omega }_{aperture}\phantom{\rule{2.6108pt}{0ex}}Are{a}_{surface}={\Omega }_{surface}\phantom{\rule{2.6108pt}{0ex}}Are{a}_{aperture}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}},$

which shows that these viewpoints are equivalent.]

If we assume a wavelength of 550 nm, the corresponding number of photo-electrons released in the detector in time $t=1\phantom{\rule{2.6108pt}{0ex}}s$ is

$\begin{array}{llll}\hfill {N}_{electrons}=& {P}_{detector}\phantom{\rule{2.6108pt}{0ex}}t\phantom{\rule{2.6108pt}{0ex}}QE\phantom{\rule{2.6108pt}{0ex}}\frac{\lambda }{hc}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \left(2.7×1{0}^{-15}J\phantom{\rule{2.6108pt}{0ex}}{s}^{-1}\right)\left(1\phantom{\rule{2.6108pt}{0ex}}s\right)\phantom{\rule{2.6108pt}{0ex}}\left(0.9\right)\phantom{\rule{2.6108pt}{0ex}}\frac{550×1{0}^{-9}\phantom{\rule{2.6108pt}{0ex}}m}{\left(6.63×1{0}^{-34}\phantom{\rule{2.6108pt}{0ex}}J\phantom{\rule{0.3em}{0ex}}s\right)\phantom{\rule{2.6108pt}{0ex}}\left(3×1{0}^{8}\phantom{\rule{2.6108pt}{0ex}}m\phantom{\rule{2.6108pt}{0ex}}{s}^{-1}\right)}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill \approx & 6800\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ where $h$ is Planck’s constant and $c$ is the speed of light.

These 6800 electrons are for the total TOA radiance in this green band. As previously noted, the TOA radiance is typically 90% atmospheric path radiance, in which case only 680 of these electrons correspond to the water-leaving radiance from a square meter of the sea surface.

Kepler’s third law of planetary motion and Newton’s law of gravity give the relation between a satellite’s orbital period $T$ and the radius $r$ of its orbit:

 ${T}^{2}=\frac{4{\pi }^{2}{r}^{3}}{GM},$

where $G=6.67×1{0}^{-11}\phantom{\rule{2.6108pt}{0ex}}{m}^{3}\phantom{\rule{2.6108pt}{0ex}}k{g}^{-1}\phantom{\rule{2.6108pt}{0ex}}{s}^{-2}$ is Newton’s gravitational constant and $M=5.97×1{0}^{24}\phantom{\rule{2.6108pt}{0ex}}kg$ is the mass of the earth. For a satellite at altitude 705 km above the earth (whose mean radius is 6731 km), this gives a period of 6385 s. This corresponds to a speed of $v=\left(2\pi r\right)∕T=7318\phantom{\rule{2.6108pt}{0ex}}m\phantom{\rule{2.6108pt}{0ex}}{s}^{-1}$ relative to the ground. The time required for the satellite to travel 1 m is then $1.37×1{0}^{-4}\phantom{\rule{2.6108pt}{0ex}}s$. For 1 m spatial resolution, this short exposure time reduces the number of detected photons by a factor of $1.37×1{0}^{-4}$ compared to the 1 s collection time computed above. Thus the number of water-leaving photons collected during the time the satellite passes over the $1\phantom{\rule{2.6108pt}{0ex}}{m}^{2}$ area is only 0.093. Alas, the actual situation is even worse because the sensor is not collecting light from just one pixel, but from perhaps 1000 pixels as the sensor either scans back and forth or rotates to observe a wide swath to either side of the satellite nadir point. Thus, in the $1{0}^{-4}\phantom{\rule{2.6108pt}{0ex}}s$ time required for the satellite to move forward by 1 m, the sensor must collect photons from 1000 pixels, reducing the number for each to roughly 0.0001 photoelectrons per pixel. Collecting only 0.0001 photoelectrons per pixel would not yield a very good image.

These physical and orbital constraints show one reason why orbiting satellites do not obtain meter-scale ocean color imagery: There simply are not enough photons leaving the ocean surface from a square meter of area to form an image. Many more photons must be collected, and there are several ways to do this:

• View a larger surface area, which both increases the number of photons leaving the surface and allows for longer integration times.
• View the surface area for a longer time, e.g., from a geostationary satellite that can stare at the same point for very long times (but a geostationary satellite has an altitude of 36,000 km, which makes the solid angle much smaller).
• Get closer to the surface, e.g. by using an airborne sensor ﬂying at a few kilometers above the sea surface. This greatly increases the solid angle of the sensor and allows for longer integration times for a slowly ﬂying aircraft.
• Increase the bandwidth.
• Increase the aperture of the receiving optics.
• Use multiple detector elements to observe the same ground pixel nearly simultaneously, either on the same or successive scans, and then combine the photons collected from the diﬀerent sensors.

Suppose, as is typical of ocean color sensors such as MODIS, that we image a $1\phantom{\rule{2.6108pt}{0ex}}k{m}^{2}$ area of ocean surface. This increases the number of photons leaving the imaged pixel by a factor of $1{0}^{6}$ and the integration time by a factor of $1{0}^{3}$ (the increased time for the satellite to travel 1 km rather than 1 m). The sensor now collects 93,000 photoelectrons from the water-leaving radiance. The total number of TOA photoelectrons, including the atmospheric path radiance, would be roughly ten times as large, about $1{0}^{6}$. In practice, this number will be less because of the duty-cycle time of the sensor: a rotating sensor will be viewing the ocean only about one-third of the time, and a scanning sensor will require time to stop and start each scan. Nevertheless, we can still collect roughly $1{0}^{6}$ photoelectrons for each pixel, of which $1{0}^{5}$ correspond to water-leaving radiance.

The signal to noise ratio, SNR, is given in general by

 $SNR=\frac{{N}_{PE}}{\sqrt{{\mathsc{𝒩}}_{PE}^{2}+{\mathsc{𝒩}}_{DC}^{2}+{\mathsc{𝒩}}_{RO}^{2}+{\mathsc{𝒩}}_{QN}^{2}}}\phantom{\rule{2.6108pt}{0ex}},$

where the signal ${N}_{PE}$ is the number of photoelectrons counted. The terms in the denominator represent the noise terms for the photoelectrons (PE), dark current (DC), sensor read-out (RO), and quantization (QN). Dark current noise results from the spontaneous emission of photoelectrons within the sensor, read-out noise comes from the sensor’s analog front-end electronics when the collected photoelectrons are read from the sensor, and quantization noise comes from uncertainty when the analog signal is digitized. (Writing the total noise as the square root of the sum of the individual noise terms squared assumes that the individual noise processes are uncorrelated.) The number of photoelectrons counted during a given time interval is described by a Poisson probability distribution, for which the noise is the standard deviation of the distribution, which in turn equals the number ${N}_{PE}$ of photoelectrons counted. The emission of dark current photoelectrons is also a Poisson process. The SNR can thus be written

 $SNR=\frac{{N}_{PE}}{\sqrt{{N}_{PE}+{N}_{DC}+{\mathsc{𝒩}}_{RO}^{2}+{\mathsc{𝒩}}_{QN}^{2}}}\phantom{\rule{2.6108pt}{0ex}}.$

Assuming no dark current or other noise, this gives an SNR of

 $SNR=\frac{1{0}^{6}}{\sqrt{1{0}^{6}}}=1000.$

for the number of photoelectrons estimated above. The actual SNR will be somewhat less due to dark current, read out, and quantization noise. The MODIS sensor has SNR requirements of 750 to 1000 for various bands, so we have achieved the approximate number of photons needed at the sensor.

These simple estimates illustrate the severe constraints on the design of any ocean color sensor. In practice, the engineering of a satellite ocean color sensor requires great sophistication to achieve the needed SNR.

#### Comments for Counting Photons:

Loading Conversation


text:
Looking for Algebra worksheets?
Check out our pre-made Algebra worksheets!
 Tweet

# Number Properties Questions - All Grades

You can create printable tests and worksheets from these Number Properties questions! Select one or more questions using the checkboxes above each question. Then click the add selected questions to a test button before moving to another page.

1 2 3 4 ... 9
Which property is illustrated by the equation $ax + ay = a(x + y)$?
1. commutative
2. associative
3. distributive
4. multiplicative inverse
Using the commutative property, what is $2xx3xx8$ equal to?
1. $2xx6$
2. $2-3-8$
3. $2+3+8$
4. $3xx2xx8$
Which multiplication property does the equation below represent?
$7*(9*2)= (7*9)*2$
1. Commutative Property
2. Identity Property
3. Associative Property
4. Zero Property
Choose the property described below:

$8+7=7+8$
1. Reflexive Property
2. Symmetric Property
3. Transitive Property
Which of these shows the commutative property?
1. 1+3=4 and 4+3=1
2. 1+2=3 and 2+1=3
Grade 1 Number Properties CCSS: 1.OA.B.3
Which is the related addition fact for $4 + 8 = 12$?
1. $6 + 6 = 12$
2. $7 + 5 = 12$
3. $8 + 4 = 12$
4. $9 + 3 = 12$
Grade 1 Number Properties CCSS: 1.OA.B.3
Which is the related addition fact for $5 + 4 = 9$?
1. $1 + 8 = 9$
2. $2 + 7 = 9$
3. $3 + 6 = 9$
4. $4 + 5 = 9$
Which of the following is an example of the Closure Property of Multiplication?
1. 72 + 51 = 51 + 72
2. 65 + 0 = 65
3. 6 x 4 = 24
4. (2 + 4) + 5 = 2 + (4 + 5)
Grade 6 Number Properties CCSS: 6.EE.A.3
Using the distributive property, what is $7(x + 9)$ equal to?
1. $(7x) + (7xx9)$
2. $63x$
3. $(7x) + (9x)$
Grade 3 Number Properties CCSS: 3.OA.B.5
Using the distributive property, what is 7(8 + 2) equal to?
1. 70
2. (7 x 8) + (7 x 2)
3. (8 x 7) + (8 x 2)
Grade 3 Number Properties CCSS: 3.OA.B.5
Grade 3 Number Properties CCSS: 3.OA.B.5
Using the distributive property, what is 2(9 + 7) equal to?
1. (2 x 9) + (2 x 7)
2. (9 x 2) + (9 x 7)
3. 32
Grade 3 Number Properties CCSS: 3.OA.B.5
$2 + (3 + 5) = (2 + 3) + 5$Which property applies?
1. Identity
2. Commutative
3. Associative
4. Distributive
$b * 4$ is equal to which of the following expressions?
1. $b + 4$
2. $(b + b + 4)$
3. $4 * b$
4. $b - 4$


text:
Rheological properties of halloysite soil slurry: a case study of weathered tephra involved in a shallow landslide triggered by the 2018 Eastern Iburi earthquake in Hokkaido, Japan – Earth, Planets and Space

May 11, 2022

Difference in yielding behaviors of slurries formed from dried and moist soil

The present experiments revealed that the yielding and fluidization behavior of slurries of weathered tephra depend greatly on several factors. As mentioned above, the yield stress is greatly affected by the slurry’s water content; despite this, comparison with the liquidity index shows that the yield stresses of both types of slurry are generally lower than those expected from empirical relations (Fig. 5b). One reason for this is that the yield stress for constructing the empirical curves was determined by a dynamic method, i.e., interpolation of a stress–strain rate curve to obtain a shear stress in the limit of zero shear rate, which is commonly higher than the yield stress determined by the static and direct method adopted here.

Creep tests by Carrière et al. (2018) on several clay soils at different IL states determined the critical shear stress at which the slurry viscosity bifurcates with increase of strain. Despite the clay slurries being prepared with dried samples, the stress values agreed well with the empirical relationship described above. Therefore, the low stress values obtained in this study may reflect the characteristics of the present sample, not the experimental method or sample preparation procedure.

An important finding of the present experiment is that the rheological properties of the slurry differ significantly between two different preparation methods (i.e., using dried or moist soil). Several studies have shown that the material properties of soils are altered by heating and drying (Casagrande 1932; Terzaghi et al. 1996; Huvaj and Uyeturk 2018). In general, after oven drying, the Atterberg limits such as liquid limit and plastic limit decrease by several tens of percent. Huvaj and Uyeturk (2018) studied the effect of drying on the Atterberg limit of pyroclastic soils containing halloysite, finding that the limits of samples dried at 60 °C or 110 °C were 1% to 30% lower than those of moist samples; they considered this likely owing to aggregation of the soil fabric or irreversible changes in the properties of halloysite after drying. Figure 5a shows a difference of approximately 0.03 to 0.04 in the scale of φ for similar yield stress values between the two slurries, corresponding to a difference of ~ 40% of the water content ratio. This is equivalent to a ~ 20% decrease in water content (as a proportion of the total moisture content), which is roughly consistent with the findings of Huvaj and Uyeturk (2018). Such differences led the estimated yield stress of the two slurries at the natural moisture state to vary by several orders of magnitude, from 2.2 kPa to several pascals. Kameda and Okamoto (2021) assumed Bingham flow in a numerical model of a landslide triggered by the Eastern Iburi earthquake. The modeling constrained the yield stress of the collapsed soils in a saturated state to be ~ 1.5 kPa, which is of the same order as that of the slurry containing moist soil (i.e., 2.2 kPa), suggesting that such slurry is more likely to represent the soil state when the landslide occurred. In fact, in situ monitoring of the tephra after the earthquake indicated sustained high saturation (> 90%) even during a period of no precipitation (Wang et al. 2021). Similar tephra containing halloysite is likely to be widely distributed in volcanic regions worldwide. If it is sufficiently dried to the extent of the dried soil used in this study and then saturated by rainfall, its yield stress may be reduced significantly, as demonstrated in Fig. 4. It is also expected that the soil strength in such a state will decrease with time rather than immediately after hydration; therefore, these cases will show a time difference between precipitation and the onset of coseismic ground motion due to an earthquake, which will be an important factor in determining slope stability during the earthquake.

Dependence of slurry properties on pH

Our results demonstrate a dependence of yield stress on slurry pH at a given water content. Although the two types of slurry showed differing sensitivity to pH, the yield stress generally decreased under alkaline conditions and increased under acidic conditions (Fig. 6a), as predicted previously from measurements of zeta potential for samples with a finer size fraction (Kameda 2021). Comparison of the two types of slurry showed slurry made of dried soil had a more moderate dependence on pH than that made of moist soil. Zeta potential measurements in the present study used samples of the same size fraction as those used for the rheological tests; the potential showed a near-steady increase from a large negative value (− 40 mV) to near-zero with decreasing pH (Fig. 6b). A large negative zeta potential can enhance colloidal stability, and thus reduce the shear strength of slurry (Plaza et al. 2018). More quantitatively, the following model of yield stress–DLVO (Derjaguin–Landau–Verwey–Overbeeck) force and interaction energy have often been applied to clay slurries based on a constant surface potential for interactions between spherical particles (Scales et al. 1998):

$${tau }_{y}approx frac{{phi }^{2}}{a}left(frac{A}{12{{D}_{0}}^{2}}-2pi {zeta }^{2}frac{{kappa mathrm{e}}^{-kappa {D}_{0}}}{1+{mathrm{e}}^{-kappa {D}_{0}}}right),$$

(6)

where a is the particle radius, A is the Hamaker constant, D0 is the surface separation distance between interacting particles, and (kappa) is the reciprocal of the Debye length (double-layer thickness). This model predicts a linear correlation between the yield stress and the square of zeta potential, and some clay slurries have been found to follow this correlation over a wide pH range (e.g., Au et al. 2014; Au and Leong 2016). As mentioned above, data for the slurry of dried soil show an almost linear correlation between yield stress and the square of zeta potential (Fig. 6c), indicating that the above yield stress–DLVO force model is likely applicable to the pH range tested here. In contrast, the results for the slurry of moist soil did not easily fit a single straight line. This may indicate either a change in the configuration of individual clay particles at different pH conditions (Shankar et al. 2010) or the addition of another attractive force, especially under acidic conditions, where the yield stress increased dramatically (Fig. 6a).

Although the microstructural mechanism for such behavior remains uncertain, the observed pH dependence of yield stress confirmed the prediction of Kameda (2021). Based on a simple mixing experiment of rainwater and tephra sample, Kameda (2021) raised the possibility of a change in pH of pore fluid in the tephra layer during rainfall, with pH likely increasing as the solid/liquid ratio decreases. If so, the chemical environment may have changed due to the rainfall before the earthquake in the present case, reducing the yield stress and thus the stability of the slope, which may have been a factor contributing to slope failure during the earthquake.

Fluidization of slurry under oscillatory strain

At the time of the Eastern Iburi earthquake, the slope was subjected to strong shaking with a maximum seismic intensity of 7 (JMA intensity scale) near the epicenter, which is thought to have directly triggered the slope collapse. To examine the effect of such periodic strain acting on the slope’s weathered tephra, dynamic viscoelasticity tests were conducted at frequencies of 0.5 to 10 Hz (Figs. 7 and 8). The results indicate that both types of initially solid-like slurry can liquify within the tested frequency range. However, the slurry of moist soil had a higher resistance to oscillatory strain than the other slurry, as expected from the vane test, and its fluidizing point occurred at a higher moisture content than in the case of the slurry of dried soil (Figs. 7d and 8d). In all cases, the strain amplitude of the solid–liquid transition (i.e., crossover strain) increased as the water content decreased or the frequency decreased. The slurry of moist soil, which is probably closer to the state of the soil at the time of the earthquake, had its crossover strain greatly increase when the applied frequency was less than 2 Hz, which may have been caused by the increase of storage modulus G′ above a shear strain of ~ 10% (Fig. 8a, b). In this study, experiments at even lower ratios of water content, such as the natural moisture state, could not be performed, because the samples cracked during shearing. However, extrapolation of the experimental data indicates that the phase transition is expected under conditions of natural water content at small strain conditions when the frequency exceeds 5 Hz (Fig. 8d).

Previous geotechnical experiments using the same tephra materials indicated that halloysite-bearing tephra is more susceptible to liquefaction than other horizons under cyclic loading at a frequency of 1.5 Hz (Li et al. 2020). On the other hand, ground motion with strain occurring in several pulses at higher frequency (~ 3 Hz or more) during the earthquake was recorded at several seismological stations near the epicenter (Wang et al. 2019). Recent continuous seismic observation after the Eastern Iburi earthquake revealed that the characteristic frequency of the slope oscillation associated with seismic activity ranged between 5 and 7 Hz (Wang et al. 2020). The experimental results indicate that such oscillatory strain could easily fluidize the tephra layer, while at lower frequencies it may exhibit some resistance to fluidization. Although this study tested sieved soil whose original texture was completely broken, the obtained rheological properties may be a key to understanding the factors contributing to the stability of the slope during the long period after tephra deposition and the mechanism of eventual collapse due to seismic shaking.


text:
# 3D homogenous coordinates versus 4D homogenous coordinates

We know that any 2D point $$(x, y)$$which represents as 3D homogeneous coordinates is of the form $$(x, y, 1)$$ which is the points of projective plane $$P^2.$$

If I use the same concepts for 3D points $$(x, y, z)$$ which represents as 4D homogeneous coordinates is of the form $$(x, y, z, 1).$$

My question is $$(x, y, z, 1)$$ could be the points of 3D projective space?

My second question is always homogenous coordinates are the points of projective space ?

N. B. -- $$P^2$$ is projective plane which contains the points of $$\mathbb{R^2}$$($$\mathbb{R^2}$$ points in $$P^2$$ can be represents as $$(x, y, 1)$$ ) and ideal points, is of the form $$(x, y, 0)$$ which is called points at infinity.

• Is $P^2$ supposed to be the real projective plane $RP^2$? Nov 22 at 0:26
• @lightxbulb explicitly $P^2=\mathbb{R^2}$+ ideal points. Nov 22 at 3:14
• Is the $+$ supposed to be union? Define "ideal points". Better yet - post the reference that you are using which employs this $P^2$ notation. Nov 22 at 6:32
• @lightxbulb ideal points is of form (x, y, 0) which is point at infinity. Nov 22 at 7:01
• The medium article explains your confusion. They define the term projective plane to be the plane $(x,y,1)$. It is true by their definition. This is usually not how the term is used as it refers to the real projective plane which is not just $(x,y,1)$. I could not find the definition of the medium article's author anywhere else (including the book that they cite). But by their definition $(x,y,z,1)$ would correspond to their definition of projective plane in 3D. I advise against using this terminology as I believe it is nonstandard and confusing (conflicts with the original term). Nov 22 at 10:00


text:
## CryptoDB

### Shabsi Walfish

#### Publications

Year
Venue
Title
2010
EUROCRYPT
2009
TCC
2008
CRYPTO
2007
TCC
2007
TCC
2006
TCC
2006
EPRINT
Cryptographic protocols are often designed and analyzed under some trusted setup assumptions, namely in settings where the participants have access to global information that is trusted to have some basic security properties. However, current modeling of security in the presence of such setup falls short of providing the expected security guarantees. A quintessential example of this phenomenon is the deniability concern: there exist natural protocols that meet the strongest known composable security notions, and are still vulnerable to bad interactions with rogue protocols that use the same setup. We extend the notion of universally composable (UC) security in a way that re-establishes its original intuitive guarantee even for protocols that use globally available setup. The new formulation prevents bad interactions even with adaptively chosen protocols that use the same setup. In particular, it guarantees deniability. While for protocols that use no setup the proposed requirements are the same as in traditional UC security, for protocols that use global setup the proposed requirements are significantly stronger. In fact, realizing Zero Knowledge or commitment becomes provably impossible, even in the Common Reference String model. Still, we propose reasonable alternative setup assumptions and protocols that allow realizing practically any cryptographic task under standard hardness assumptions even against adaptive corruptions.
2004
EPRINT
We present a new, elegant composition method for joint signature and encryption, also referred to as signcryption. The new method, which we call *Padding-based Parallel Signcryption* (PbPS), builds an efficient signcryption scheme from any family of trapdoor permutations, such as RSA. Each user U generates a single public/secret key pair f_U/f^{-1}_U used for both sending and receiving the data. To signcrypt a message m to a recipient with key f_{rcv}, a sender with key f_{snd} efficiently transforms m into a pair {w|s}, and simply sends { f_{rcv}(w) | f^{-1}_{snd}(s) }. PbPS enjoys many attractive properties: simplicity, efficiency, generality, parallelism of encrypting''/signing'', optimal exact security, flexible and ad-hoc key management, key reuse for sending/receiving data, optimally-low message expansion, long message and associated data support, and, finally, complete compatibility with the PKCS#1 infrastructure. The pairs {w|s} sufficient for the security of PbPS are called *universal two-padding schemes*. Using one round of the Feistel transform, we give a very general construction of such schemes. Interestingly, we notice that all popular padding schemes with message recovery used for plain signature or encryption, such as OAEP, OAEP+, PSS-R, and scramble all, encrypt small'', naturally consist of two pieces {w|s}. Quite remarkably, we show that all such pairs become special cases of our construction. As a result, we find a natural generalization of all conventional padding schemes, and show that any such padding can be used for signcryption with PbPS. However, none of such paddings gives optimal message bandwidth. For that purpose and of independent interest, we define a new hybrid'' between PSS-R and OAEP, which we call *Probabilistic Signature-Encryption Padding* (PSEP). We recommend using PbPS with PSEP to achieve the most flexible and secure signcryption scheme up-to-date. To justify this point, we provide a detailed practical comparison of PbPS/PSEP with other previously-proposed signcryption candidates.


text:
# Solving L1-Minimization problem with [R]NN

I want to solve $$\min \:\: \vert \vert x \vert \vert_1 \:\:\:\: s.t. \:\: Ax=y$$ using a Neural Network or Recurrent Neural Network.

I found only this paper, but there $x$ has to be greater or equal zero. Do you know any other good papers, tutorials or other projects about how to solve this problem using a network?

• This may not be the best way, but if you substitute $x1 - x2$ for $x$, where $x1$ and $x2$ are constrained to be non-negative, then you get a problem equivalent to your original, and involving only variables which are constrained to be non-negative. In particular, $x$ is now unconstrained. Jun 26 '17 at 18:41
• @whuber There's more than one way to skin a cat. In any event, a neural network is certainly not the most computationally efficient or robust way of solving a Linear Programming problem. So hopefully this whole exercise is just "for fun". Jun 26 '17 at 19:54
• Use software designed for Linear Programming, which implements a Simplex or Interior Point solver. You can not just pick up a book, any book, and hope to reproduce anything anywhere near that quality or robustness. BTW, if your problem has special structure, such as network structure, that can be exploited to increase the efficiency. Jun 26 '17 at 20:27
• @N8_Coder My point is that if you really want to solve a Linear programming problem, use software designed for it - don't use a neural network except for educational and experimental purposes. Jun 27 '17 at 19:22
• @N8_Coder You would be able to comment if you merge your accounts so that the system understood you were the asker rather than some random person. Do NOT post comments as answers. Jun 28 '17 at 6:56


text:
# Problems with a switch

## Recommended Posts

Carolina    176
Hi there! Well,I´m in trouble with college assigments again...I have to program something like a poor version of Excel,dealing with charts and cells containing values and formulae.Everything was going fine until I got to the point where I have to make a switch command,with various cases inside it.I have a char variable named Command that holds the value of any command the user inputs,resulting in chart modificatons. This variable named Command can be operators like +,-,* and /,and I´d like to make a case for when Command is one of these operators. I tried something like this: case '+': { code inside break; } But it didn´t work...However,case 'N',which is one of the possible cases too,works perfectly.what´s the secret here?=) Thanks in advance!

##### Share on other sites
GHertZ    122
i havent tried this out, but take a look at an ascii chart and put in the number for each sign that you want to use.

##### Share on other sites
There's no need to use an ASCII chart, that's just doing work the hard way that the compiler should do for you. When you use '+', it should subsitute the ASCII value of the + character automatically.

Truth be told, Carolina, you didn't give us enough information to help you with this one. switch / case works fine, and based on your wording there could be an infinite number of things going wrong for you.

The basic format should go something like this:

switch(character){case '+':    // code    break;case '-':    // code    break;case '/':    // code    break;default:    // code    break;}

If that's not working for you, it's simply something wrong in your code in a part you didn't describe here.

##### Share on other sites
NoMonkey    290
I suggest you put a break point in before your switch and check the value of Command each time a user presses the key. You may need to do what GHertZ suggested and look up the value you recieve in an ASCII table depending on how you are storing it. If Command is a char, then it should show you the actual letter otherwise, if it is an int, you will get the ASCII value.

Hope this helps!

##### Share on other sites
Carolina    176
Thanks for the help,I solved my problem... But now I have another: I am using a function named getche(),which is very useful because then I don´t need to press <enter> every time I input numbers or characters in my program.I use Dev-C++ compiler and I don´t even need to include libraries other than stdio.h and stdlib.h for it to work:it just does.
However,I need to submit the program to my professor and he uses gcc.Now I discovered that gcc doesn´t recognize getche()! Is there a way I can make gcc recognize it,like a library or something?Or is there a function that works similarly to getche()?

Thanks for everything,

Carol

##### Share on other sites
Zeikstraal    122
Have you tried adding #include <conio.h>?

##### Share on other sites
Carolina    176
Yeah,conio.h works in Dev-C++,but not in gcc...That´s why I´m thinking if maybe there isn´t some library in gcc that supports getche().I was searching at the Net and found there is one named ncurses tha supports getch(),but getch() and getche() aren´t the same...

##### Share on other sites
Zahlman    1682
Have you checked the behaviour of ncurses "getch"? If all it does is input one character and leave it on the screen, then you should be able to "fake it" by using the ncurses getch, then other ncurses functions to go back to the appropriate screen location and output a space. (Just outputting backspace followed by space might do it, too.)


text:
CBSE 10th Mathematics | Coordinate Geometry | Area of A Triangle

Coordinate Geometry | Area of A Triangle

The area of a DABC with vertices $\displaystyle A\left( {{x}_{1}},{{y}_{1}} \right),B\left( {{x}_{2}},{{y}_{2}} \right)$ and $\displaystyle C\left( {{x}_{3}},{{y}_{3}} \right)$ is given by area $\displaystyle \left( \Delta ABC \right)=\left| \frac{1}{2}\left\{ {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right) \right\} \right|.$

Proof :

Let $\displaystyle A\left( {{x}_{1}},{{y}_{1}} \right),B\left( {{x}_{2}},{{y}_{2}} \right)$ and $\displaystyle C\left( {{x}_{3}},{{y}_{3}} \right)$ be the vertices of the given DABC.

Draw AP, BM and CN perpendiculars to the x-axis.

Then, $\displaystyle MP=\left( {{x}_{1}}-{{x}_{2}} \right),PN\left( {{x}_{3}}-{{x}_{1}} \right)$

and MN = $\displaystyle \left( {{x}_{3}}-{{x}_{2}} \right)$

\     area of DABC

= ar (trap. BMPA) + ar (trap.. APNC) ar (trap. BMNC)

= $\displaystyle \left\{ \frac{1}{2}\left( AP+BM \right)\times MP \right\}+\left\{ \frac{1}{2}\left( AP+CN \right)\times PN \right\}-\left\{ \frac{1}{2}\left( BM+CN \right)\times MN \right\}$

= $\displaystyle \frac{1}{2}\left( {{y}_{1}}+{{y}_{2}} \right)\left( {{x}_{1}}-{{x}_{2}} \right)+\frac{1}{2}\left( {{y}_{1}}+{{y}_{3}} \right)\left( {{x}_{3}}-{{x}_{1}} \right)-\frac{1}{2}\left( {{y}_{2}}+{{y}_{3}} \right)\left( {{x}_{3}}-{{x}_{2}} \right)$

= $\displaystyle \frac{1}{2}\left[ {{x}_{1}}\left( {{y}_{1}}+{{y}_{2}}-{{y}_{1}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{2}}+{{y}_{3}}-{{y}_{1}}-{{y}_{2}} \right)+{{x}_{3}}\left( {{y}_{1}}+{{y}_{3}}-{{y}_{2}}-{{y}_{3}} \right) \right]$

$\displaystyle =\frac{1}{2}\left\{ {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right) \right\}.$

Since the area is never negative, we have

area $\displaystyle \left( \Delta ABC \right)=\frac{1}{2}\left| \left\{ {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right) \right\} \right|$

Condition for Collinearity of Three Points

Let the given points be $\displaystyle A\left( {{x}_{1}},{{y}_{1}} \right),B\left( {{x}_{2}},{{y}_{2}} \right)$ and $\displaystyle \left( {{x}_{3}},{{y}_{3}} \right)$.

Then, A, B and C are collinear

Þ area of DABC = 0

Þ  $\displaystyle \frac{1}{2}.\left\{ {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right) \right\}=0$

Þ $\displaystyle {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right)=0$

Solution:

The area of the triangle formed by the vertices A(1, –1), B(–4, 6) and C(–3, –5) is given by

$\displaystyle \frac{1}{2}[1(6+5)+(-4)(-5+1)+(-3)(-1-6)]$ = $\displaystyle \frac{1}{2}(11+16+21)=24$

So, the area of the triangle is 24 square units.

Solution:

Since the given points are collinear, the area of the triangle will be 0,

Þ    $\displaystyle \frac{1}{2}[2(k+3)+4(-3-3)+6(3-k)]=0$

Þ    $\displaystyle \frac{1}{2}[2k+6-24+18-6k]=0$

Þ     $\displaystyle \frac{1}{2}(-4k)=0$

Therefore, k = 0

Solution:

By joining B to D, two triangles ABD and BCD are formed.

Now the area of DABD

=

=$\displaystyle \frac{1}{2}(50+8+48)=\frac{106}{2}=53$ square units

Also, the area of DBCD

= $\displaystyle \frac{1}{2}[-4(-6-5)-1(5+5)+4(-5+6)]$

= $\displaystyle \frac{1}{2}(44-10+4)=19$ square units

So, the area of quadrilateral ABCD = 53 + 19 = 72 square units.

Solution:

Let D be the area of the triangle formed by the points (a, b + c), (b, c + a) and (c, a + b)

Þ    $\displaystyle \Delta =\frac{1}{2}\left| \left\{ a(c+a)+b(a+b)+c(b+c) \right\}-\left\{ b(b+c)+c(c+a)+a(a+b) \right\} \right|$

Þ    $\displaystyle \Delta =\frac{1}{2}\left| (ac+{{a}^{2}}+ab+{{b}^{2}}+bc+{{c}^{2}})-({{b}^{2}}+bc+{{c}^{2}}+ca+{{a}^{2}}+ab) \right|$

Þ    D = 0

Since area is zero \ the given points are collinear.

Solution:

The vertices of the given DABC are $\displaystyle A\left( k,2k \right),B\left( -2,6 \right)$ and $\displaystyle C\left( 3,1 \right)$.

\
$\displaystyle {{x}_{1}}=k,\,{{y}_{1}}=2k,{{x}_{2}}=-2,{{y}_{2}}=6$ and $\displaystyle {{x}_{3}}=3,{{y}_{3}}=1$

Area of $\displaystyle \Delta ABC=\frac{1}{2}\left| \left\{ {{x}_{1}}\left( {{y}_{2}}-{{y}_{3}} \right)+{{x}_{2}}\left( {{y}_{3}}-{{y}_{1}} \right)+{{x}_{3}}\left( {{y}_{1}}-{{y}_{2}} \right) \right\} \right|$

=$\displaystyle \frac{1}{2}\left| k\left( 6-1 \right)-2\left( 1-2k \right)+3\left( 2k-6 \right) \right|$

=$\displaystyle \frac{1}{2}\left| 5k-2+4k+6k-18 \right|$

= $\displaystyle \frac{1}{2}\left| 15k-20 \right|$ sq units.

But, area of DABC = 5 sq units     [given]

\    $\displaystyle \frac{1}{2}\left| 15k-20 \right|=5\,\,\,\,\,\,\Rightarrow \,\,\,\,\,\left| 15k-20 \right|=10$

Þ
$\displaystyle \left( 15k-20 \right)=10$ or $\displaystyle \left( 15k-20 \right)=-10$

Þ
$\displaystyle k=2$ or $\displaystyle k=\frac{2}{3}$.

Hence,          $\displaystyle k=2$ or $\displaystyle k=\frac{2}{3}$

This site uses Akismet to reduce spam. Learn how your comment data is processed.


text:
• Subject:

...

• Topic:

...

The earth's magnetic field at a certain place has a horizontal component 0.3 Gauss and the total strength 0.5 Gauss. The angle of dip is

(a) ${\mathrm{tan}}^{-1}\left(\frac{3}{4}\right)$                      (b) ${\mathrm{sin}}^{-1}\left(\frac{3}{4}\right)$

(c) ${\mathrm{tan}}^{-1}\left(\frac{4}{3}\right)$                      (d) ${\mathrm{sin}}^{-1}\left(\frac{3}{5}\right)$


text:
# Newton’s method I: Quadratic convergence rate

The next couple posts will focus on our favorite second order method: Newton’s method. I’ve been going through them, partly as “review” (in quotes because I ended up learning a lot of new things) and partly to develop some intuition as to when acceleration-by-leveraging-second-ordery-info might actually help. This first post will be super short, and focus on analysis of Newton’s method under two strong-convexity-like conditions.

Consider the unconstrained minimization of $f(x)$ where $f$ is everywhere twice differentiable. Then Newton’s method iterates as $$x^{(t+1)} = x^{(t)} – (\nabla^2 f(x^{(t)}))^{-1}\nabla f(x^{(t)}).$$ Unlike gradient descent, there is no “step size”; the curvature information is already incorporated in the Hessian inverse. Also, as a responsible optimist, I must point out that in practice we never take the Hessian inverse, but instead compute the matrix-inverse-vector-product using some iterative method, like conjugate gradient.

Suppose the Hessian is $M$-smooth $$\|\nabla^2 f(x)-\nabla^2 f(y)\|_2 \leq M\|x-y\|_2.\quad (\star)$$ We all learn Newton’s method as having a locally quadratic convergence rate. The usual question is then: what is the global convergence rate? Well, actually, if you stare and squint at the proofs given in Nesterov’s works, you see the answer is actually right in front of you.

• If the function is convex but not strongly convex, then Newton’s method may not work at all. In fact, for any point $x^{(t)}$, if $\nabla^2 f(x^{(t)})$ has 0 eigenvalues, you can either take a damped Newton step (via line search) or switch to gradient steps. So, unsurprisingly, the rate here should be $O(1/t)$, the same as in gradient descent.
• If the function is $\mu$-strongly convex and $\mu > 0$, then Newton’s method has a locally quadratic convergence rate, and the proof is like 3 lines.
• If the function is $\mu$-strongly convex at the optimum, then we still cannot say anything globally, but we can say that the method has local quadratic convergence. The proof is a bit longer, and the point at which quadratic convergence kicks in could be much slower, though.

In other words, it’s an all-or-nothing game. Either I get quadratic convergence, or I am better off switching to gradient descent (at least until we’re close enough to the quadratic regime).

The proof for this convergence is actually really cute, once you know all the pieces. Use the integral form for $\nabla f(x^{(t)})$: $$\nabla f(x^{(t)}) -\nabla f(x^*)= \int_0^1 \nabla^2 f(x^* – \tau (x^{(t)}-x^*) ) (x^{(t)}-x^*) d\tau.$$ Then, since $\nabla f(x^*) = 0$,

\begin{eqnarray*}
x^{(t+1)}-x^* &=& x^{(t)} -( \nabla^2 f(x^{(t)}))^{-1} \nabla f(x^{(t)})\\
&=& x^{(t)} – ( \nabla^2 f(x^{(t)}))^{-1} \left( \int_0^1 \nabla^2 f(x^* – \tau (x^{(t)}-x^*) ) (x^{(t)}-x^*) d\tau\right)\\
&=&( \nabla^2 f(x^{(t)}))^{-1}\left(\int_0^1\nabla^2 f(x^{(t)}) – \nabla^2 f(x^* – \tau (x^{(t)}-x^*) d\tau \right)(x^{(t)}-x^*).
\end{eqnarray*}

Taking norms,

\begin{eqnarray*}
\|x^{(t+1)}-x^*\|_2 &\leq & \| \nabla^2 f(x^{(t)}))^{-1}\|_2\left(\int_0^1 M\|x^*-\tau (x^{(t)}-x^*)-x^*\|_2d\tau\right)\|x^{(t)}-x^*\|_2 \\
&=& \| \nabla^2 f(x^{(t)}))^{-1}\|_2\;\underbrace{(\int_0^1 \tau d\tau)}_{=1/2}\; M \|x^{(t)}-x^*\|_2^2 \\
&=& \frac{M}{2}\| \nabla^2 f(x^{(t)}))^{-1}\|_2 \|x^{(t)}-x^*\|_2^2.
\end{eqnarray*}

Now the question is, what is an upper bound on $\| \nabla^2 f(x^{(t)}))^{-1}\|_2$? Well, if $f$ is $\mu$-strongly convex everywhere, then $\| \nabla^2 f(x^{(t)}))^{-1}\|_2 \leq 1/\mu$ everywhere, and we can say globally that

$$\|x^{(t+1)}-x^*\|_2 \leq \frac{M}{2\mu} \|x^{(t)}-x^*\|_2^2.$$

Now, if $\| \|x^{(t)}-x^*\|_2 < \frac{2\mu}{M}$, this implies contraction, and henceforth, for $t’ > t$, there will be quadratic convergence. If that hasn’t yet occurred, a simple line search add-on should make sure that a damped Newton’s step makes progress, globally.

If, on the other hand, we loosen our requirements to simply that $\lambda_{\min}(\nabla^2 f(x^*)) \geq \mu$, then by ($\star$), $\lambda_{\min}(\nabla^2 f(x^{(t)}) \geq \mu – M\|x^{(t)}-x^*\|_2$, and

$$\|x^{(t+1)}-x^*\|_2 \leq \frac{M}{2(\mu-M\|x^{(t)}-x^*\|_2^2)} \|x^{(t)}-x^*\|_2^2\leq \frac{M}{2\mu(1-\theta))} \|x^{(t)}-x^*\|_2^2$$

if $\|x^{(t)}-x^*\|_2 \leq \frac{\mu\theta}{M}$ for some $\theta < 1$. Nesterov picks $\theta = 2/3$, but anything will do.

I don’t think Nesterov himself was very satisfied with this rate; in particular, the second case where we’re adding a new parameter we can’t measure ($\theta$) on top of the $L$-continuous Hessian and $\mu$-strongly convex, it is hard to list exactly which functions are satisfying these conditions. It seems like this was the motivation behind the next two topics: Newton with cubic regularization for generalized functions, and Newton’s method for self-concordant functions.

## Reference

[1] Nesterov, Y. Lectures on Convex Optimization. 2018, Chapter 1.2.

Edit July 12 2002: Thanks Mark Schmidt for pointing out that there is a serious flaw in my statement that Newton’s method can converge quadratically, globally. We can guarantee globally that $\|x^{(t+1)}-x^*\|_2 \leq \frac{M}{2\mu}\| \nabla^2 f(x^{(t)}))^{-1}\|_2 \|x^{(t)}-x^*\|_2^2$, but unless $\|x^{(t)}-x^*\|_2 < \frac{2\mu}{M}$, this does not imply contraction. In other words, even if $f$ has an $L$-continuous Hessian and is $\mu$-strongly convex, we still need to use line search and damped Newton in the initial phases until the locally quadratic convergence rate kicks in. tl;dr, no global quadratic convergence (sorry!). This has been amended in the post.


text:
COSAM » Events » 2016 » October » Applied Mathematics Seminar

 Applied Mathematics Seminar Time: Oct 28, 2016 (02:00 PM) Location: Parker Hall 328 Details: Speaker: Hans-Werner van Wyk Title: A novel finite difference method for the fractional Laplacian Abstract: In this talk we discuss a novel finite difference method to discretize the fractional Laplacian operator over a bounded domain. Our approach is based on splitting the hypersingular kernel function in its integral representation into two parts, allowing some of the singularity to be absorbed by the operator’s input function $$u$$. We then approximate the resulting weakly singular integral by a weighted trapezoidal rule. Our method is simple to implement and has a higher accuracy than the existing finite difference approximations. For smooth enough $$u$$ and a judicious choice of the splitting parameter, we obtain second order accuracy, uniformly for any diffusion order.

Last updated: 10/27/2016


text:
RevEcoR – Reverse ecology linking organism(s) and their environments

Background

Ecology is the scientific With the rapid and inexpensive next-generation sequencing technologies The structure of complex biological systems reflects not only their function but also into the habitates in which they envoled and are adapted to. Reverse Ecology- an emerging new frontier in Evolutionary Sysmtems biology at the interface of computational biology, genomics and environmental science, which uses population genomics to study ecology with no a priori assumptions about the organism(s) under consideration. This term was suggested by Li et al. during a conference on ecological genomics in Christchurch, New Zealand. It facilitates the translation of high through genomic data into large scale ecological data, and utilizes system-based method to obtain novel insights of pooly characterized microorganisms and relationships between microorgasnims or their environments in a superorganism.Traditional approach, however, can only applied to a small scale and for ralatively well-studied systems.

Introduction

This manual is a brief introduction to structure, funcitons and usage of RevEcoR package. The RevEcoR package implements the reverse ecology framework. It can be used for reconstruction of metabolic environmnets using a cross-species analysis, idenditifying the set of compounds and predicting the species interactions on a large scale with a graph-based algorithm mentioned in Borenstein et al.. The reverse ecology framework which takes advances of system biology and genomic metabolic modeling, aims to predict the ecological traits of poorly studied microorganisms, their interactions with other microorganisms, and the ecology of microbial communities from system-level analysis of complex biological networks.

Two softwares (NetSeed and NetCooperate), have developed by Borenstein’s lab for studying the ecology of microorganisms. Howevever, neither of them supports metabolic network reconstruction of species, and both are limited in small scale analysis (two species). The RevEcoR package provides an interface to microbiome reverse ecology analysis on a large scale. The main features of this package consists of several steps:

1. Downding metabolic networks data from KEGG database: all species metabolic network data could be downloaded from the KEGG PATHWAY datase with the KEGG REST API (application programming interface).

2. Reconstruction metabolic networks of all species: a directed graph whose nodes represent compounds and whose edges represent reactions linking substrates to products.

3. Identify seed set of a specific organism: each metabolic network was decomposed into its strongly conneted components (SCC) using Kasaraju’s algorithm. The SCC forms a directed graph whose nodes are the components and whose edges are the orginal edges in the graph that connect nodes in two different components. Then, detectting the seed set with this SCCs.

4. Host-microbe and microbe-microbe cooperation analysis

Installation

RevEcoR is free available on CRAN. You can install the latest released version from CRAN as following:

install.packages("RevEcoR")

or the latest development version from github. To install packages from GitHub, you first need install the devtools package on your system with install.packages("devtools"). Note that devtools sometimes needs some extra non-R software on your system – more specifically, an Rtools download for Windows or Xcode for OS X. There’s more information about devtools here.

if (!require(devtools)
install.packages("devtools")
devtools::install_github("yiluheihei/RevEcoR")

After installation, you can load RevEcoR into current workspace by typing or pasting the following codes:

library(RevEcoR)

Downloding the metabolic data of a specific organism

Both the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Integrated Microbial Genomes database (IMG) [7] collect complete high-quality genome sequences and metagenome sequences with a comprehensive set of publicly available bacterial, archaeal, eukaryotic, and phage genomes, as well as engineered, environmental and host associated metagenome samples. All the sequences and functional annotation profile can be download directly.

The KEGG database provides a REST-style API compared with IMG database. This package provides function getOrgMetabolicData to download the specific organism metabolic data from KEGG database, and return a list where each element consists of three elements: reaction, substrate and product. Data in IMG, we can only download it manually. In addition, users can annotate their private genomic data with KO terms on IMG systems or in their local machine to obtain the annotation profile.

## download sample metabolic data from remote KEGG database
buc <- getOrgMetabolicData("buc")
data(kegg_buc)
head(buc)

### Reconstruction of organism metabolic network

The crucial premise of Reverse Ecology is that genomic information can be converted into ecological information. Metabolism linking miroorgasnim with the biochemical environmnet acrossing compound exchanges: import of exogenous compounds and impact the composition of its environment via secretion of other compounds. Thus, metabolic network is the reflection of interation between organisms and their environment.

Graph-based representation of metabolic reactions where nodes represent compounds and edges represent reactions is a common tool in analyzing and studing metabolic networks. A directed edge from node a to b indicates that compound a is a substrate in some reaction that produces compound b.

Once the metabolic data is obtained, reconstructGsMN could be used to reconstruction the metabolic network of a specific organism.

## species in KEGG
buc.net <- reconstructGsMN(kegg_buc, RefData = NULL)
igraph::print.igraph(buc.net)
## IGRAPH DN-- 259 681 --
## + attr: name (v/c)
## + edges (vertex names):
##  [1] cpd:C00022->cpd:C00002 cpd:C00022->cpd:C00008 cpd:C00022->cpd:C00074
##  [4] cpd:C00022->cpd:C00024 cpd:C00022->cpd:C06010 cpd:C00022->cpd:C00109
##  [7] cpd:C00022->cpd:C00118 cpd:C00022->cpd:C00900 cpd:C00022->cpd:C05125
## [10] cpd:C00022->cpd:C06006 cpd:C00022->cpd:C11437 cpd:C00068->cpd:C01081
## [13] cpd:C00068->cpd:C05125 cpd:C00068->cpd:C05381 cpd:C04332->cpd:C00255
## [16] cpd:C04332->cpd:C04732 cpd:C04332->cpd:C15556 cpd:C00255->cpd:C04332
## [19] cpd:C00255->cpd:C00061 cpd:C00255->cpd:C04732 cpd:C00003->cpd:C00006
## [22] cpd:C00003->cpd:C00455 cpd:C00003->cpd:C00857 cpd:C00006->cpd:C00003
## + ... omitted several edges
igraph::plot.igraph(buc.net, vertex.label=NA, vertex.size=5, edge.arrow.size=0.1)

## ko annotation profile species detected in a human microbiome in IMG (not in KEGG)
annodir <- system.file("extdata/koanno.tab",package = "RevEcoR")
data(RefDbcache)
g2 <- reconstructGsMN(metabolic.data, RefData = RefDbcache)

Identify seed set of metabolic network

As the interactions with the environment was reflected in the metabolic networks, these networks could be used not only to infer metabolic function but alse to obtain insights into the growth environments in which the species evolved.

Apparantly, organisms can survive in a wide range of enviromnets and may activate only a subset of the pathways in the network of each environment, using a different set of exogenously acquired compounds (termed seed set). The seed set of the network is defined as the minimal set of compounds in the network that allows the synthesis of all other compouds, and can serve as a good proxy for its environment and can be conceived as the essential and effective biochemical environment.

getSeedSets was used to detect the seed set of a metabolic network which returns a seedset-class object. Futhermore, some methods was supported for seedset-class, e.g length, nonseeds. For more details on seedset-class, see ?seedset-class.It can help us to get the compound that organisms are exogenously acquired compouds from the environment, representing the organism’s nutritional profile. This algorithm is based on a fast method Kasaraju algorithm for SCC decomposition which is implementaed in Kasaraju. For more details, see ?getSeedSets and ?KasarajuSCC.

## seed set prediction
seed.set <- getSeedSets(buc.net, 0.2)
show(seed.set)
## Object of class seedset
##    IGRAPH: -- 259 nodes 681 edges --
##   seedset length  9
head(seed.set@seeds)
## [[1]]
## [1] "cpd:C15532"
##
## [[2]]
## [1] "cpd:C15811"
##
## [[3]]
## [1] "cpd:C05729"
##
## [[4]]
## [1] "cpd:C00993"
##
## [[5]]
## [1] "cpd:C00031"
##
## [[6]]
## [1] "cpd:C00392"
## The node colored with red represents the species' seed set
nodes  <- igraph::V(buc.net)$name seeds <- unlist(seed.set@seeds) seed.index <- match(seeds,nodes) node.color <- rep("SkyBlue2",length(nodes)) node.color[seed.index] <- "red" igraph::plot.igraph(buc.net, vertex.label=NA, vertex.size=5, edge.arrow.size=0.1, vertex.color = node.color) Predict species interactions The topology of metabolic networks can provide insight not only into the metabolic process that accur within each species, but also into interactions between different species. Here we provides calculateCoopreationIndex using three cooperation index: competition index, coplementarity index measure the microbe-microbe co-occurrence pattern. More details, see ?calculateCoopreationIndex. # ptr metabolic network data(kegg_ptr) ##ptr.net <- reconstructGsMN(getOrgMetabolicData("ptr")) ptr.net <- reconstructGsMN(kegg_ptr) # cooperation analysis between buc and ptr cooperation.index <- calculateCooperationIndex(buc.net,ptr.net) cooperation.index  ##$competition.index
##            [,1]      [,2]
## [1,] 1.00000000 0.1818182
## [2,] 0.03571429 1.0000000
##
## $complementarity.index ## [,1] [,2] ## [1,] 0.00000000 0.4545455 ## [2,] 0.07142857 0.0000000 Blow we applied RevEcoR to predict interactions among several human oral microbiota species whose co-occurrence patterns have previously been described. Our seven sample oral species metabolic dataset was downloaded from the IMG server; it consists of the following Genomes Online Database (GOD) IMG identification numbers: Gc0016386, Gi07614, Gi00264, Gc00809, Gc00643, Gi03876, and Gi07289. The interactions among the seven species is calculated as follwing: The metabolic networks of seven species was reconstructed first: ##metabolic network reconstruction of these seven species net <- lapply(anno.species, reconstructGsMN) Then caculateCooperationIndex is used to calculated the species interactions among the seven species. ## caculate interactions among vious species interactions <- calculateCooperationIndex(net, p = TRUE) ## competition index$competition.index
Aa        Ao        Fn        Pg        Sg        So        Va
Aa 1.0000000 0.4736842 0.3157895 0.2280702 0.4210526 0.4385965 0.2456140
Ao 0.4736842 1.0000000 0.3684211 0.3333333 0.4736842 0.4736842 0.2456140
Fn 0.5000000 0.5833333 1.0000000 0.4166667 0.5833333 0.5555556 0.4166667
Pg 0.4193548 0.6129032 0.4838710 1.0000000 0.6129032 0.5161290 0.3870968
Sg 0.5454545 0.6136364 0.4772727 0.4318182 1.0000000 0.9090909 0.3863636
So 0.5813953 0.6046512 0.4651163 0.3720930 0.9302326 1.0000000 0.3953488
Va 0.4827586 0.4827586 0.5172414 0.4137931 0.5862069 0.5862069 1.0000000
## p value of competition index
$competition.index.p Aa Ao Fn Pg Sg So Va Aa 0.000 0.001 0.001 0.001 0.001 0.001 0.001 Ao 0.001 0.000 0.001 0.001 0.001 0.001 0.001 Fn 0.001 0.001 0.000 0.001 0.001 0.001 0.001 Pg 0.001 0.001 0.001 0.000 0.001 0.001 0.001 Sg 0.001 0.001 0.001 0.001 0.000 0.001 0.001 So 0.001 0.001 0.001 0.001 0.001 0.000 0.001 Va 0.001 0.001 0.001 0.001 0.001 0.001 0.000 ## complementarity index$complementarity.index
Aa        Ao        Fn         Pg        Sg         So        Va
Aa 0.0000000 0.1052632 0.1228070 0.07017544 0.0877193 0.08771930 0.1228070
Ao 0.1403509 0.0000000 0.1403509 0.07017544 0.1228070 0.12280702 0.1403509
Fn 0.1944444 0.1666667 0.0000000 0.16666667 0.1111111 0.11111111 0.1388889
Pg 0.2258065 0.2258065 0.1612903 0.00000000 0.1612903 0.19354839 0.2258065
Sg 0.2272727 0.1818182 0.1590909 0.09090909 0.0000000 0.04545455 0.1590909
So 0.1860465 0.1395349 0.1860465 0.09302326 0.0000000 0.00000000 0.1395349
Va 0.2068966 0.1724138 0.1379310 0.17241379 0.1379310 0.13793103 0.0000000
## p value of complementarity index
$complementarity.index.p Aa Ao Fn Pg Sg So Va Aa 0.000 0.001 0.001 0.001 0.001 0.001 0.001 Ao 0.001 0.000 0.001 0.001 0.001 0.001 0.001 Fn 0.001 0.001 0.000 0.001 0.001 0.001 0.001 Pg 0.001 0.001 0.001 0.000 0.001 0.001 0.001 Sg 0.001 0.001 0.001 0.001 0.000 0.001 0.001 So 0.001 0.001 0.001 0.001 0.001 0.000 0.001 Va 0.001 0.001 0.001 0.001 0.001 0.001 0.000 We found that Streptococcus oralis (So) and Streptococcus gordonii (Sg) had the lowest complementarity index (0.04 and 0.00, respectively) and the highest competition index (0.91 and 0.93, respectively) among all pairs. Appanrently, the p value of all interactions are biologcial/statistical significance. This indicates that these two species are antagonistic, which is consistent with previous findings. Comparing predicted interactions and co-occurrences To further evaluate the predicted interactions of RevEcoR on a large scale, we integrated functions mentioned above to investigate species interactions in the gut microbiome. We focused on a list of 116 prevalent gut species, whose genome sequence is available in IMG database and sequence coverage is more than 1% in at least one metagenomic sample of 124 individuals. Genome annotation profiles of this 116 species was collected from IMG database and was used to calculated the interactions (competition and complementarity index) for all pairs of species. 1. Predicting species interactions in gut microbiome For each species, we download the list of genes mapped to the Kyoto Encyclopedia of Genes and Genomes orthologous groups (KOs) was downloaded with a in-house R script. This data, which was used to reconstruct the metabolic network of each species, have been saved as gut_microbiome.rda in subdirectory data of RevEcoR. We can load it as the following code: data(gut_microbiome) ## summary(gut_microbiome)  Then, it was used to reconstuct the metabolic network, predict the seed set, and finally predict the pairs of interactions between different species: gut.nets <- lapply(gut_microbiome,reconstructGsMN) seed.sets <- lapply(gut.nets,getSeedSets) ## Since calculation is on large scale, species interactions prediction may take several hours gut.interactions <- calculateCooperationIndex(gut.nets) competition.index <- gut.interactions$competition.index
complementarity.index <- gut.interactions\$complementarity.index

2. Obtaining co-occurrence scores in gut microbiome

Specially, it will help us to predict whether species compete with one another tend to co-occur or to exclude by comparing the species interactions and co-occurrences. We obtained co-occurrence scores directly from Levy R’s research, and saved it as occurence.tab in subdirectory inst/extdata of RevEcoR . Co-occurrence score was calculated based on species abundances across all samples and measured by the Jaccard similarity index.

occurrence.score <- read.delim(system.file("extdata/occurrence.tab",
package = "RevEcoR"),stringsAsFactors = FALSE, quote = "")

3. Comparing the interactions and co-occurrences

The Jaccard index measures similarity between finite sample sets, and is defined as the size of the intersection divided by the size of the union of the sample sets. Thus, co-occurrence scores are generally symmetric whereas competition index or complementarity index are not according to their definitions.

A symmetric version of two interaction index was gnerated by replacing each element of the interaction indices with the mean of each value and that transpose value.

competition.index <- (competition.index + t(competition.index))/2
complementarity.index <- (complementarity.index + t(complementarity.index))/2

Subsequently, the Spearman correlation between the co-occurrence scores and the two interaction indices was calculated. A permutation-based Mantel test, which is commonly used in ecology, was used to determin the significance of this correlation.

## upper triangles, which is used to calculate the correlation
competition.upper <- competition.index[upper.tri(competition.index)]
occurrence.upper <- occurrence.score[upper.tri(occurrence.score)]
complementarity.upper <- complementarity.index[upper.tri(complementarity.index)]

## calculate the spearman correlation betwwen co-occurrence scores and two
## interactions indices
competition.cor <- cor(competition.upper,occurrence.upper,method="spearman")
complementarity.cor <- cor(complementarity.upper,occurrence.upper,method="spearman")

## permutation-based mantel test. Random permutation the co-occurance score
## 10000 times, P value is the fraction of correlations as high as or higher
## than the original
if (require(magrittr)){
null.stat <- replicate(10000,
sample(1:116) %>% occurrence.score[.,.] %>%
.[upper.tri(.)]
)
competition.null <- cor(competition.upper,null.stat)
complementarity.null <- cor(complementarity.upper,null.stat)
length(which(competition.null >= competition.cor)) ## 0 p.competition < 0.00001
length(which(complementarity.null <= complementarity.cor)) ## 0 p.complementarity< 0.00001
}

We found that competition index is significant positively correlated with co-occurrence (cor = 0.261, P < 10-4, Mantel correlation test), whereas the complementarity index is significant negatively correlated with co-occurrence (cor = -0.259, P < 10-4, Mantel correlation test). This suggests that competition is liked to be the key factor to promote the assembly of gut microorganisms

sessionInfo

The version number of R and packages loaded for generating the vignette were:

sessionInfo()
## R version 3.2.4 Revised (2016-03-16 r70336)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 10240)
##
## locale:
## [1] LC_COLLATE=C
## [2] LC_CTYPE=Chinese (Simplified)_China.936
## [3] LC_MONETARY=Chinese (Simplified)_China.936
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.936
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
## [1] RevEcoR_0.99.3 knitr_1.12.3
##
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.4     lattice_0.20-33 XML_3.98-1.4    gtools_3.5.0
##  [5] digest_0.6.9    plyr_1.8.3      grid_3.2.4      formatR_1.3
##  [9] magrittr_1.5    evaluate_0.8.3  stringi_1.0-1   Matrix_1.2-4
## [13] rmarkdown_0.9.5 tools_3.2.4     stringr_1.0.0   purrr_0.2.1
## [17] igraph_1.0.1    yaml_2.1.13     htmltools_0.3.5